Numerical and topological description of contact areas at different size scales for the contact of rough solid surfaces

Abstract

This study introduces a novel heuristic phenomenological model for analyzing the evolution of contact areas on rough surface. Contrasting with traditional methods, it employs a cut-off threshold approach to track numerical and topological metrics across different deformation stages. The model quantifies contact area distributions, nested sub-regions, and self-affine parameters, revealing universal trends across scales spanning nanometers to kilometers. Metrics for synthetically generated isotropic surfaces with Hurst exponents H = 2.5 and 3.5 correlate closely with those from AFM and SEM experimental datasets, respectively. In addition, the model has been tested on NASA's SRTM datasets. Cross-correlation demonstrate significant similarities in numerical and topological metrics across diverse measurement techniques, surface types, and scales, highlighting the method's robustness and calibration-free scale invariance. This approach bridges gaps in multiscale tribological analysis, offering deeper insights into frictional transitions and surface interactions. Beyond tribology and materials science, this general approach enables fundamental characterization of surface morphology as such, making it applicable to diverse fields including geomorphology, biomimetics, and nanotechnology.

---

New concepts

We introduce a scale-invariant computational approach to characterize contact areas between rough surfaces, leveraging topological and geometric descriptors such as Hausdorff dimension, lacunarity, and object nesting. By applying a threshold slicing model to diverse surface data – from nanoscale AFM data to meter-scale Earth's geodetic topographies data – we reveal universal trends in contact geometry during simulated frictional transitions. Our method transforms surface relief data into a binarized multiscale representation, enabling automated extraction of numerical and topological features that remain consistent across material types, measurement techniques, and size scales (10⁻⁹ to 2000 m). These descriptors offer a new pathway to integrate topological data analysis into AI-driven surface characterization pipelines. The observed scale-invariant behavior suggests a fundamental statistical underpinning of frictional contact evolution, paving the way for predictive modeling of surface interactions. This work bridges classical tribology with non-equilibrium physics and machine learning algorithms, offering a versatile tool for autonomous materials discovery and the engineering of frictional interfaces.

Introduction

The significance of our work is underscored by recent advancements in nanoarchitectured materials, such as crumpled graphene oxide (GO) membranes, which demonstrate how tailored nanostructures can achieve exceptional performance in gas separation.¹ These membranes exhibit a unique combination of high permeability and selectivity due to their anisotropic nanochannels, highlighting the critical role of controlled nanoarchitecture in overcoming traditional material limitations. Similarly, our topological approach to surface analysis bridges multiscale tribological behavior with universal metrics, offering a calibration-free, scale-invariant framework. By leveraging hierarchical surface features and fractal geometry, we provide a robust methodology that transcends conventional asperity-based models, enabling precise characterization of contact mechanics across diverse applications – from nanoscale friction to geomorphology. This aligns with the broader paradigm shift in materials science, where engineered nanostructures unlock unprecedented functionality, as exemplified by GO membranes, reinforcing the transformative potential of our topological model.

Dry friction is a universal phenomenon, yet deriving characteristics of frictional behaviour from fundamental physical principles remains challenging.^2,3 Both static and sliding friction are typically mediated through individual asperities, which interact between surfaces. As a collective phenomenon, friction is characterized by significantly reducing the degrees of freedom associated with individual asperity contacts to a limited set of multiscale degrees of freedom.⁴

There has been a long-standing debate in tribology about how to reconcile elastic contact of rough surfaces with Amontons’ first law. Although the classical Hertz solution⁵ for a single elastic contact predicts a nonlinear dependence of the area on the load (A ∝ W^2/3), Greenwood and Williamson (GW)⁶ showed that for real surfaces with asperities distributed over the height, the contact area almost becomes linearly dependent on the load. Their model, which considers the surface as a set of spherical protrusions subjected to elastic deformation, was further developed in.^7–16 An important contribution was made by Archard,¹⁷ who demonstrated that a multi-level hierarchy of asperities (protrusions on protrusions) leads to a power-law dependence close to linear (A ∝ W^N, where N → 1 as the number of levels increases). His ideas, which anticipated fractal concepts,¹⁸ explain why even a purely elastic contact can yield a nearly linear A(W) relationship. Subsequent studies of the atomic friction, including those of Persson,¹⁹ have confirmed that this proportionality holds in the presence of adhesion. Thus, the current consensus is that elastic contact of rough surfaces does indeed obey Amontons’ law.²⁰ While successfully explaining Amontons’ law through emergent linear load-contact area relationships,²¹ the GW model's assumption of monoscale asperities limits its applicability to real multiscale surfaces. This limitation motivated Majumdar and Bhushan's fractal reinterpretation,²² where the height-truncation paradigm was extended to multilevel surface “islands”. The advancement naturally incorporated surface hierarchy effects that Archard had anticipated theoretically, but which became mathematically tractable only through fractal concepts.^23,24 The evolutionary path from GW statistics to fractal geometry reflects tribology's ongoing effort to reconcile contact mechanics with increasingly sophisticated surface characterization data.

The height-cutoff approach in statistical models reduces complex surface interactions to discrete microcontacts occurring at the most prominent asperities. This formalism reveals that the real contact area (A) exhibits near-linear dependence on normal force (F_N) when surfaces contain sufficient asperity density. The underlying mechanism involves two competing effects: as load increases, existing contacts expand while new asperities enter the contact zone, with the decreasing stiffness of individual contacts being compensated by their growing number. This dynamic maintains approximate proportionality A ∝ F_N, providing the foundation for Amontons’ friction law. The Bowden-Tabor model offers an extreme case where purely plastic deformation enforces strict linearity (A = F_N/H) through constant hardness pressure.²⁵ Real contact conditions typically show transitional behavior: initial plastic crushing of sharp asperities during run-in reduces contact pressures to elastic levels, after which the system maintains quasi-linear area-load dependence due to asperity height distribution.²⁶ Consequently, both elastic and plastic regimes demonstrate empirically observed linear trends in friction force versus load, differing primarily in their resultant friction coefficients after run-in stabilization.

While successful, the GW model relies on several simplifying assumptions. First, it treats asperities as non-interacting, which holds for small compressions but fails when neighboring contacts merge under significant loads. Nayak demonstrated that contact count peaks and then decreases due to coalescence, a phenomenon unaccounted for in the original formulation.^27–29 Second, the model assumes uniform asperity curvature, whereas real surfaces exhibit Rayleigh-distributed curvatures (Whitehouse & Archard, 1970s).³⁰ Third, its parameters are scale-dependent – higher resolution measurements reveal smaller asperities with sharper curvatures, necessitating spectral integration limits. Most critically, the basic model neglects adhesion and plastic deformation. Later extensions (e.g., CEB model, 1987;³¹ Kogut-Etzion, 2002³²) introduced plasticity criteria while retaining the height-truncation approach, but still ignore contact interactions and macroscopic surface geometry. These limitations persist in most statistical asperity models.

The Mandelbrot's fractal geometry in the 1970s–80s emphasized inherent self-affine properties across multiple scales.³³ The scaling behavior, first rigorously demonstrated by Sayles and Thomas (1978), persists across approximately ten orders of magnitude until reaching material-dependent cutoffs determined by microstructural features like grain size or manufacturing signatures.^34–36 Beyond this characteristic roll-off scale, the surface statistics transition from fractal scaling to different behavior regimes. The essential insight is that real surfaces exhibit this hierarchical, self-affine geometry until encountering fundamental physical limits of the material or measuring system.

Fractal approaches fundamentally transform our understanding of surface contact by modeling roughness as a continuous height function z(x,y) characterized by spectral properties rather than discrete asperities. The foundational work of Majumdar and Bhushan (1990) demonstrated that intersecting fractal surfaces with horizontal planes produces contact areas following power-law distributions, achieved through their innovative use of the Weierstrass function.^37,38 Their approach effectively translated multiscale roughness into contact statistics by representing each contact region with equivalent smooth protrusions. These developments built upon Archard's pioneering hierarchical concept from the 1950s, which revealed how adding microasperities progressively modifies the contact area-load relationship from the classical Hertzian A ∼ F^2/3 toward linearity, approaching A ∼ F^8/9 with one additional level and A ∼ F^26/27 with two levels.²³ Archard's work remarkably presaged fractal concepts by showing how infinite hierarchies naturally produce near-linear contact area dependencies.³⁹ Contemporary fractal theories, exemplified by Persson's (2001) model, employ spectral density formulations to solve contact problems across all relevant scales.^11,12 While predicting nonlinear behavior at minimal pressures, Persson's framework confirms the approximately linear A(F) relationship for elastic non-adhesive contacts that aligns with both experimental observations and GW model predictions. This approach's mathematical formulation as a diffusion-like equation provides unique capabilities to simultaneously predict contact areas, pressure distributions, and interfacial gaps – a significant advancement over traditional asperity-based models.

Fractal analysis reveals the scale-dependent nature of contact area, where apparent continuous contact patches decompose into smaller islands with increasing magnification. For ideal self-similar surfaces, this hierarchical structure suggests vanishing true contact area, though real materials impose physical limits through either plastic deformation at small scales or atomic discreteness. This framework explains the resolution-dependence of measured roughness parameters – each tenfold increase in magnification typically reveals an order-of-magnitude increase in surface slope and curvature, as demonstrated by Jacobs et al.'s (2018) nano diamond studies showing growth in measured roughness from millimeter to nanometer scales.⁴⁰ The review by Jacobs and Pastewka (2022) suggests that surface topography can be viewed as a material parameter to establish the structure–properties relationships.⁴¹

Among new techniques are various data science methods, such as topological data analysis (TDA), machine learning (ML) and artificial neural networks (ANN).⁴² The development of materials on the Frontier of materials science requires a new methodological framework for describing the relief of real surfaces on nanoscale.^43–45

The emerging field of triboinformatics bridges classical tribology with modern data science, employing machine learning to extract complex patterns from tribological datasets.⁴⁶ Recent work demonstrates how these data-driven approaches can predict wear rates and friction coefficients by analyzing subtle topographic patterns imperceptible to traditional parameters.⁴⁷ The strength of triboinformatics lies in its ability to process multi-scale surface information without relying on pre-defined physical models.

Recent applications demonstrate TDA's superiority over conventional roughness parameters.⁴⁸ Studies by Yesilli et al.⁴⁹ reveal that persistence-based metrics can distinguish surfaces with different wear histories and predict functional performance more effectively than traditional R_a or R_q values. The approach bridges the gap between discrete contact mechanics and continuous surface characterization, enabling more accurate predictions of tribological behavior across scales.

These approaches analyze complete topographic signatures rather than relying on simplified parameters, capturing crucial multiscale interactions between surface features and contact behavior neural networks trained on Hertzian contact solutions enable instant predictions for arbitrary rough surfaces.^50,51 These developments are creating a new paradigm of “digital twins” in tribology, where comprehensive surface characterization (combining spectral, topological, and morphological data) enables precise, real-time simulation of contact behavior, wear processes, and lubrication performance. This represents a fundamental shift from conventional roughness parameters to multidimensional surface descriptors in predictive modeling.^52–57

Understanding surface roughness is the key to understanding multiple-asperity friction mechanisms.⁵⁸ Usually, mechanical contact occurs on the tops of the asperities so that a significant portion of the interface between two contacting solids is in the “separation“ state, while only a small fraction is in the contact state. The local normal, σ(x, y), and shear, τ(x, y), stresses act at the interface. At certain contact points, the frictional condition relating the normal and shear stress components, τ = μσ (where μ is the coefficient of friction), may be satisfied, while at other points τ < μσ. Consequently, locally, the contact at a point can be in the slip or the stick state. Global sliding is possible when the entire contact zone is in the slip state. The initiation of the global sliding with increasing shear load force applied to the solid is related to surface roughness properties.

Consider a simplified model of contact of a smooth rigid surface vs. a rough elastic-plastic surface characterized by a certain random surface profile z(x,y). While it is difficult to evaluate the map of the slip zones at a frictional interface, one can obtain an experimental height profile of a surface, z(x,y). When surfaces are pressed together by the normal load force W, in the absence of shear force, the asperities of a rough surface deform with the stress being roughly proportional to the height of the asperity, σ(x,y) ∼ z(x,y). Given this approximate proportionality, one can argue that the condition of the local transition to slip, τ(x,y) = μσ(x,y) in the presence of both normal force W and the shear force F, could be replaced by a similar condition τf = K z(x,y), where τf is the frictional strength of the interface, and K is a proportionality coefficient. Here, the profile map z(x,y) corresponds to the case when only the normal load force W is present. In this way, the study of the transition cab be replaced by the study of the function z(x,y).

The abrupt and rapid stick-slip transition has features similar to the critical behavior.⁵⁹ Two objects in a state of contact with the illustrated forces are shown in Fig. 1A, an illustration of the evolution of the shape of the contact pads depending on the applied normal force is shown in Fig. 1B. Due to the presence of asperities, with increasing shear load, the individual contacts are more likely to transform from the stick to the slip state, i.e., the state when the critical condition for the shear-to-normal stress ratio is locally satisfied. The stress state at every point of contact depends both on the shear load at that point and, due to the elastic relaxation, on the states of neighboring points thus forming slip zones of correlated neighboring points. The “islands” of the slip zones grow with the shear load until they merge together forming a network of slip zones. This makes the transition to the global sliding of the two bodies in contact similar to classical problems of percolation.

Fig. 1 A schematic showing various concepts related to the stick-slip transition. (A) The similarity of frictional sliding of two blocks subjected to the normal load, W, and tangential load, F; (B) typical stages of a 2D percolation-type behaviour. With increasing normal load, the size of slip zone spots (black) increases. A transition to the global sliding is expected when the slip zone forms a simply connected area.

A rough textured solid surface is a complex object; therefore, it is often difficult to completely characterize surface texture with one or several parameters. Different roughness parameters are needed for various applications. Surface texture parameters are standardized to characterize surfaces based on both their roughness and broader morphology. International Organization for Standardization, 251786⁶⁰ and 42877,⁶¹ Japanese Industrial Standards B 06018,⁶² and the German Institute for Standardization DIN 47629,⁶³ Russian GOST 1156-78⁶⁴ offers standards for a comprehensive framework for assessing surface characteristics. Within these standards, a wide range of varied surface texture parameters is defined, encompassing both profile-based and area-based measurements. However, this is not enough to describe gas–liquid–solid interactions, for example, wetting of rough surfaces or superhydrophobic phenomena, where other parameters, such as the Wenzel roughness factor are used.⁶⁵ Moreover, different surfaces can be described by the same value of roughness parameters. For example, the same value of mean roughness which only characterizes the height of surface asperities may correspond to surfaces with very different horizontal surface morphology. A textured surface is characterized by a various number of parameters that can be represented as belonging to certain high-dimensional data space.⁶⁶ TDA is used for dimensionality reduction of complex higher-dimensional surface texture data and inferring relevant features and properties of the surfaces. In our previous works, among the surfaces of brass samples, structured by cavitation bubbles, layer-by-layer polyelectrolyte assemblies, WO₃ thin films, biopolymer thin film and wrinkling in graphene oxide membranes, examples of complex data sets in terms of roughness topological and geometrical characteristics were found.^51,68–72

There are many methods of characterizing the texture of the surface: optical and mechanical profilometry, atomic force microscopy (AFM) with different modes, scanning electron microscopy (SEM), optical microscopy with different modes, and replicas and comparators (comparison standards) with optical devices. Among these methods, the most accurate is the AFM method. It is most often used in scientific research to describe the topography of a surface. However, the AFM method requires much time to obtain results (unlike SEM), a fine selection of the type of measurement mode, a highly qualified operator, and long-term interpretation of results. The disadvantages of the AFM method also include susceptibility to electromagnetic, acoustic, and vibration noise.

It is known that developed surfaces occur at a variety of scales: from micro to macro scales. For example, the surface of the Earth can also be described in terms of surface roughness. Our new approach are used to analyze digital data, making it universal regardless of scale. In this work, we implement a new approach to analyze the contact area for each step of tangential cutting-off (threshold) of the surface along the height and show the universality of the new methodology for the study of data received using AFM and SEM techniques, synthetic maps of height, moreover, geodetic data areas of Earth relief taken via NASA's Space Radar Topography Mission (SRTM). Using geodetic data demonstrates the scale universality of the proposed approach. In this study, we apply a novel mathematical methods of surface roughness analysis to rough surfaces at three different scales: nanoscale (AFM), microscale (SEM), macroscale (geodesic data of Earth surface), and scale-less synthetic isotropic surfaces to investigate common characteristics of rough surfaces at these scales.

Materials and methods

Atomic force microscopy datasets and measurements

To demonstrate the operation of the algorithm and the computational results, 6 AFM datasets from our previous studies were selected. Due to the fact that all the data sets used in the work—the results of atomic force microscopy—have already been presented in our other articles, we will provide in this article only a brief description of each of the data sets, refraining from a detailed description of the methods for obtaining these objects. For a detailed introduction to the materials and methods used to make these objects, we kindly ask the reader to refer to the articles listed in the reference section, where extensive information on experimental procedures techniques can be found. First dataset is the results of semiconductor wafer Si with smooth relief. Next, two of them relate to polyionic assemblies (polyelectrolytes multilayers PEM), containing 4 bilayers of polyethyleneimine (PEI)/poly (4-styrene sulfonate) sodium salt (PSS) bilayers on a smooth silicon wafer.⁶⁹ Forth one is ultrasonically threated brass surface.⁶⁷ One more AFM data is modulated magnetron sputtered WO₃ thin film,⁶⁸ and the last one is bio-based spin-coated poly (3-hydroxybutirate) (P3HB).⁷⁰ Detailed information regarding the preparation methods for these samples can be found in the corresponding references, where the protocols and conditions used are thoroughly described, providing a comprehensive understanding of the sample preparation process.

The topography of all samples was measured by scanning probe microscope NTEGRA Aura (NT-MDT, Russia) and HA_NC probes with tip diameter of approximately 20 nm. AFM was employed to scrutinize the sample surface in both semi-contact and contact scanning modes. The semicontact mode was chosen as primary mode for its minimal invasiveness to the samples, operating frequency was approximately 140 kHz and an average probe force constant was about 3.5 N m⁻¹. Image resolution was 256 × 256 points, corresponding to a point-to-point spacing of approximately 20 nm over a scanning area of 5 × 5 μm². For the P3HB relief data set, the data set was cropped from 512 × 512 to 256 × 256 points. The scanning speed was maintained at approximately 8 μm s⁻¹. The study was carried out in air conditions with an average temperature range of 21–23 °C and a relative humidity range of 30–40%. To mitigate the impact of vibration noise, both passive and active vibration protections were employed during measurements. For precise alignment of the probe and selection of the AFM measurement area, an optical registration system was utilized. AFM dataset, statistics, and results of modelling provided in Fig. S1–S4 and Table S1. One example of AFM image – surface of polyelectrolyte multilayer (PEI/PSS)₄ is shown in Fig. 2A.

Fig. 2 (A) (PEI/PSS)₄ polyelectrolyte assembly relief (AFM); (B) ZnO microcrystals relief (SEM); (C) Grand Canyon relief (SRTM). Synthetic reliefs with Hurst exponent (D) H = 2.5; (E) H = 3.5.

Scanning electron microscopy datasets and measurements

SEM imaging was performed using a VEGA 3 scanning electron microscope (TESCAN, Czech Republic). Secondary electron emission imaging was employed to analyze the surface morphology of the samples. For SEM sample preparation, a 20 μL particles water suspension was deposited onto a silicon substrate using a pipette dispenser and then dried under standard conditions. For ZnO particles, imaging was performed at a working distance of 4.5 mm, an accelerating voltage of 15 kV, and a magnification of 15.1k×. For the NiO particles, the working distance was 3.71 mm, the accelerating voltage 22 kV, and the magnification 22.2k×. Similarly, for the CuO particles, the working distance was set to 3.9 mm, the accelerating voltage to 22 kV, and the magnification to 8.01k×. For the MXene particles, the working distance was 10.28 mm, with an accelerating voltage of 5 kV and a magnification of 2.18k×. These imaging conditions were selected to optimize the visualization of the surface features of the respective particle types while maintaining clarity and minimizing charging artifacts. No additional filtering was applied to the SEM data.

The hydrothermal synthesis of zinc oxide (ZnO), nickel oxide (NiO), and copper oxide (CuO) microcrystals was carried out using a similar procedure for all three metal oxides microcrystals. The alkali solution for each metal was prepared by dissolving 1.487 g of zinc nitrate hexahydrate [Zn(NO₃)₂·6H₂O], 1.454 g of nickel nitrate hexahydrate[Ni(NO₃)₂·6H₂O], or 1.347 g of copper nitrate tetrahydrate [Cu(NO₃)₂·3H₂O] in deionized water, and then adding 4 g of NaOH on fill deionized water to form a 10 mL solution with a molar concentration of 0.005 M for each metal precursor. For each solution, the molar ratio of metal ions to hydroxide ions was maintained at 1 : 20, with [Zn²⁺], [Ni²⁺], or [Cu²⁺] = 5 mM and [OH⁻] = 0.1 M. Next, 3 mL of the prepared precursor solution was mixed with 6 mL of ethylenediamine (C₂H₄(NH₂)₂), 5 mL of deionized H₂O, and 30 mL of ethyl alcohol (C₂H₅OH). The mixture was then pretreated under an ultrasonic bath for 40 min to ensure uniform mixing and dispersion. Afterward, the solution was transferred to a Teflon-lined autoclave, and the hydrothermal synthesis was conducted at 180 °C for 20 h in an electric oven. Following the reaction, the crystalline products were washed five times with deionized water in a centrifuge to achieve a final pH of 7. The study also utilized MXene particles, which were synthesized following the method described in the relevant publication.⁶⁹ SEM dataset and statistics provided in Fig. S5–S8 and Table S2. One example of SEM image is shown in Fig. 2B and will be referred to hereinafter as ZnO.

NASA's Shuttle Radar Topography Mission datasets

Geodetic data for Earth's topography were obtained from NASA's Shuttle Radar Topography Mission (SRTM) using interferometric synthetic aperture radar (InSAR) technology utilizing publicly available datasets. The data were selected with consideration of Earth's curvature, particularly the variation in longitude degree length at different latitudes, to ensure accurate topographic representation. Four distinct datasets were selected for analysis, each representing a specific geographic region. The first dataset corresponds to Mount Ararat, with latitudinal coordinates ranging from 39.526° to 39.795°N, and longitudinal coordinates from 44.206° to 44.556°E. The second dataset pertains to Mount Fuji, with latitudes between 35.229° and 35.499°N, and longitudes between 138.568° and 138.898°E. The third dataset focuses on the Karelia region, with latitudes ranging from 61.812° to 62.081°N, and longitudes from 32.308° to 32.881°E. Finally, the fourth dataset covers the section of the Grand Canyon, with latitudinal coordinates between 36.268° and 36.537°N, and longitudinal coordinates from 112.705° to 112.371°W. Each dataset is represented as a 964 × 964 square matrix, covering an area of 30 × 30 km with ≈31 m per pixel line resolution, thus providing high-resolution topographic data for each region.

To avoid cluttering up the main text of the article, it provides examples of results for polyionic assembly PEI/PSS No. 1, Grand Canyon region, SEM of ZnO mycrocrystals, and two synthetic maps of height with Hurst exponent H = 2.5 and 3.5. The remaining results are presented in the SI. A standard description of surfaces in terms of roughness, typical linear surface profiles, and histograms of height distribution are also provided in the SI. SRTM dataset and statistic provided in Fig. S9–S12 and Table S3. One example of image is shown in Fig. 2C.

Synthetic generation of fractal-like random isotropic surfaces

To model the stochastic roughness of surfaces across multiple spatial scales, we implemented a method for generating synthetic surfaces based on spectral synthesis with controlled statistical and spectral properties. The generated surfaces emulate self-affine isotropic topographies characterized by power-law spectral decay and prescribed height statistics.⁶⁶

The procedure proceeds as follows. First, two independent random matrices of size (2N + 1) × (2N + 1) are generated: a Gaussian noise matrix G_i,j sampled from a normal distribution N(μ, σ), and a uniform noise matrix H_i,j sampled from a uniform distribution U(a, b). These matrices are used to modulate the amplitude and phase components of a synthetic surface in the frequency domain. Each spatial frequency component (m, n) ∈ [−N, N]² is assigned a complex-valued spectral coefficient according to:

where H > 0 is the prescribed spectral exponent controlling the decay of the power spectrum (related to the Hurst exponent via H = 3 − D, with D being the fractal dimension), and θ^m,n ∼ U(0,1) is an independent random phase. An inverse two-dimensional discrete Fourier transform is then applied to the matrix of coefficients C_m,n, yielding a complex-valued spatial field. The absolute value of this field defines the surface height map:

To ensure that the surface roughness matches a prescribed standard deviation σ_target, the resulting surface is normalized as:

z′(x,y) = z(x,y)·σ_target/σ(z)

where σ(z) denotes the standard deviation of the generated surface. This method produces statistically isotropic, self-affine rough surfaces with controlled spectral content and height distribution. It is particularly suited for numerical studies of contact mechanics, multiscale roughness effects, and stochastic tribological phenomena.

To evaluate the behavior of the proposed method across different spectral regimes, we generated a series of synthetic surfaces with systematically varied spectral exponent values. Specifically, a total of 10 surfaces were synthesized with spectral exponents ranging from 0.5 to 5.0 in increments of 0.5. The grid size was set to 128, resulting in square matrices of size 256 × 256 due to internal zero-padding and mirroring required by the Fourier-based generation method. All surfaces were generated with a fixed random seed for reproducibility. A Gaussian phase noise with standard deviation σ = 0.1 was applied. To introduce phase irregularity, a secondary uniform noise matrix h ∼ U(0,1), which modulated the phase of the Fourier coefficients. This enhanced the stochastic diversity in the synthesized surfaces while maintaining the target spectral envelope. All surfaces were normalized to have zero mean and unit standard deviation, ensuring consistency in amplitude scaling across the spectral series. These synthetic surfaces serve as a controlled testbed for evaluating metric sensitivity to spectral content.

Synthetic surface topographies were generated on square grids of size n × n (n = 256) in arbitrary spatial units. The surfaces are defined as two-dimensional scalar fields representing height values over a discretized domain with unit spacing in both lateral directions (Δx = Δy = 1 [a.u.]). No physical scale was imposed, allowing flexible rescaling depending on the application context (e.g., contact simulations, fractal analysis). The vertical scale (height) was normalized to ensure a prescribed root-mean-square roughness, also expressed in arbitrary units. Generated dataset and statistic provided in Fig. S13 and Table S4.

Fig. 2 shows the surface reliefs, the results of applying the software to which will be discussed in the main text of the article: samples of the (PEI/PSS)₄ polyelectrolyte assembly area #2 (AFM), ZnO microcrystallites (SEM), and the Grand Canyon (SRTM), and synthetic data with spectral exponent 2.5 (Fig. 2D) and 3.5 (Fig. 2E). The surfaces have significant differences between themselves in their morphology and texture (Fig. S1, S5, S9, S13 and Tables S1–S4), and also differ significantly in roughness parameters: deviation from the mean plane, standard deviation, and the difference in height between the deepest depression and the highest protrusion (roughness parameters are presented in Tables S1–S3).

Computational methods

In this study, our objective is to systematically compute the changes in the contact geometry between two surfaces, as illustrated in Fig. 3A.

Fig. 3 Schemes of the model, thought experiment: (A) the mechanism of contact between two textured surfaces under the application of normal and tangential forces, illustrating the conditions of static or sliding friction; (B) a conceptual thought experiment modeling the contact between a textured, absolutely ductile body and an ideally smooth, absolutely rigid body; (C) and (D) the principle of algorithm operation demonstrated using a surface relief height map: dividing the height map with an equidistant plane and binarizing regions below the plane to zero and above the plane to one.

Within this contact configuration, four distinct forces act on the system: two normal forces compressing the bodies against each other, and two tangential forces contributing to the initiation of sliding motion. To simplify the model, we assume that the surface under study is an absolutely ductile body, while the second object is completely rigid and has a perfectly smooth surface (Fig. 3B). It is important to note that our model does not explicitly address the displacement or trajectory of material resulting from the removal of surface layers along the height axis. Instead, our primary focus is on analyzing the evolution of the contact shape and topology, as well as developing software to simulate this process using relief data of any shape.

The determination of the shape and type of contact areas is performed by constructing equidistant planes along the height axis (the step size between these planes can vary, but in this study, we use 1000 equidistant planes between the maximum and minimum relief heights for all surfaces). The distance between the equidistant planes in relative units or # of threshold plane – thousandths of the maximum height of the relief – will later be used in the work as a value on which the number of objects of the contact pads, their topological and fractal properties depend. At each step, the portion of the surface relief above the plane is considered the contact area, as illustrated in Fig. 3C. This process involves binarizing the surface into a binary image where “1” represents contact pads, and “0” represents areas without contact. The method is demonstrated in Fig. 3C using an example of AFM data for a (PEI/PSS)₄ suraface.

To obtain a set of slices, raw data, which is presented as square matrices, is processed into arrays of slice values ranging from 0 to 1 with a given step. Each resulting slice is binarized using the Bradley-Roth algorithm,⁷³ chosen for its ability to smooth noise while preserving surface details and for its adaptability to any values in the height matrix. This algorithm meets our requirements, as the binarization boundary is defined following the Wellner scheme.⁷⁴

The characteristics calculated using the proposed method (Fig. 4) include the number and area of first-order lands (black pixels, or “0” against a background of white pixels, or “1”) and holes (white pixels, or “1”, against a background of black pixels, or “0”). Second-order lands (black “0” on a white “1” background over black “0”) and holes (white “1” on a black “0” background over white “1”) are also identified and analyzed. The method calculates additional metrics for these contact structures, including the Hausdorff dimension and lacunarity, both of which can be computed using box-counting techniques. Fixed grid and sliding box techniques to calculate Hausdorff dimension and lacunarity available to use through presenting in this paper software.

Fig. 4 Scheme of 1st and 2nd orders lands and holes (A) distribution of fractal dimensions on slices, binarized distribution of lands and holes on various slides (% of maximum height) (B) and principle of “lands” and “holes” identification in the program of data processing (C).

In this study, calculations are performed for contact pads of the second degree of nesting. For some surfaces, contact pads with higher degrees of nesting are observed, but their influence on the transition from rest to sliding is assumed to be negligible. The number and area of contact pads can be computed using either the 4-connected method or the 8-connected method, also referred to as the von Neumann and Moore neighborhoods, respectively, in cellular automata theory. The 4-connected neighborhood considers only those pixels connected by their sides as a single contact pad, whereas the 8-connected neighborhood includes pixels connected by both sides and corners as part of the same object. The paper presents the results of calculations using von Neumann connectivity. Calculations using various types of connectivity are also available thanks to the software provided in the work.

A detailed description of the developed algorithm is provided in the SI as pseudocode. The full implementation, released under the MIT license, is accessible via the https://github.com/ShockOfWave/Fractal-Analisys.

Results and discussion

Conceptually, the cut-off threshold method emulates the emergence of submerged topography during water recession: as the virtual ‘water level’ progressively lowers (decreasing threshold plane), isolated contact “islands” (binary regions) nucleate, expand, and eventually coalesce – revealing the intrinsic hierarchical landscape of the surface. Analogous to post-diluvial terrain exposure, our threshold-driven analysis exposes nested topological features across scales: each reduction of the cutting plane “drains” the computational domain, unmasking self-similar contact structures as latent morphology becomes discretely measurable.

Initially, the model was designed specifically for AFM data analysis. However, during validation, we identified its potential applicability to a broader range of imaging techniques, including SEM and SRTM. Remarkably, the numerical and topological dependencies on the cutting plane number exhibited consistent patterns across all datasets – including synthetic data – despite their distinct physical origins and scale. This uniformity suggests that the model captures fundamental numerical and topological invariants.

Numerical features, quantity of 1st and 2nd holes and lands vs. # of threshold plane shows in Fig. 5A–D correspondingly (in semi-logarithmic scale). The relationship between numerical features – specifically, the counts of 1st and 2nd-order holes and lands – and the threshold plane number is illustrated in Fig. 6 (Fig. 6A and C depict the trends for holes, while Fig. 6B and D present the corresponding dependencies for lands).

Fig. 5 Number of objects for binarized images on number of the threshold plane for: (A) 1st order holes; (B) 2nd order holes; (C) 1st order lands; (D) 2nd order lands; for: (PEI/PSS)4 #1 (AFM), ZnO (SEM), Grand Canyon (SRTM), and two synthetic reliefs with H = 2.5 and 3.5.

Fig. 6 (A) 1st order holes area; (B) 2nd order holes area; (C) ratio of the areas of the 1st to 2nd holes; (D) ratio derivative; for: (PEI/PSS)₄ #1 (AFM), ZnO (SEM), Grand Canyon (SRTM), and two synthetic reliefs with H = 2.5 and 3.5.

Visual analysis of the dependencies shows their qualitative similarity for all the studied relief types. The number of holes of the first kind demonstrates an increase in the range up to 350 cutting plane numbers, after which a decrease is observed: for SEM up to values of about 10⁴ holes and lands counts, for synthetic H = 3.5 up to ∼10³, for H = 2.5 and AFM up to 10¹, for SRTM up to 10⁰. Characteristic areas of determining the dependencies start from 0.8 for SEM, from 0.13 for synthetic H = 2.5, from 0.21 for AFM and H = 3.5, and from 0.24 for SRTM. Similar scales of values are characteristic of other numerical parameters. Of particular interest is the appearance of a significant number of holes of the second kind in a narrow range between 250 and 400, indicating a critical transition in topological properties. We hypothesize that this range corresponds to the stick-slip transition in our friction model, where connected slip pathways percolate through the contact interface. Notably, the synthetic data exhibit distinct behavioral patterns depending on the roughness parameter H. A clear correlation emerges between surface types: synthetic relief with H = 3.5 closely mirrors SEM data in its topological dependencies, while H = 2.5 shows greater affinity with AFM measurements. This divergence suggests that the Hurst exponent H – a critical parameter governing surface roughness – plays a pivotal role in determining the topological characteristics of these structures.

Notably, the synthetic data exhibit distinct behavioral patterns depending on the roughness parameter H. A clear correlation emerges between surface types: synthetic relief with H = 3.5 closely mirrors SEM data in its topological dependencies, while H = 2.5 shows greater affinity with AFM measurements. This divergence suggests that H – a critical parameter governing surface roughness – plays a pivotal role in determining the topological characteristics of these structures. These findings underscore the sensitivity of synthetic models to input parameters. Even modest variations in H induce substantial topological shifts, effectively transitioning the system's behavior between regimes characteristic of different experimental techniques.

This analysis also reveals an implicit normalization within the algorithm. While SEM data lack explicit height units (representing secondary electron detector signals), synthetic data are dimensionless in lateral measurements. In contrast, AFM and SRTM data are defined in both height and lateral dimensions (nm μm⁻¹ and m km⁻¹, respectively), preserving a characteristic aspect ratio of ∼10³ between lateral and vertical scales. Despite these dimensional differences, all numerical metrics exhibit similar behavioral patterns. However, a critical distinction emerges for SRTM results: their valid range is limited to # threshold planes 222–575. This limited operational range suggests noize and uneven distribution of brightness across the image. The datasets exhibit significant variation in their inherent spatial dimensionality. AFM data consists of 256 × 256-point height matrices, while SEM provides higher resolution 1024 × 1024-point matrices. SRTM data shows an intermediate dimensionality of 964 × 964 matrices, with synthetic surface reliefs generated at 256 × 256 matrices resolution. Remarkably, despite these substantial differences in initial data resolution, our unified processing algorithm demonstrates robust scale invariance. This is evidenced by both the consistent behavior of derived numerical metrics and the preservation of fundamental topological patterns across all datasets. Quantitatively, the maximum number of first-order topological features is inherently constrained by matrix dimensions. For 4-connectivity analysis, each isolated “island” requires a minimum 3 × 3-pixel domain (9 pixels), theoretically limiting feature counts to ∼N²/9, where N is the linear matrix dimension. The method's particular effectiveness in handling the non-standard 964 × 964 points SRTM data – while maintaining complete consistency with results from other techniques – strongly suggests its capacity to extract truly invariant topological features, independent of both spatial scale considerations and specific measurement methodology artifacts. Quantitative analysis confirms that observed first-order feature counts (Table 1) are orders of magnitude below theoretical maxima derived from matrix dimensionality. For 4-connectivity processing, the upper bound is constrained to ∼N²/9 where N = linear pixel dimension, representing the minimum 3 × 3-pixel domain required for isolated feature identification for 4-connected method. Crucially, our empirical results demonstrate profound margin between observed and theoretical limits. These observations collectively reinforce our algorithm's capacity to extract genuine multiscale invariants beyond dimensional artifacts. Crucially, while grid dimensions constrain object counts, they also influence area metrics. At critical threshold planes (#300–400), aggregated areas of first-order holes exhibit resolution-dependent saturation: SRTM/SEM datasets plateau near 10⁵ pixels², AFM/synthetic surfaces (both H = 2.5 and H = 3.5) converge to ∼10⁴ pixels². This divergence confirms that absolute geometric properties remain scale-bound, whereas topological descriptors.

Table 1 Observed and theoretical maxima of first-order features

Dataset	Matrix size	Theoretical 1st order holes max	Observed 1st order holes max	Margin factor ×
AFM	256 × 256	7.280	∼10^1	∼10^2
SEM	1024 × 1024	116.508	∼10^4	∼10^1
SRTM	964 × 964	103.298	∼10^0	∼10^5
Synthetic	256 × 256	7.280	∼10^3 (H = 3.5)	∼10^1
			∼10^2 (H = 2.5)	∼10^2

---

Analysis shows that there is absence of grid saturation artifacts: no dataset approaches pixel-density limits, validating topological significance over numerical constraints. We observe scale-independent object density – feature counts reflect intrinsic surface complexity rather than measurement resolution (e.g., SEM's 16× higher resolution than AFM yields only 1000× more features). SRTM's results exceptional behavior – Geological surfaces exhibit sparse feature distributions (≤1 object), contrasting with synthetic/experimental data, yet remain fully compatible within the method's operational framework. Notably, synthetic surfaces with H = 3.5 show 100× higher feature density than AFM despite identical matrix size, confirming Hurst exponent's role in governing topological complexity. These observations collectively reinforce our algorithm's capacity to extract genuine multiscale invariants beyond dimensional artifacts.

These observations collectively reinforce our algorithm's capacity to extract genuine multiscale invariants beyond dimensional artifacts. Crucially, while grid dimensions constrain object counts, they also influence area metrics. At critical threshold planes (#300–400), aggregated areas of first-order holes exhibit resolution-dependent saturation: SRTM/SEM datasets plateau near 10⁵ pixels², AFM/synthetic surfaces (both H = 2.5 and H = 3.5) converge to ∼10⁴ pixels². This divergence confirms that absolute geometric properties remain scale-bound, whereas topological descriptors (connectivity, nesting hierarchy) maintain invariance across resolutions.

In this study, the following interpretation of the Pearson correlation coefficient is adopted: values above 0.9 are considered to indicate a high correlation, in the range of 0.7–0.9 – a significant correlation, while values below or equal to 0.5 indicate the absence of a significant correlation. Fig. 7 presents a comparative analysis of the correlation dependencies between the quantitative characteristics of 1st and 2nd holes and islands. The results demonstrate significant correlation relationships between different types of data. For the first-order holes, a high correlation is observed between the PEM (AFM) and ZnO (SEM) data with a coefficient of 0.84, as well as between SEM (ZnO) and the synthetic relief H = 2.5 (0.87). The correlation between the Canyon and synthetic relief with H = 3.5 is 0.84. In the case of second-order holes, the correlation between SEM and synthetic relief H = 2.5, as well as between Canyon and synthetic relief H = 3.5, reaches 0.68, which is close to the threshold of significant correlation (0.7). The analysis of the 1st-order lands reveals even more pronounced correlation relationships: the AFM, SEM, and SRTM data demonstrate a strong correlation with the synthetic surfaces. In particular, the correlation between SEM, SRTM, and the synthetic relief H = 2.5 and H = 3.5 synthetic surfaces approach high values. Of particular note is the exceptionally strong correlation (0.96) between the SRTM data and synthetic H = 3.5, while the correlation with H = 2.5 is significantly lower (0.63). These results indicate that the degree of correlation significantly depends on both the type of topological objects analyzed (holes/lands), the order in which they are considered, and the parameters of the synthetic surfaces. The observed patterns may reflect fundamental features of the formation of surface structures with different research methods and synthesis parameters.

Fig. 7 Cross-correlation between: (A) 1st holes; (B) 2nd holes; (C) 1st lands; (D) 2nd lands; for: PEM (AFM), ZnO (SEM), grand canyon (SRTM), and two synthetic reliefs with H = 2.5 and 3.5.

To quantitatively assess the similarity in morphological evolution across different surface types, pairwise Pearson correlation coefficients were computed for topological and geometrical descriptors as functions of threshold level. The examined variables included binary indicators (1st holes, 2nd holes, 1st lands, 2nd lands), the areas of black regions on a white background (other holes) and white regions on a black background (inner holes), normalized ratios of island count to landmass (outer area/inner area), their derivatives with respect to the threshold level (outer/inner derivative), as well as the Hausdorff dimension at each level and lacunarity (Fig. 8 and 9C).

Fig. 8 Cross-correlation between: (A) 1st holes area; (B) 2nd holes area; (C) 1st/2nd areas; (D) 1st/2nd areas deriative; for: PEM (AFM), ZnO (SEM), Grand Canyon (SRTM), and two synthetic reliefs with H = 2.5 and 3.5.

The correlation matrices revealed unexpected regularities: despite differences in material origin, imaging modality, and scale (AFM, SEM, SRTM, synthetic), many threshold-dependent variables exhibited high correlations across samples. In particular, surfaces derived from different physical processes (AFM, SEM, SRTM) often showed correlation values exceeding 0.9 for individual variables. This supports the hypothesis that threshold-induced morphological transitions are governed by scale-invariant principles that remain stable across at least six orders of spatial magnitude. For the variable holes (white regions on black background), a moderate level of correlation was observed across samples, with a mean value of 0.63 and a standard deviation of approximately 0.31. Within the AFM group, the metric remained consistently high, enabling its use within a single modality. However, consistency dropped sharply when comparing AFM, SEM, and SRTM. Negative correlations (as low as −0.11) between certain pairs indicate inverted behavior of the metric under threshold variation, likely due to differences in microrelief density or specific features of the binarization algorithm.

A similar structure was observed for the 2nd holes variable (nested white regions surrounded by black), although intergroup correlations were even lower (mean ∼0.53). This reflects the metric's sensitivity to nested topological features influenced by texture variations and noise in the data. For 1st lands, intra-group consistency remained high within the AFM and SEM groups but dropped when compared to synthetic surfaces, likely due to differing topological nature of outer contours. SRTM occupied an intermediate position between physical and synthetic samples. A similar structure was found for 2nd lands, albeit with less pronounced extreme values, suggesting this descriptor is relatively stable across scales. The inner and outer holes area variables showed high consistency within groups, especially in AFM and SRTM. Correlations between AFM and SEM were also relatively high, particularly for samples with pronounced multiscale structures. Exceptions were observed for synthetic surfaces, especially at low threshold levels, where sharp drops in consistency occurred.

The outer/inner area showed high reproducibility within the AFM. SEM and SRTM surfaces exhibited coherent trends with noticeable within-group differentiation. SRTM showed high correlation with AFM and SEM, confirming structural affinity between natural and synthetic morphologies in terms of component distribution. The average correlation for this metric was 0.75. The derivative of outer/inner area derivative with respect to threshold emphasized the gradients of change in the spatial ratio of islands to land. It proved especially sensitive to differences in morphological dynamics, with an average correlation of about 0.63. The least consistent pairs (e.g., Si with H = 0.5 and 1.5) highlighted sharp mismatches in early segmentation stages. At the same time, high correlations were observed between MXene and ZnO, and between Lakes and Canyon, confirming structural similarity in segmentation dynamics.

The Hausdorff dimension and lacunarity (Fig. 9), evaluated at each threshold level, demonstrated the highest degree of universality. Intra-group consistency reached 0.993 for SRTM surfaces, 0.990 for SRTM, 0.978 for AFM, and 0.973 for SEM. Inter-group correlations were particularly strong between AFM and SRTM (0.976), SEM and SRTM (0.965), and SPTM and SRTM (0.969). This indicates that fractal dimension may serve as a reliable, universal descriptor of structural complexity, invariant to imaging method and scale. Overall, the cross-correlation analysis demonstrated that with appropriate normalization and descriptor selection, it is possible to uncover universal morphological patterns that persist across imaging scales, modalities, and surface types. Of particular importance are fractal dimension and the derived ratio-based metrics, which exhibit the greatest resilience to methodological variation and can be recommended as reference features for structural classification of surfaces.

Fig. 9 (A) Hausdorff dimension; (B) Lacunarity; (C) cross-correlation between Hausdorff dimension.

Cumulative distribution plots are shown in Fig. 10. In order to confirm the uniformity of the distribution law of islands/holes on sections, a comparison of the topography of various data of micro- and macro-dimensional surfaces obtained by different methods (AFM, SEM, SRTM, synthetic data) was performed. Fig. S3, S7, S11 and S15 shows a comparison of 3 characteristic data for topographies of different scales obtained by different methods, which show a good match of the nature of the curve's distribution. It should be noted that the data collection method was different in nature (van der Waals forces – AFM, work function – SEM, height of point – SRTM), but led to the same distributions. It should be noted that in some specific cases (absolutely smooth inclined surfaces, local delta functions, etc.), rarely encountered in nature, this dependence may differ from the observed one. It is worth noting that the 2nd order lands distribution differs from the others and is more similar to the standard Gaussian. This is confirmed for different samples at different scales obtained by different types of measure techniques. AFM on the Hausdorff topology, lacunarity and holes/lands ratios differ greatly for the Si sample, which may be due to a larger number of islands and holes compared to other samples. According to SEM data, nickel and copper oxides have more “flattened” distributions, which can be explained by small deviations in the height of their relief. This is confirmed by the distribution of linear surface profiles. In this regard, smooth areas with a smaller number of holes/islands can be distinguished on these samples, which is reflected in the corresponding hole/island distribution curves.

Fig. 10 Cumulative distribution plots comparing different statistical distributions for the analysis of internal and external holes and lands. Each subplot corresponds to a specific distribution: (top left to bottom right) normal, lognormal, half-normal, inverse Gaussian, gamma, Chi-square, inverse Chi-square, F-ratio, and Weibull distributions. The curves represent the empirical data for all calculated objects.

According to the data on geographical objects, it is clear that for the Lakes and Canyon datasets, the distribution also has a “flattened” character, which is logically explained by the presence of relatively flat zones of the same height (plateau and flat lake zones). The analysis revealed that the distributions of the number of objects in equidistant relief planes exhibit general similarity and appear to follow a universal statistical law. Our attempts to identify this law showed that the most common distribution models do not adequately describe the observed dependencies. This suggests that the distributions may be complex, possibly representing linear combination of several simpler ones.

Universal looks of dependencies were unexpected, especially considering the consistent dependencies observed in the number of objects, regardless of the object of study, the instrumental method, or the scale range (from 10⁻⁹ up to 2000 m). We aim to demonstrate that these distributions share a common origin. Potential approaches to address this include: calculating PCA, cross 1-Wasserstein distances or cross Kullback–Leibler divergence between distributions. It can be theoretically described through reformulating the problem in probabilistic terms (e.g., “calculating the probability of N holes forming in a map at the next iteration” or analyzing the equivalence to “finding the number of zeros in a shifted random value”). If the original landscape's distribution function is symmetric, this symmetry should hold across all such cases. Other possibilities include applying standard tests to reduce the degrees of freedom in the distributions. Our team views the challenge of identifying the underlying distribution type as an exciting opportunity. However, resolving this issue lies beyond the scope of this study. We plan to address the universality of dependencies comprehensively in future research. It is also planned to present the influence of the choice of the number of cutoff planes, 4- and 8-connectivity, box- and sliding-box counting methods on the simulation results.

One of the most striking results of this study is the emergence of universal dependencies in topological descriptors across different length scales, from nanometer-scale AFM data to kilometer-scale geodetic measurements. This suggests that despite differences in material properties, external conditions, and measurement techniques, the evolution of contact areas follows similar geometric and topological rules. Such universality hints at an underlying statistical or physical principle governing the transition from static to sliding friction, reinforcing the idea that surface interactions can be understood through a topological lens rather than solely relying on classical mechanical models. Another important aspect of our findings is their potential to refine existing models of frictional contact. While conventional approaches describe friction in terms of asperity interactions and stress distributions, our analysis shows that topological complexity – quantified through measures such as Hausdorff dimension and lacunarity – may serve as an independent predictor of frictional behavior. This could lead to a paradigm shift in tribology, where topological invariants provide a more generalized framework for studying friction across different materials and scales. Future studies should explore whether these topological characteristics influence friction coefficients directly and how they can be incorporated into predictive models for practical applications.

Analysis of AFM and SEM data on the surfaces of nanostructured materials, as well as data SRTM, allowed us to gain a deeper understanding of the geometric and topological characteristics of the studied surfaces at different size scale levels. The results of our study demonstrated that the proposed metrics provide valuable information about the geometric and topological properties of the surfaces and the contact area. It is especially noteworthy that the dependences of these metrics on the degree of surface deformation remain universal in the range of scales from 10⁻⁹ to 10³ m and do not depend on the type of surface or the method of obtaining relief data. It was found that the studied surfaces have similar topological characteristics during the transition from static friction to sliding friction. Thus, we can conclude that these metrics carry a significant amount of information and can be useful for further research in the field of tribology and materials science. The findings of this study extend beyond traditional tribological analysis and offer significant insights into non-equilibrium physics. The observation of universal dependencies in topological metrics across scales suggests that frictional processes may be governed by fundamental principles akin to phase transitions and critical phenomena. In non-equilibrium systems, where fluctuations and transient states play a crucial role, the ability to characterize contact evolution with robust topological measures could help elucidate the dynamics of self-organizing systems. This framework not only contributes to a deeper understanding of the stick-slip behaviour but also establishes a bridge between frictional instabilities and broader non-equilibrium processes observed in complex systems. Furthermore, the application of our threshold model combined with topological data analysis has far-reaching implications for materials science and surface interactions. By revealing universal characteristics across disparate scales – from nanostructured films to geological formations–this approach provides a new tool for predicting and tailoring surface properties. Such insights are invaluable in the design of advanced materials where controlled friction and adhesion are critical, such as in microelectromechanical systems and nanotechnology. Ultimately, this work lays the groundwork for developing a unified theoretical framework that can inform the engineering of surfaces with optimized performance in a wide range of applications, thereby reinforcing the importance of interdisciplinary research in advancing both fundamental science and technological innovation.

Conclusions

The model presented in this work outperforms classical models (e.g., Greenwood–Williamson (GW)) by integrating topological data analysis (TDA) and machine learning, eliminating reliance on idealized asperity assumptions. The proposed model demonstrates consistent topological and numerical metrics across scales from nanometers (AFM) to meters (SRTM), revealing universal trends in contact area evolution. Key metrics include synthetic data with H = 2.5 closely matched AFM results, while H = 3.5 aligned with SEM data, highlighting its role in predicting surface behavior, and high cross-correlation coefficients (>0.9) in Hausdorff dimension and lacunarity across diverse datasets confirm scale-independent fractal properties.

The model identifies critical thresholds (e.g., 250–400 # of threshold planes) where second-order contact areas emerge, signaling percolation-like transitions from static to sliding friction. Fractal dimensions are ranged from 1.5 to 2.0, with lacunarity values varying by surface type, quantifying multiscale roughness. The method's invariance to measurement techniques (AFM/SEM/SRTM) and resolution (256 × 256 to 964 × 964 pixels) enables broad applications in tribology, materials science, and geophysics. This framework bridges gap between theory and real-world surface interactions, offering predictive power for frictional behaviour across scales.

Thus our work advances beyond prior efforts by introducing a universal topological framework that quantifies surface interactions without empirical calibration, validated across 6 orders of magnitude (nm–m). Unlike traditional methods reliant on idealized asperities (e.g., Greenwood–Williamson), our model captures emergent multiscale behaviors through fractal and persistence metrics, akin to how crumpled GO membranes¹ exploit anisotropic nanoarchitecture to defy permeability-selectivity trade-offs. The innovation lies in (1) identifying scale-invariant descriptors (e.g., Hausdorff dimension) and (2) linking stick-slip transitions to percolation thresholds via topological criticality. This paradigm shift enables predictive design of surfaces for applications ranging from MEMS to geophysics, aligning with the broader materials science movement where tailored nanostructures (e.g., GO wrinkles, fractal contacts) enable unprecedented functionality.

Author contributions

ASA and MVZ: conceptualization, methodology, investigation (experiments and measurements), formal analysis, validation, visualization, writing – original draft. TAA and VIM: software development. PVG: data curation. DAK, MN, DVA, EVS: supervision, funding acquisition, resources, project administration, writing – review & editing. The authors read and approved the final manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

Geodetic data for Earth's topography were obtained from NASA's Shuttle Radar Topography Mission (SRTM) using interferometric synthetic aperture radar (InSAR) technology and are publicly available.

All datasets (AFM, SEM, SRTM, synthetic) are availible at https://github.com/aleks269/mathor_datasets/

The software developed for this study is available under the MIT license and can be accessed at https://github.com/ShockOfWave/Fractal-Analisys. Random surface synthesator developed for this study is available under the MIT licence and can be accessed at https://github.com/aleks269/rand_surf_gen_wm

Supplementary information includes AFM, SEM, SRTM, and synthetic relief images, histograms, typical linear roughness profiles, roughness statistics, results of the proposed method, pairwise cross-correlations across the entire dataset, and pseudocode listings of the software. See DOI: https://doi.org/10.1039/d5mh01327e

Acknowledgements

The research was carried out within the project No. FSER-2025-0026. The authors express their deepest gratitude to Dr. Elizaveta Anastasova for her invaluable contribution to this work in part of assistance with transition metal oxide particles synthesis and for conducting the study of these particles using SEM, which greatly enriched our study. DVA acknowledges the infrastructural support from the Ministry of Education (Singapore) under the Research Centre of Excellence program (grant EDUN C-33-18-279-V12, the Institute for Functional Intelligent Materials (I-FIM)).

Notes and references

1. P. Zhang, Q. Wang, Y. Zhang, M. Lin, X. Zhou, A. David, A. Ustyuzhanin, M. Chen, M. I. Katsnelson, M. Trubyanov, K. S. Novoselov and D. V. Andreeva, Strain-induced crumpling of graphene oxide lamellas to achieve fast and selective transport of H₂ and CO₂, Nat. Nanotechnol., 2025 DOI:10.1038/s41565-025-01971-8.
2. B. Feeny, A. Guran, N. Hinrichs and K. Popp, A Historical Review on Dry Friction and Stick-Slip Phenomena, Appl. Mech. Rev., 1998, 51(5), 321.
3. G. G. Adams, Self-excited oscillations of two elastic half-spaces sliding with a constant coefficient of friction, J. Appl. Mech., 1995, 62, 867–872.
4. M. Nosonovsky and V. Mortazavi, Friction-induced vibrations and self-organization: mechanics and non-equilibrium thermodynamics of sliding contact, CRC Press/Taylor & Francis, 2013.
5. H. Hertz, Über die Berührung fester elastischer Körper, J. Reine Angew. Math., 1881, 92, 156.
6. J. A. Greenwood and J. B. P. Williamson, Contact of nominally flat surfaces, Proc. R. Soc. London, Ser. A, 1966, 295(1442), 300–319.
7. D. J. Whitehouse and J. F. Archard, The properties of random surfaces of significance in their contact, Proc. R. Soc. London, Ser. A, 1970, 316, 97–121.
8. A. W. Bush, R. D. Gibson and T. R. Thomas, The elastic contact of a rough surface, Wear, 1975, 35, 87–111.
9. J. A. Greenwood and J. H. Tripp, The elastic contact of rough spheres, J. Appl. Mech., 1967, 34, 153–159.
10. J. A. Greenwood and J. H. Tripp, The contact of two nominally flat rough surfaces, Proc. Inst. Mech. Eng., 1970, 185(1), 625–633.
11. B. N. J. Persson, Elastoplastic Contact between Randomly Rough Surfaces, Phys. Rev. Lett., 2001, 87, 116101.
12. B. N. J. Persson, Contact Mechanics for Randomly Rough Surfaces, Surf. Sci. Rep., 2006, 61, 201–227.
13. J. Joe, M. Scaraggi and J. Barber, Effect of fine-scale roughness on the tractions between contacting bodies, Tribol. Int., 2017, 111, 52–56.
14. J. A. Greenwood and J. J. Wu, Surface Roughness and Contact: An Apology, Meccanica, 2001, 36, 617–630.
15. J. R. Barber, Multiscale Surfaces and Amontons’ Law of Friction, Tribol. Lett., 2013, 49, 539–543.
16. L. Afferrante, G. Carbone and G. Demelio, Interacting and Coalescing Hertzian Aperities: A New Multiasperity Model, Wear, 2012, 278–279, 28–33.
17. J. F. Archard and W. Hirst, An Examination of a Mild Wear Process, Proc. R. Soc. London, Ser. A, 1957, 238, 515–528.
18. B. R. Mandelbrot, The Fractal Geometry of Nature, W.H. Freeman and Company, New York, NY, USA, 1977.
19. B. N. J. Persson, I. M. Sivebaek, V. N. Samoilov, K. Zhao, A. I. Volokitin and Z. Zhang, On the origin of Amonton's friction law, J. Phys.:Condens. Matter, 2008, 20(39), 395006.
20. R. I. Taylor, Rough Surface Contact Modelling–A Review, Lubricants, 2022, 10(5), 98.
21. J. Gao, W. D. Luedtke, D. Gourdon, M. Ruths, J. N. Israelachvili and U. Landman, Frictional Forces and Amonton's Law: From the Molecular to the Microscopic Scale, J. Phys. Chem. B, 2004, 108, 3410–3425.
22. A. Majumdar and B. Bhushan, Role of fractal geometry in roughness characterization and contact mechanics of surfaces, ASME J. Tribol., 1990, 112, 205–216.
23. J. F. Archard, Elastic Deformation and the Laws of Friction, Proc. R. Soc. London, Ser. A, 1957, 243, 190–205.
24. J. F. Archard, Single Contacts and Multiple Encounters, J. Appl. Phys., 1961, 32, 1420.
25. F. P. Bowden and D. Tabor, The Friction and Lubrication of Solids; Oxford Classics Series, Clarendon Press, Oxford, UK, 2001.
26. P. J. Blau, On the Nature of Running-In, Tribol. Int., 2008, 38, 1007–1012.
27. P. R. Nayak, Random Process Model of Rough Surfaces, J. Lubr. Technol., 1971, 93, 398–407.
28. P. R. Nayak, Some Aspects of Surface Roughness Measurement, Wear, 1973, 26, 165–174.
29. P. R. Nayak, Random Process Model of Rough Surfaces in Plastic Contact, Wear, 1973, 26, 305–333.
30. D. J. Whitehouse and J. F. Archard, The Properties of Random Surfaces of Significance in Their Contact, Proc. R. Soc. London, Ser. A, 1970, 316, 97–121.
31. W. R. Chang, I. Etsion and D. B. Bogy, An elastic-plastic model for the contact of rough surfaces, J. Tribol., 1987, 109(2), 257–263.
32. L. Kogut and I. Etsion, Elastic-plastic contact analysis of a sphere and a rigid flat, J. Appl. Mech., 2002, 69(5), 657–662.
33. B. B. Mandelbrot, The Fractal Geometry of Nature, W. H. Freeman and Company, New York, 1977.
34. R. S. Sayles and T. R. Thomas, Measurements of the Statistical Microgeometry of Engineering Surfaces, J. Lubr. Technol., 1979, 101, 409–417.
35. R. S. Sayles and T. R. Thomas, Surface Topography as a Nonstationary Random Process, Nature, 1978, 271, 431–434.
36. T. R. Thomas and R. S. Sayles, Some Problems in the Tribology of Rough Surfaces, Tribol. Int., 1978, 11, 163–168.
37. A. Majumdar and B. Bhushan, Role of Fractal Geometry in Roughness Characterization and Contact Mechanics of Surfaces, ASME J. Tribol., 1991, 112, 205–261.
38. A. Majumdar and B. Bhushan, Fractal Model of Elastic-Plastic Contact between Rough Surfaces, ASME J. Tribol., 1991, 113, 1–11.
39. R. L. Jackson, An Analytical Solution to an Archard-Type Fractal Rough Surface Contact Model, Tribol. Trans., 2010, 53, 543–553.
40. A. Gujrati, S. R. Khanal, L. Pastewka and T. D. B. Jacobs, Combining TEM, AFM, and Profilometry for Quantitative Topography Characterization across All Scales, ACS Appl. Mater. Interfaces, 2018, 10, 26169–29178.
41. T. D. B. Jacobs and L. Pastewka, Surface topography as a material parameter, MRS Bull., 2022, 47, 1205–1210.
42. M. Nosonovsky and G. G. Adams, Vibration and stability of frictional sliding of two elastic bodies with a wavy contact interface, J. Appl. Mech., 2004, 71, 154.
43. L. Chen, S. Zhang, Y. Duan, X. Song, M. Chang, W. Feng and Y. Chen, Silicon-containing nanomedicine and biomaterials: materials chemistry, multi-dimensional design, and biomedical application, Chem. Soc. Rev., 2024, 53, 1167.
44. D. Tian, Y. Song and L. Jiang, Patterning of controllable surface wettability for printing techniques, Chem. Soc. Rev., 2013, 42, 5184.
45. H. K. Raut, V. A. Ganesh, A. S. Nair and S. Ramakrishna, Anti-reflective coatings: A critical, in-depth review, Energy Environ. Sci., 2011, 4, 3779.
46. M. S. Hassan and M. Nosonovsky, Triboinformatics: machine learning algorithms and data topology methods for tribology, Surf. Innov., 2022, 10(45), 229–242.
47. N. Yin, Z. Xing, K. He and Z. Zhang, Tribo-informatics approaches in tribology research: A review, Friction, 2023, 11, 1–22.
48. J. F. Senge, A. Heydari Astaraee, P. Dłotko, S. Bagherifard and W. A. Bosbach, Extending conventional surface roughness ISO parameters using topological data analysis for shot peened surfaces, Sci. Rep., 2022, 12, 5538.
49. H. Canot, P. Durand, E. Frénod, B. Hassoune-Rhabbour, V. Nassiet and O. Tramis, Topological data analysis for roughness surfaces of bonding assembly, Discrete Contin. Dyn. Syst. - S, 2025, 18(7), 1743–1766.
50. M. Nosonovsky and A. S. Aglikov, Triboinformatics: Machine Learning Methods for Frictional Instabilities, Facta Univ. Ser. Mech. Eng., 2024, 22(3), 423.
51. M. S. Hasan and M. Nosonovsky, Topological data analysis for friction modeling, Europhys. Lett., 2021, 135(5), 56001.
52. M. V. Zhukov, M. S. Hassan, P. Nesterov, M. Sabbouh, O. Burdulenko, E. V. Skorb and M. Nosonovsky, Topological Data Analysis of Nanoscale Roughness in Brass Samples, ACS Appl. Mater. Interfaces, 2022, 14, 2351–2359.
53. A. T. S. Ireddy, F. D. E. Ghorabe, E. I. Shishatskaya, G. A. Ryltseva, A. E. Dudaev, D. A. Kozodaev, M. Nosonovsky, E. V. Skorb and P. S. Zun, Benchmarking Unsupervised Clustering Algorithms for Atomic Force Microscopy Data on Polyhydroxyalkanoate Films, ACS Omega, 2024, 9, 21595–21611.
54. I. D. Chetinel, A. A. Botnar, A. S. Novikov, E. A. Muraveva, A. T. S. Ireddy, P. S. Zun, S. E. Solovieva, I. S. Antipin, E. V. Skorb and A. A. Muravev, Colloid J., 2025, 87(2), 236–245.
55. A. A. Botnar, O. P. Novikov, O. A. Korepanov, E. A. Muraveva, D. A. Kozodaev, A. S. Novikov, M. Nosonovsky, E. V. Skorb and A. A. Muravev, Crystallization Control of Anionic Thiacalixarenes on Silicon Surface Coated with Cationic Poly(ethyleneimine), Langmuir, 2024, 40(46), 24634–24643.
56. S. Shityakov, A. S. Aglikov, E. V. Skorb and M. Nosonovsky, Voronoi Entropy as a Ligand Molecular Descriptor of Protein−Ligand Interactions, ACS Omega, 2023, 8, 46190–46196.
57. A. S. Aglikov, M. V. Zhukov, T. A. Aliev, D. A. Kozodaev, M. Nosonovsky and E. V. Skorb, New metrics for describing atomic force microscopy data of nanostructured surfaces through topological data analysis, Appl. Surf. Sci., 2023, 670, 160640.
58. A. I. Vakis, V. A. Yastrebov, J. Scheibert, C. Minfray, L. Nicola, D. Dini, A. Almqvist, M. Paggi, S. Lee, G. Limbert, J. F. Molinari, G. Anciaux, R. Aghababaei, S. E. Restrepo, A. Papangelo, A. Cammarata, P. Nicolini, C. Putignano, G. Carbone, S. Stupkiewicz and M. Ciavarella, Modeling and simulation in tribology across scales: an overview, Tribol. Int., 2018, 125, 169.
59. M. Urbakh, J. Klafter, D. Gourdon and J. Israelachvili, The nonlinear nature of friction, Nature, 2004, 430(6999), 525.
60. International Organization for Standardization. ISO Norm 25178: Geometrical Product Specifications (GPS) – Surface Texture: Areal. Switzerland, 2016.
61. International Organization for Standardization. ISO Norm 4287: Geometrical Product Specifications (GPS) – Surface Texture: Profile Method. Switzerland, 1997.
62. Japanese Standards Association. JIS B 0601 Geometrical Product Specifications (GPS) – Surface Texture: Profile Method–Terms, Definitions and Surface Texture Parameters. 2017.
63. German Institute for Standardization. DIN 4762:1989-01 Surface Roughness Terminology Surface & Its Parameters. 1989.
64. GOST 25142-82.
65. M. Nosonovsky and B. Bhushan, Energy transitions in superhydrophobicity: low adhesion, easy flow and bouncing, J. Phys.:Condens. Matter, 2008, 20(39), 395005.
66. D. Saupe, Algorithms for random fractals, The Science of Fractal Images, Springer, 1988, pp. 71–113.
67. M. V. Zhukov, A. S. Aglikov, M. Sabboukh, D. A. Kozodaev, T. A. Aliev, S. A. Ulasevich, M. Nosonovsky and E. V. Skorb, AFM-Topological Data Analysis of Brass after Ultrasonic Surface Modification, ACS Appl. Eng. Mater., 2023, 1(8), 2084.
68. M. Sychov, A. Eruzin, A. Semenova, P. Katashev, S. Mjakin, M. V. Zhukov, A. Aglikov, M. Nosonovsky and E. V. Skorb, Deposition of Nanostructured Tungsten Oxide Layers by a New Method: Periodic Modulation of the Deposition Angle, Langmuir, 2023, 39(35), 12336.
69. A. S. Aglikov, T. A. Aliev, M. V. Zhukov, A. A. Nikitina, E. Smirnov, D. A. Kozodaev, M. I. Nosonovsky and E. V. Skorb, Topological Data Analysis of Nanoscale Roughness of Layer-by-Layer Polyelectrolyte Samples Using Machine Learning, ACS Appl. Electron. Mater., 2023, 5(12), 6955.
70. F. D. E. Ghorabe, A. Aglikov, A. S. Novikov, M. Nosonovsky, G. A. Ryltseva, A. E. Dudaev, N. G. Menzianov, E. V. Skorb and E. I. Shishatskaya, Topography hierarchy of biocompatible polyhydroxyalkanoate film, RSC Adv., 2024, 14, 19603.
71. A. Zelenina, E. V. Skorb, D. V. Andreeva and N. Orekhov, Ultrafast water diffusion along the interface between oxidized and pristine regions in graphene oxide: Reactive molecular dynamics study, Comput. Mater. Sci., 2025, 247, 113461.
72. M. Chen, Q. Wang, M. Trubyanov, K. Yang, A. S. Aglikov, G. Qi, E. V. Skorb, K. S. Novoselov and D. V. Andreeva, Large-scale self-assembly of anisotropic graphene oxide films via blade coating: sustainable design and stimuli-responsive performance for biomimicry, Mater. Des., 2023, 233, 112205.
73. D. Bradley and G. Roth, Adaptive thresholding using the integral image, J. Graph. Tools, 2007, 12(2), 13.
74. P. D. Wellner, Adaptive Thresholding for the DigitalDesk, EuroPARC Technical Report, 1993, EPC-93-110, 404.

---

Footnote

† These authors contributed equally to this work.

---