Probabilistic, Fractal, and Related Techniques for Analysis of Engineering Surfaces

In many engineering fields surface topography is of crucial importance solving problems of friction and other problems of tribology. A review of mathematical approaches for description of topography of engineering surfaces is presented. Firstly, we give a brief introduction to some of statistical parameters used for description of surface roughness. It is argued that although some of these parameters may be quite useful for specific engineering problems, a set of finite numbers of parameters cannot describe contact properties of rough surfaces. Then we discuss various models of surface roughness based on Gaussian models of the asperity heights. The results of application of various modern tests of normality for checking whether the distribution of the asperity heights is Gaussian, are presented. Further fractal models of roughness are discussed. Using fractal parametric-homogeneous (PH) surfaces, it is demonstrated that tribological properties of a rough surface cannot be characterized just by the fractal dimension of the surface. It is also shown that models based solely on the power-spectral density function (PSDF) are quite similar to fractal models and these models do not reflect tribological properties of surfaces. In particular, it is demonstrated that different profiles may have the same PSDF.


INTRODUCTION
The paper deals mainly with surfaces used in engineering practice that will be referred to as engineering surfaces. It is known that all engineering surfaces are rough (see e.g., Archard et al., 1975;Whitehouse, 2011) and therefore, contact between engineering surfaces is realized by a number of contact spots (see e.g., Zhuravlev, 1940Zhuravlev, , 2007Holm, 1941;Goryacheva, 1998;Borodich, 2007). If the surface profile z(x) is studied using Fourier decomposition, and the term 'roughness' is attributed to the short wavelength shapes, while the long wavelength shapes are referred to as "waviness" of the surface (see e.g., Morales-Espejel et al., 2000;Borodich and Bianchi, 2013). If the waviness is extracted from the surface profile then the rough surface may be considered as nominally flat (see e.g., Greenwood and Williamson, 1966). Roughness of engineering surfaces is a crucial factor for performance of tribological components. The energy dissipation during sliding of dry engineering surfaces and correspondingly, the friction are enormously influenced by the surface profile (see e.g., Borodich and Savencu, 2017). Here we present a critical review of some popular statistical, fractal and related techniques for modeling and analysis of the surface roughness.
One of the first attempts to employ statistical methods for description of surface roughness was presented by Abbott and Firestone (1933) who calculated the cumulative distribution function of the surface heights. In tribology this parameter is called the Abbott-Firestone curve or the bearing area curve. Independently, Zhuravlev (1940) employed this parameter in his statistical model of contact between rough surfaces that were represented as collection of spherical protuberances having identical radii. He explained that the number of contacting spheres a specific height increases as the surfaces approach each other. In the English language literature this model is usually attributed to Greenwood and Williamson (1966).
After the bearing area curve parameter was introduced, there was a period that can be referred to as "the parameter rash" because a huge number of statistical parameters of roughness were introduced (Whitehouse, 1982). These characteristics were related to both the vertical distribution of heights and the horizontal distribution of the rough profiles (Nowicki, 1985). The next step in surface roughness characterization was the idea that it may be modeled using theory of random processes. In 1953 Linnik and Khusu presented a seminar talk where they suggested to use graphs of a stationary Gaussian random process in order to describe surface roughness, see Linnik and Khusu (1954a) for detail, as well as Linnik and Khusu (1954b) and Khusu et al. (1975). Linnik and Khusu (1954a) suggested to study the following correlation function for the Gaussian random process where K(0) and α are some parameters of the roughness. Independently, the same idea was introduced later by Whitehouse and Archard (1970). They presented an absolutely correct statement that if a profile z(x) of a random rough surface is Gaussian then it can be fully described by a distribution of asperity heights and the correlation (auto-correlation) function of the process R. The auto-correlation is defined as If one takes the Fourier transform of R(δ) then the power spectrum G(ω) or the power-spectral density function (PSDF) is obtained. If the signal frequency is denoted as ω then the PSDF is defined as Developing the random signal approach, Sayles and Thomas (1978) presented experimental relations between wavelength and the scaled power spectral density for many different surfaces. They argued that the scaled spectral density functions of many surface profiles can be approximately presented as G(ω) = 2π /ω 2 . Sayles and Thomas (1978) referred to as the surface topothesy. As Dr. Sayles said to one of the authors (FB), they never claimed that the real surfaces are fractal; in fact the fractal terminology to surface roughness description was triggered by Berry and Hannay (1978) who presented a comment to Sayles and Thomas (1978) paper where they claimed that geometric properties of rough surfaces can be characterized by a new concept 'fractal' that was described in detail by Mandelbrot (1977). Another important step in the promotion of the fractal approach to surface roughness description was the studies of the Weierstrass-Mandelbrot fractal function by Berry and Lewis (1980). Later the Weierstrass type functions were used by many researchers as a model of rough surfaces (see e.g., Roques-Carmes et al., 1988;Majumdar and Bhushan, 1990). For some period of time, the fractal models became very popular, there were even statements that "fractals are everywhere." Fractal approach to surface topography were so popular that one could say that it became an "emperor" of many research areas. Speaking about fractal approaches in fracture mechanics Borodich (1999) argued that instead of careful presentations of the state-of-theart, the papers dedicated to fractal analysis are often based on repetition of common myths about fractals (we call this as vulgar fractal approaches). Hence, in most of the papers dedicated to fractal approaches to fracture and surface topography, the state-of-the-street is ruling. Borodich (2002) listed examples introduced earlier by Borodich and Onishchenko (1993) and Borodich (1993) and reminded that "the fractal dimensions alone cannot characterize the features of contact." The same situation is related to the papers dedicated to fractal approaches to surface topography, nevertheless there is still a stream of papers based on vulgar interpretation of the fractal approaches. One could mention here a statement by Mandelbrot (1998) that "fractals are not a panacea; they are not everywhere." Using examples introduced by Borodich and his co-authors, we will show that quite often there is no meaning in fractal analysis of surface roughness.
Nowadays another tendency is quite popular, namely to describe rough surfaces using solely the PSDFs of surface topography. We argue that these papers are in essence an attempt to resurrect the fractal approach. Indeed, these papers contain usually a mixture of correct statements related to Gaussian processes and wrong statements based on attempts to extend the power-spectral analysis to non-Gaussian surfaces. In addition, these papers suffer often by employment of the ill-defined terms, such as the Hurst exponent, non-accurate statements about selfaffinity of surface roughness and vulgar interpretation of fractal models. We will show that the power-spectral analysis applied to non-Gaussian surfaces is a kind of reformulation of the vulgar fractal approach. Using an analogy to Andersen's tail about new clothes produced by cunning weavers, we can say that attempts to model surfaces solely by the use of the PSDF of its roughness are "the emperor's new clothes." The paper is organized as follows. In section 2 probabilistic characteristics of rough surfaces are discussed. In section 3 we consider some approaches to modeling of surface roughness using graphs of random processes that in turn, assume that the asperity heights are normally distributed or they employ similar assumptions that involve normal distributions. We give a brief description of statistical methods employed for checking normality of distributions and some results of application of these tests to roughness of engineering surfaces. In section 4 we discuss briefly the fractal approaches to surface roughness, in particular using a kind of parametric-homogeneous functions having fractal graphs, e.g., the Weierstrass-Mandelbrot functions, we clarify some common misinterpretations of statements by Berry and Lewis (1980). Finally in section 5 we discuss the models of surface roughness based solely on properties of the autocorrelation function or its Fourier transform (the PDSF). We argue that these models do not reflect tribological properties of surfaces. In particular, it is demonstrated that rather different profiles may have the same inclination of the PSDF in logarithmic coordinates, or they may have even the same PSDF. Thus, there is very small scientific value (if any) of a number of papers that model surface roughness only by its PSDF.

PROBABILISTIC CHARACTERISTICS OF SURFACE ROUGHNESS
If one considers a nominally flat surface and a plane perpendicular to the surface, then the surface profile is the crosssection between the plane and the rough surface. The rough profile may be presented as graph of a function z(x). Letz denote the mean profile line. If the origin level of the height measurements is taken atz then Let us mention here several popular height parameters: R max is the maximum height of the profile z(x) defined on an interval [−L, L], R a is the arithmetical mean deviation of the surface, and the root mean square (rms) height R q or σ 2 that is the square root of the mean square deviation with respect to the mean profile linē z = 0. The mathematical expressions for these parameters are given by where n is the number of points of measurements on the interval and z(x i ) is the measured height at the interval point x i . In addition, the arithmetic mean height R z is often used for practical applications. This parameter may be calculated as the average distance between the five highest picks and the five lowest points of the profile, i.e., One can introduce the density probability function φ(z). This function shows the probability that the height z(x) at a surface point x is between z and z + dz. Then, the expressions for R a and σ 2 in (3) is written as It is natural that the roughness parameters depend on the scale of considerations. For example, if z i = z(x i ) are measured with a stylus steps h then one can calculate the curvature of the profile peaks κ(x i ) However, it was found that the mean curvature varies depending on the sampling intervals (Greenwood, 1992). There was a hope that the fractal dimension could provide a scale independent parameter of surface roughness. It will be discussed later that actually this assumption was not justified. Kragelsky (1948) published one of the first papers where he claimed that the roughness heights distribution is Gaussian (normal). Although Zhuravlev (1940) presented the general expressions of his model for arbitrary bearing area curve, his example that employed linear dependence as an approximation of the bearing area curve was criticized by Kragelsky (1948) who wrote that the normal distribution of heights should be used.
We will not list here a number of parameters used in literature on tribology to describe the surface roughness. Various attempts to describe surface topography using several statistical parameters of surface roughness are described in detail in many books and papers, see e.g., Khusu et al. (1975), Nowicki (1985), and Whitehouse (2011). Some of these parameters are useful but most are not (Whitehouse, 1982).
An example of a very useful parameter is the Abbott-Firestone curve of the surface heights (bearing area curve) or the cumulative distribution function (z) For example, the (z) was used by Zhuravlev (1940). Let us demonstrate that for rough surfaces their contact properties are correlated with (z) that is equal to the length of a horizontal slice of the surface profile at the level h. Sometimes this curve can be used to estimate the force acting of a rough solid penetrating into an elastic foundation. Indeed, it if the characteristic size of contact region between a blunt punch and a thin elastic layer is larger than the layer thickness then the leading term of the asymptotic solution may be modeled as the Fuss-Winkler foundation (see e.g., Aleksandrov, 1963;Borodich et al., 2019b;Erbaş et al., 2019; and references therein). The Fuss foundation can be represented by a collection of independent springs attached to a rigid flat or as an elastic "mattress" (Winkler, 1867;Johnson, 1985) or as a punch acting on a liquid layer (Fus, 1801). Of course, there are restrictions on the use of the Fuss-Winkler foundation (Johnson, 1985;Popov, 2010). As Fus (more often his surname is written as Fuss) noted himself: "when a crumbly surface has a rigid substrate and does not have such depth as it should be for penetration of wheels into it by laws of hydrostatic, therefore the angle GOR" (the angle between the wheel and the surface flat) "should not be determined just by actual penetration depth but rather by the depth of the rigid support" (Fus, 1801). The Fuss-Winkler foundation was employed to model contact between an elastic foundation and (i) a nominally flat C B profile Mosolov, 1991, 1992); (ii) the rough hierarchical multi-level profile Onishchenko, 1993, 1999), and (iii) both nominally concave and convex fractal parametric-homogeneous punches (Borodich, 1998b).
The roughness parameters used may be formally divided in the following groups (Nowicki, 1985): (i) various parameters related to the shape of asperities, e.g., the rms curvature of the asperities; (ii) various parameters related to the asperity heights; (iii) horizontal parameters; (iv) parameters connected with amplitude of the asperities and their spatial extend, like the high spot count.

GAUSSIAN RANDOM PROCESSES AS MODELS OF ROUGH SURFACES
As it has been mentioned above, the surface roughness may be considered as graphs of random processes. The question if the process is Gaussian (normal) or not is crucial for estimation of validity of the models. Indeed, the overwhelming majority of papers on surface topography use statistical models of surface roughness based on explicit or implicit assumption of normality of roughness heights, see e.g., Fuller and Tabor (1975), Galanov (2011), and a discussion by Borodich et al. (2016).

Tests for Normality of Surface Roughness
There are many tests of normality (Thode, 2002). Each of these tests provides a quantitative estimation of proximity between the theoretical Gaussian distribution and an observed sample of measurements by producing the so-called p-value. The estimations are based on a particular test statistic. Borodich et al. (2016Borodich et al. ( , 2019a and Pepelyshev et al. (2018) used the most popular tests of normality to check normality of roughness of various surfaces. The afollowing tests were employed: Kolmogorov-Smirnov (KS), Anderson-Darling (AD), Cramer-von Mises (CVM), Shapiro-Wilk (SW), Shapiro-Francia (SF), Lilliefors (LF), and the Pearson. The p-value is a number that characterizes for the observed measurements, the significance at the scale [0, 1] that the normality hypothesis is true. It is possible to nominate the acceptable significance level. In our tests it was 5%. If the p-value is more than this level then the hypothesis should be accepted. This subject is a current challenging task in contact mechanics description and comprehension of nanorelated phenomena as highlighted by Carpick (2018).
Let us mention here a model of dry friction developed by Borodich and Savencu (2017). In this model molecular and chemical interactions are mainly connected to the nano-scale asperities, while mechanical interlocking between surfaces are connected to the micro-scale asperities. According to this model, for proper modeling of friction, one needs to get data about the surface roughness at both atomic/nano and at micro-scales.
Using the above mentioned tests of normality some data obtained for grinding surfaces was tested and results of the normality tests were negative at both nano and microscales (Borodich et al., 2016).
Then the test of normality were applied to surfaces of the epoxy resin replicas of polishing papers of various nominal asperity sizes. A white light interferometer (Zygo NewView 6000; Zygo Corporation, Middlefield, CT, USA) at a magnification of 50 was used for characterization of the surface roughness. The normality tests showed that the height distribution of the surfaces of nominal 0.3 and 1 µm are Gaussian (Pepelyshev et al., 2018). Finally, normality of roughness of carbon-based coatings deposited by direct current pulsed magnetron sputtering at two different substrate bias voltages was checked. The same as in the case of grinding surfaces, the roughness was measured at atomic/nano scales by AFM (Atomic Force Microscopy), while a profilometer was used for measuring of micro-asperity heights (Borodich et al., 2019b). It was found that surfaces at micrometer scale are normal. It is interesting to note that the AFM measurements with the 117 nm steps showed that the roughness of surfaces was Gaussian, while the AFM measurements with the 10 nm steps showed that the assumption of normality of roughness is not satisfied. This means that the use of the above mentioned statistical models of contact between nominally flat surfaces are justified only up to nanoscale, while their applicability at atomic and few nanometers scales may be questionable.

Stochastic Models of Surface Roughness
The normality of roughness is related to vertical distribution of profile heights. However, it is also very important to consider the horizontal distribution of the asperities. As it has been noted by Maugis (2000), two profiles may have the same peak height (local extrema) distributions, but the roughness may be rather different in the horizontal extension. If process is Gaussian then the horizontal properties may be fully described by the autocorrelation (correlation) function R of the process (2). Instead of R, one can use the PSDF (3), its Fourier transform G(ω). Khusu et al. (1975) presented results of very detailed studies of Gaussian processes in applications to surface roughness. In addition to (1), they considered three other correlation functions. An extended list of various correlation functions was presented by Dette et al. (2015). However, we would like to emphasize again that the above results are effective if the roughness is Gaussian. If the roughness is not normal then the properties of the sample paths are not fully determined by the mean and covariance functions (see e.g., Ghosal and Van der Vaart, 2017).
It is known (see e.g., Aldous, 1989 andBibby et al., 2005) that a non-Gaussian processes can be generated by the mean-reversing stochastic differential equation where x ≥ 0 and W B (x) is the standard Brownian motion (Wiener process). Choosing the appropriate value of the parameter µ and the function σ (·), we obtain a certain height distribution of the process X(x), while the auto-correlation function is ρ(x) = e −θ |x| and the power spectra is G(ω) = 2 π θ/(θ 2 + ω 2 ) for any choice of µ and σ (·). For example, the height distribution of X(x) has the gamma density with parameters λ and β if µ = β/λ and σ (X) = √ 2θ X/λ and the beta density with parameters α and β if µ = α/(α + β) and σ (X) = 2θ X(1 − X)/(α + β). Figure 1 shows simulated profiles with various height distributions and the same autocorrelation function R(x) = e −θ |x| .
In general, using (7), one can construct non-Gaussian processes with flexible auto-correlation function as the sum where X 1 (x), . . . , X m (x) are independent and each X i (x) is the solution of the stochastic differential equation Bibby et al. (2005). Such the process h(x) has the auto-correlation function and the power spectra The above examples support the statement that the mean and covariance functions of a random non-Gaussian process do not determine the finite-dimensional distributions of a random process (Gusak et al., 2012).

FRACTALS APPROACHES TO SURFACE TOPOGRAPHY
Concepts of fractal, fractal geometry and fractal dimension (FD) were introduced by Mandelbrot (1975). However, he did not give any definition of fractals. Later he defined fractal sets in a metric space saying "a fractal will be defined as a set for which the Hausdorff-Besicovitch dimension strictly exceeds the topological dimension" (Mandelbrot, 1977). Finally he withdrew the definition and suggested to use the term "without a pedantic definition" (Mandelbrot, 1983). This was the reason for a comment by Greenwood (1992): "Mandelbrot is somewhat reluctant to define 'fractals' or 'fractal dimension' preferring to offer examples." A popular description of many self-similar sets studied by fractal geometry, such as the Cantor staircase, von Koch snowflake, Sierpinski carpet, Menger sponge, along with Peano curve and other sets may be found in the book for school children published by Vilenkin (1968). Definitely Mandelbrot (1975) wrote his book under influence of Vilenkin's book. Indeed, at least 20% of Figures presented by Mandelbrot (1975) had analogs in the book by Vilenkin (1968). The above mentioned sets are out the scope of mathematical programmes for schools and for non-mathematical specializations of universities because the sets are more irregular than sets considered in common courses on Euclidean geometry. To encourage children for studies of set theory, Vilenkin used many funny stories and terms, e.g., he used the term the Devil staircase to describe the Cantor staircase. The Devil staircase term was later used and popularized by Mandelbrot (1983). Because there is no generally accepted definition, we will define fractals as objects with a non-integer FD. Evidently, the term FD must be defined separately. One of the authors (FB) and his co-authors have already published several reviews related to fractals, in particular a review related to the use of fractal concepts in fracture mechanics (Borodich, 1999). There are reviews of application of fractal ideas in contact problems (Borodich and Onishchenko, 1999;Borodich, 2013c) and several articles about the use of fractal concepts in tribology (Borodich, 2013a,b;Borodich and Evans, 2013). Borodich introduced the C B profile that is also known as the Cantor set model (Warren and Krajcinovic, 1996), and as the Cantor-Borodich profile, structure or fractal, see e.g., Abuzeid andEberhard (2007), Soldatenkov (2015), and Thielen et al. (2016). At least two kind of contact problems can be solved for a punch described by the C B profile Mosolov, 1991, 1992). The model was developed further by Borodich and Onishchenko (1993), Warren and Krajcinovic (1996), Plesha and Ni (2001), and others. Although Archard (1957) introduced the idea of hierarchical structure of roughness, there were no papers developing the idea of "bump on bump" structure of roughness until Onishchenko (1993, 1999) introduced the multilevel hierarchical model of roughness. Now this idea of hierarchy of surface structures is quite popular due to discoveries of Gorb and his co-workers along with other experts in biological objects. In particular, the hierarchical structures were discovered in adhesive setae of geckos (Gao et al., 2005) and spiders (Schaber et al., 2019); waterrepellent coatings of whip-spiders (Wolff et al., 2016); attachment discs of spiders (Wolff et al., 2015), and super-black snake skin (Spinner et al., 2013). Applications of fractal concepts to surface roughness were also discussed by Borodich and Galanov (2002) and by Borodich et al. (2016). Thus, here we will just remind some basic features of the fractal approach to surface roughness.

Mathematical and Physical Fractals
Unfortunately, many papers dedicated to applications of fractals do not provide definitions of used terminology. Often such papers are just a mess of vague discussions and nonjustified statements. To clarify the subject, we split fractals into mathematical and physical (natural) fractals (see e.g., Borodich, 1997). Both mathematical and physical fractals use the concepts of a cover. This means that the object (set) is covered by cubes of size at most or equal to δ. Fractal geometry deals with mathematical fractals. The mathematical methods of fractal geometry are described in many books and papers (see e.g., Falconer, 1990;Tricot, 1995) where various FD are studied in application to mathematical objects. Various FDs are used in the studies, mainly the Hausdorff dimension and box-counting dimensions usually attributed to Minkowski, Bouligand, Pontrjagin, Schnirelman, and Kolmogorov. These FDs can be calculated by taking the limits when δ → 0. In addition, fractal geometry term is also quite often applied loosely to a collection of semi-empirical or empirical methods for estimations of the FDs of objects.
If real world or numerically simulated objects demonstrate the power-law of the number-radius relation then these objects are physical fractals. The power-law of the number-radius relation is is the number of cubes of size δ used to cover the object, D is a FD of the object, * and δ * are the upper and lower cut-offs of the physical fractal law, respectively. The former relation of (8) is used when the cover size δ is varied and the object size R is fixed, while in the latter relation is used when the cover size δ is fixed and the object size R is varied. In this latter case R * and r * are the upper and lower cut-offs. The the slope of linear approximation of ln(N(δ * )) against ln(R) is used to estimate the value of D. The main distinction between these kinds of fractals is the following: the physical fractal behavior (8) is observed on a bounded region of scales only, while to study mathematical fractals one has to scale of consideration to zero limit.

Self-Similar Sets
Introduction to fractals starts often by presenting self-similar sets mentioned above (see examples given by Vilenkin (1968) and Mandelbrot, 1975). Due to such examples, there is a myth that all fractals are self-similar. However, self-similar sets are a very specific kind fractals. In general, self-similarity is not related to mathematical fractals. Their scaling properties are based on scaling of fractal measure or quasi-measure (see a discussion by Feng, 2010, andBorodich, 2019), while for physical fractals their scaling properties are reflected by the relation (8).

Some Specific Features of Mathematical Fractals
Speaking about mathematical fractals, one has to specify the FD used. It is wrong to say that all FDs are equal to each other. For example, the Hausdorff dimension dim H S of a set S may be not equal to its box-counting dimension dim B S. However it is known that dim H S ≤ dim B S.
The statement that all nowhere differentiable curves are fractals is wrong. For example, Borodich (1998a) introduced the Moscow University function U M that is nowhere differentiable, while it is a non-fractal curve and each finite subinterval of [1, 2] contains an infinite number of reduced copies of the whole function, i.e., it is a self-similar curve.
A mathematical fractal curve has an infinite length. Borodich (1997) formulated the following paradox in application to fractal fracture: if the crack surface is modeled as a mathematical fractal and the concept of the Griffith surface energy is used, then the propagation of the fractal crack is impossible.
If a surface is smooth then it is not a mathematical fractal (see e.g., Falconer, 1990). Even if a mathematical fractal curve is continuous everywhere, it is nowherehas differentiable. Therefore, it is often very difficult to formulate a boundaryvalue problem for solids having fractal boundary (see a discussion by Borodich and Volovikov, 2000). Indeed, such nowhere differentiable profile do not have even normal vector.

Some Specific Features of Physical Fractals
If the FD is specified then it is convenient to use the fractional part of the FD D * . Then the FDs of fractal profiles and surfaces are 1 + D * and 2 + D * , respectively. Sayles and Thomas (1978) presented the experimental relations between normalized PSDF and wavelength in logarithmic coordinates as a very impressive united line for many different surfaces. One could think that the relation spans from micrometers to many meters. However, Berry and Hannay (1978) noted that for each individual type of topography the span was over much smaller ranges. They argued that these experimental results could be represented as a united line lg(G) ∼ lg(ω) due to a specific normalization of the results. Further, researchers from Jerusalem presented results of an analysis of data published by Physical reviews journals and showed that in average the physical fractals span about 1.5 orders of magnitude . Mandelbrot (1998) argued that the published limited-range examples of power law correlations may be explained as unfortunate side effects of enthusiasm, imperfectly controlled by refereeing. However, the Jerusalem group disagreed with his arguments and believed that the limited-range relations for natural fractal are the dominant fractals observed in Nature . However, if the value of the FD is stable for less than two or three decades then fractal concept is not useful, as it was stated by Whitehouse (2001).
In general, the FD values obtained by various practical methods are not reliable (see a discussion by Borodich and Evans, 2013). For example, the power spectrum method is often employed for estimation of FD values of self-similar or self-affine fractals, in particular rough profiles (see e.g., Dubuc et al., 1989;Schmittbuhl et al., 1995). The method is based on a not very accurate statement when the power spectrum obeys a power law G(ω) ∼ c/ω α ... it is reasonable to expect a signal with a ω −α power spectrum to have a graph of dimension (5 − α)/2 (see Falconer, 1990). Normally the exponent α is in the range 1 < α ≤ 3.
As Whitehouse (2001) noted there is a very small spread of the FD values obtained for surfaces produced by different manufacturing processes. In addition, there is no well-established algorithm for estimations of the fractal law cut-offs. Thus, the physical significance of the fractal approach is very limited. We can add that if one can attribute the fractal scaling for a small range that spans for just 1.5 or 2 decades then fractals do not provide a scale-independent description of surface roughness.

Parametric-Homogeneity
One of the authors (FB) and his co-authors introduced different types of fractal profiles that allowed to handle rough fractal surfaces rigorously. There were introduced : (i) the C B profile Mosolov, 1991, 1992), (ii) multilevel Hierarchical profile Onishchenko, 1993, 1999), and (iii) graphs of parametric-homogeneous functions (Borodich, 1993(Borodich, , 1998a. The first two profiles were based on an iterative procedure and they has been already mentioned above, while PH-functions will be described below. Note that the general statements obtained by Borodich and his co-authors about the contact problems for all these three types of surfaces do not depend on the statements about their FDs, i.e., they are valid for both fractal and non-fractal cases. The same author (FB) introduced also the concept of parametric-homogeneity. The concept includes parametric quasi-homogeneous (PQH) transformations (in particular parametric-homogeneous (PH) transformations), corresponding functions, and PH-and PQH-sets. The PH and PQH-functions can be considered as a natural generalization of concepts of homogeneous and quasi-homogeneous functions. If the latter kind of functions is based on the use of the classical scaling with arbitrary positive scaling factor λ, PH-functions and PQH-functions are based on the use of discrete self-similarity with a fixed rescaling parameter p. Let Z be the set of integer numbers. Then one can employ the discrete group of coordinate dilations (Ŵ p αk ) apply the following PH-transformation Ŵ p αk x = (p kα 1 x 1 , ..., p kα n x n ), p > 0, k ∈ Z.

Parametric Quasi-Homogeneous Functions
The parametric-quasi-homogeneous function of degree d and parameter p with weights α = (α 1 , . . . , α l ) is denoted by B d : R l → R and it obeys the following scaling law for a positive parameter p, p = 1 and any k ∈ Z The parameter should be unique in some domain (Borodich, 1998a,b). I α 1 = ... = α l then PQH-functions is called a PH-functions. It is clear that if p is a parameter of the PHtransformation then p k can be also used as a parameter. Hence, as the parameter, the least p : p > 1 is taken. The graphs of these functions can have various properties, in particular they can be smooth or fractal.

Weierstrass-Mandelbrot Functions
Mandelbrot (1977) generalized the Weierstrass function, whose graph is continuous everywhere and differentiable nowhere, and introduced the complex value Weierstrass-Mandelbrot (W-M) function W(x) and its particular real case C(x; p) C(x; p) = ∞ n=−∞ p −βn (1 − cos(p n x)), p > 1, 0 < β < 1 (9) Frontiers in Mechanical Engineering | www.frontiersin.org where φ n are arbitrary phases. Properties of these functions were studied in detail by Berry and Lewis (1980). The box dimension of C(x; p) graphs is D = 2−β. There is no rigorous mathematical proof that its Hausdorff dimension is the same.
It is evident that C(x; p) is a particular PH-function (9) of degree β. As it has been mentioned above, the graph of C(x; p) function was used often for modeling rough profiles.
In addition, Berry and Lewis (1980) were first who studied the discrete power spectrum G(ω) of the W graphs (9) and approximated it by a continuous spectrumḠ(ω) ∼ ω −(5−2D) . This is in agreement with the above mentioned statement by Falconer (1990).

Self-Affine Scaling and the Hurst Exponent Terms
Many papers dedicated to studies of surface roughness mention the concept of self-affine fractals without any definition of the concept. However, one should be very careful using such terms as self-affine scaling. Actually, it is known from mathematical courses that mapping on a plane is self-affine if it is a oneto-one mapping whose images of any three collinear points in turn belong to a line. Under this mapping the angle between two straight lines is normally not preserved. A particular case of such mapping is the quasi-homogeneous (QH) or anisotropic coordinate dilation where α 1 and α 2 are weights, λ x = e α 1 and λ y = e α 2 are scaling factors. This means that if a coordinate system x, y was Cartesian one then the x ′ , y ′ system is also rectangular. However, papers on fractals reduce usually the term self-affine mapping to the above particular QH case.
In 1980 the concept of fractal was very novel. Perhaps for popularization of the concept, they used some explanations and statements that were not mathematically rigorous. For example, studying W-M functions, they wrote The fractal nature of W implies that repeated magnification of its graph reveals ever-finer levels of detail.... The levels of detail are self-similar under an affine scaling in which the x axis is stretched by a factor p and W axis by p 2−D [note we have changed γ and t used by Berry and Lewis (1980) to p and x, respectively]. Unfortunately, these statements of Berry and Lewis (1980) were misinterpreted and used in many papers as generally accepted and mathematically rigorous statements applicable to all so-called self-affine fractals.
There was an attempt by Mandelbrot (1986) to develop further the concept of self-affinity. He formulated several statements about self-affine fractals that were not supported by any strict definition of these fractals. In essence, his statements are: there are both local and global scaling for self-affine fractals and each version of FD for self-affine fractals has a local and a global value, separated by a crossover. However, these statements have no mathematical meaning. Schmittbuhl et al. (1995) wrote that the scaling (10) means that the heights h are homogeneous functions of the scaling factors (putting λ x = λ y ) where the homogeneity exponent H is the self-affine exponent, or the Hurst exponent. Unfortunately, the statement by Schmittbuhl et al. (1995) about the homogeneity of the fractal graphs is not very accurate, while the meaning of the Hurst exponent is illdefined. It would be better if they would say that the fractal roughness to obey a discrete dilation homogeneity, i.e., it obeys the PH-law. Indeed, the scaling properties of C(x; p) function are connected with the discrete group of coordinate dilation Although this scaling is quite often attributed to FD of the C(x; p) graph, saying that H = 1 − D * , these statements are wrong because the scaling properties are governed by the degrees d of PH-functions. (Borodich, 1998a,b) have the same values of the Hausdorff and box dimensions and absolutely different trends. On the other hand, if one uses the above self-affine terminology then b β (x, p) and b 2 (x; p) are self-affine fractals, while b 1 (x; p) is self-similar fractal because b 1 (p k x; p) = p k b 1 (x; p).
It follows from the above discussion that Hurst exponent is an ill-defined term and it is not connected to FDs of the PH-graphs. A statement by Mandelbrot (1983) is often cited. Mandelbrot (1983) wrote that coastlines are not circles, clouds are not spheres, and mountains are not cones, assuming that they have to be modeled as fractals. Borodich and Onishchenko (1999) extended Mandelbrot's statement saying that roughness of real bodies is not a mathematical fractal and all these geometrical objects: spheres, cones, circles as well as fractals are only mathematical idealizations of complex shapes of natural objects. Mathematical and physical fractals should not be confused. Certainly one can employ mathematical fractal as a possible model that reflect the power-law number-radius relation or PSDF of a natural object within a bounded interval of scales. However, the obtained problem may be a very complicated.

POWER SPECTRAL DENSITY FUNCTION APPROACHES TO ROUGH SURFACES
The last 20 years or so the PSDF approach became very popular in tribology community (see e.g., Persson, 2006). In fact, the authors of the papers reduce the properties of rough surfaces just to their PSDF G(ω). This approach was criticized by Borodich (2002). He wrote: two punches having the same fractal surface but FIGURE 2 | The power spectra of the profile h(x) = sin(2π ln(x)/ ln(p)) for x ∈ [0.0001, 1] with p = 2 (black) and p = 1.3 (red). situated either above or below the surface show usually different asymptotics in both load-displacement and load-area relations. As an example he considered C B profile and solutions presented by Borodich and Onishchenko (1993). Then Borodich et al. (2016) showed that any surface and its replica have the same power spectral density, therefore one cannot characterize contact properties of rough surfaces just by surface "spectrum." Indeed, one can consider an absolutely flat smooth surface having many sharp dents and its complementary replica (an inverted replica): a surface having many sharp asperities of the surface roughness. Evidently, they gave absolutely different contact properties (see e.g., examples by Borodich and Onishchenko, 1993).
Moreover quite often the PSDF approach is reduced to the socalled self-affine surfaces by employment of the Hurst exponent term. As we have discussed above the term is ill-defined and such papers just create "new clothes" for the vulgar fractal approach.
Using examples, we will show below that the many common statements about rough non-Gaussian surfaces are wrong or just meaningless.

Smooth Functions May Have a Power-Law PSDF
Let us show that connections of FDs and the slopes of PSDFs in logarithmic coordinates are meaningless. As an example, consider a very simple and smooth PH-function, namely a sin log-periodic function. Figure 2 shows the power spectra for h(x) = sin(2π ln(x)/ ln(p)) with x ∈ [0.0001, 1].
It is clear that the function is non-fractal, nevertheless its power spectra in logarithmic coordinates located along a straight line within a bounded range of scales.

PSDFs Are Almost the Same for Different Truncated Weierstrass-Mandelbrot Functions
Let us consider two truncated Weierstrass-Mandelbrot functions where N is large natural number, p > 1 and D ∈ (1, 2). In addition, let us consider a function x D−2 C t (x; p). If N → ∞ then C t (x; p) → C(x; p) and D is their box-dimension, however they have different trends.
Let us fix the values p = 1.5 and D = 1.5 and vary N. Let us calculate numerically the PSDFs for obtained graphs. The results are presented in logarithmic coordinates in Figure 3.
One can see that the average PSDFs are approximately the same while the functions are rather different. Now let us fix the values p = 1.5, and N = 100 and vary D and calculate numerically the PSDFs for obtained functions. The results are presented in logarithmic coordinates in Figure 4. Note that for N = 100, the graphs of C t (x; p) and x D−2 C t (x; p) are very close to graphs of PH-functions of degrees d = 2 − D and d = 0, respectively. We know that these functions have the same FD D and different trends (see a discussion above), while their PSDFs are very similar to each other.
The same procedure can be applied to the profiles A(x; p) and x D−2 A(x; p). The results are shown in Figure 5. The results lead to the same conclusion as above.

PSDF Does Not Characterize Tribological Properties of Polished Surfaces
It is known that the tribological properties of gears can be greatly extended by improved surface finish (Krantz et al., 2001). This fact was employed by Harris et al. in a series of papers dedicated to the use of modern hard carbon-based coatings. Indeed, these hard carbon-based films can remove protuberances on the steel counterparts, reducing the values of high intensity stresses and the total number the stress concentration points. On the other hand, due to polishing of the coating asperities, the rate at which the coatings abrade steel declines as a power-law of the cycle numbers, i.e., very rapidly. The coated surface become super smooth even after just 500 cycles (see SEM images by Borodich et al., 2003). Let us consider the following thought experiment. An intact surface is modeled by a Gaussian function f 1 (x) shown in Figure 1A (this assumption is in accordance with roughness measurements of intact surfaces by Pepelyshev et al., 2018 andBorodich et al., 2019a). Then the asperities of the surface have been polished away and the roughness of the polished parts are described by 0.1f 2 (t) where the profile f 2 (x) is shown in Figure 1B. Hence, the full polished surface is described by f 3 (t) = min[f 1 (x), 0.1f 2 (x)] shown in Figure 6. Figure 7 shows the power spectra for the profile f 1 (x) (black), the profile f 2 (x) (red), and the "polished" profile f 3 (x) = min[f 1 (x), 0.1f 2 (x)] (blue). These three profiles have the different outlook and tribological properties, however have the same exponent 1.9 in the power spectra, i.e., G(ω) ∼ 1/ω 1.9 .     Figure 1A and the profile f 2 (x) is shown in Figure 1B.

CONCLUSIONS
A review of mathematical approaches for description of topography of engineering surfaces has been presented. It has been shown that in spite of huge amount of parameters used to characterize the surface topography, only some of the parameters are quite useful. However, their use is rather bounded, e.g., these parameters may be helpful at meso or even microscales, however they could be useless at the nanoscale.
There are many models of random processes, however just the case of Gaussian processes is well-elaborated. Some statistical methods employed for checking normality of distributions of asperity heights are reviewed. The preliminary results showed that the intact surfaces are quite often Gaussian at both micro and nanoscales, while the grinding surfaces are not normal.
Then the fractal approaches to surface roughness are discussed. Some common misinterpretations of statements by Berry and Lewis (1980) are also discussed. Some disadvantages of the fractal approaches and commonly reported wrong statements about fractals are listed. It is argued that the practical usefulness of the fractal approaches is rather doubtful. One should not expect that the employment of a mathematical fractal model of a rough surface will provide considerable advantages. In fact, such models may cause many mathematical difficulties. Thus, a strict approach to fractal modeling may substitute a difficult problem to another more difficult than the original one. Further, dimensions of physical fractals cannot be used as scale independent parameters. One should provide proper explanations of the fractal concepts used, otherwise this could lead to misinterpretation of the results. Unfortunately, quite often the answer to the question by Jelinek et al. (1998): Is there meaning in fractal analyses?, is "No." Finally the models of surface roughness based solely on properties of the auto-correlation function or its Fourier transform (the PDSF) are discussed. It has been noted that the PDSF approach to non-Gaussian surfaces are in essence the new clothes for the vulgar fractal approaches. We argue that these PDSF models do not reflect tribological properties of surfaces.
In particular, it is demonstrated that rough profiles may have the same slopes of the PSDF in logarithmic coordinates, while they have rather different tribological properties. Thus, there is a bounded scientific value (if any) of a number of papers based on the PSDF approach.

AUTHOR CONTRIBUTIONS
FB: substantial contributions to the conception of the work, analysis, and interpretation of published data for the work. XJ: analysis and interpretation of published data for the work. AP: statistical analysis of the results and numerical simulations for the work. All authors contributed to the article and approved the submitted version.