# Electrical and Thermal Conductivity of Complex-Shaped Contact Spots

Paul Beguin, Vladislav Yastrebov

<figure><img src="image_1.png" /></figure>

## To cite this version:

Paul Beguin, Vladislav Yastrebov. Electrical and Thermal Conductivity of Complex-Shaped Contact Spots. 2024. hal-04800297

HAL Id: hal-04800297

https://hal.science/hal-04800297v1

Preprint submitted on 24 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Electrical and Thermal Conductivity of
Complex-Shaped Contact Spots

Paul Beguin & Vladislav A. Yastrebov

MINES Paris - PSL, Centre des Matériaux, CNRS UMR 7633, BP 87, 91003, Evry, France

Abstract

This paper explores the electrical and thermal conductivity of complex contact spots on the surface of a half-space. Employing an in-house Fast Boundary Element Method implementation, various complex geometries were studied. Our investigation begins with annulus contact spots to assess the impact of connectedness. We then study shape effects on "multi-petal" spots exhibiting dihedral symmetry, resembling flowers, stars, and gears. The analysis culminates with self-affine shapes, representing a multiscale generalization of the multi-petal forms. In each case, we introduce appropriate normalizations and develop phenomenological models. For multi-petal shapes, our model relies on a single geometric parameter: the normalized number of "petals". This approach inspired the form of the phenomenological model for self-affine spots, which maintains physical consistency and relies on four geometric characteristics: standard deviation, second spectral moment, Nayak parameter, and Hurst exponent. As a by product, these models enabled us to suggest flux estimations for an infinite number of petals and the fractal limit. This study represents an initial step into understanding the conductivity of complex contact interfaces, which commonly occur in the contact of rough surfaces.

**Keywords:** Conductivity, flower-shaped spots, self-affine spots, boundary element method, fractal limit

Supplementary material for this article is supplied as a separate archive available from Zenodo 10.5281/zenodo.10200997. The fast BEM code used for simulations is shared as open-source at github.com/vyastreb/HBEM and a snapshot used for all simulations in the paper is available at Software Heritage SWHID.

Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Methods</td><td>5</td></tr><tr><td>3</td><td>Conductivity of a not simply connected spot</td><td>7</td></tr><tr><td>4</td><td>Conductivity of flower-shaped spots</td><td>9</td></tr><tr><td>4.1</td><td>Conductivity results</td><td>11</td></tr><tr><td>4.2</td><td>Alternative "multi-petal" spots</td><td>14</td></tr><tr><td>5</td><td>Conductivity of self-affine spots</td><td>17</td></tr><tr><td>5.1</td><td>Geometry of spots</td><td>17</td></tr><tr><td>5.2</td><td>Geometrical characteristics</td><td>17</td></tr><tr><td>5.3</td><td>Normalization and result</td><td>20</td></tr><tr><td>5.4</td><td>Results with renormalized standard deviation</td><td>24</td></tr><tr><td>5.5</td><td>Conductivity model</td><td>26</td></tr><tr><td>5.5.1</td><td>Covariance matrix</td><td>26</td></tr><tr><td>5.5.2</td><td>Phenomenological model</td><td>27</td></tr></table>
<table><tr><td><b>6 Conclusion</b></td><td><b>29</b></td></tr><tr><td>6.1 Flower-shaped spots and other simple forms</td><td>29</td></tr><tr><td>6.2 Conductivity of self-affine spots</td><td>31</td></tr><tr><td>6.3 Contact spots between rough surfaces</td><td>31</td></tr><tr><td><b>A Geometrical characteristics of flower-, star- and gear-like shapes</b></td><td><b>32</b></td></tr><tr><td><b>B Geometrical characteristics of self-affine spots</b></td><td><b>32</b></td></tr><tr><td><b>C Physical consistency of the phenomenological model</b></td><td><b>34</b></td></tr></table>

# 1 Introduction

The study of mechanical contact plays a crucial role in numerous natural and engineering systems. Some of these systems involve sliding motion, frictional resistance, lubrication and wear, while others operate exclusively in normal contact. In the former case, the interplay of mechanical deformation, chemistry, and thermal effects due to frictional heat generation and diffusion significantly influence the overall behavior of the interface [Blok, 1963, Bowden et al., 2001, Rice, 2006]. This heat generation, coupled with intense shear deformation, can lead to phase transitions/transformation, recrystallization, and other metallurgical or chemical effects [Jacobson and Hogmark, 2009, Goldsby and Tullis, 2011, Yamashita et al., 2015], such as local welding/galling and abrasive or adhesive wear, depending on oxygen influx [Baydoun et al., 2020, Baydoun et al., 2022]. In the case of normal contact, frictional dissipation within the interface remains minimal. However, the nature of contact interfaces profoundly impacts their conductive properties for both heat and electric charge.

In both natural and engineering systems, all surfaces exhibit roughness. As a result, at relatively light contact loads, the true contact area – formed by several intimately contacting "asperities" – is smaller than the apparent contact area between solids. Apart from other mechanical and physical properties, the contact area fraction and its morphology also dictate energy transfer through the contact interface. Since the contact area fraction evolves under increasing load [Archard, 1953, Archard, 1957, Greenwood and Williamson, 1966], the interface conductivity depends on this load as well.

Surface roughness can be described by its height distribution or its moments. In the simplest case, surfaces can be described as Gaussian, but in practice, most surfaces in contact exhibit an asymmetric height distribution due to wear or residual plastic deformations induced during contact, particularly after the running-in process. An example of this asymmetry is asphalt concrete, which, due to manufacturing techniques, inherently has relatively flat plateaus and deep valleys. Although height distribution or first moments of roughness are relevant parameters for specific applications, they do not adequately describe mechanical contact between rough surfaces, where the curvature of contacting "asperities" is of primary importance. Consequently, power spectral density (PSD) and its first moments represent the key parameters in describing the elastic contact of rough surfaces [Greenwood and Williamson, 1966, Nayak, 1971, Persson, 2001, Yastrebov et al., 2017].

Aligned with PSD-based models which provide a multiscale representation of surface features, some authors proposed fractal models of surface roughness [Majumdar and Tien, 1990, Majumdar and Bhushan, 1991, Dodds and Robson, 1973], which, to some extent also ensure a suitable description of real rough surfaces in the framework of Archard's concept of "protuberances on protuberances" [Archard, 1953] however the mechanical consistency of such models dealing with elasto-plastic contact remains questionable.

In the study of rough surfaces in contact, one can adopt either statistical or deterministic approaches. The former, which includes multi-asperity models [Greenwood and Williamson, 1966, Bush et al., 1975, Greenwood, 2006, Carbone and Bottiglione, 2008] and Persson's model [Persson, 2001, Manners and Greenwood, 2006], deals with probability densities and
average statistical properties of roughness. This approach is applicable to contacts with a statistically meaningful number of contacting asperities.

These analytical and semi-analytical methods are built on certain assumptions and, as such, have limitations when it comes to predicting the rigorously formulated contact problem between rough surfaces [Manners and Greenwood, 2006, Carbone and Bottiglione, 2008, Yastrebov et al., 2017]. On the other hand, deterministic approaches, such as multi-asperity models with interactions [Afferrante et al., 2012, Yastrebov, 2019, e.g.] can effectively account for a small number of contact spots, their interaction and spatial heterogeneities in roughness. Different more advanced methods can be employed to address contact problems within a deterministic framework. The most versatile is the finite-element method [Hyun et al., 2004, Pei et al., 2005, Yastrebov et al., 2011, Gao and Bower, 2006], which can inherently handle nonlinear and heterogeneous material behavior and large deformations. However, it is computationally intensive and necessitates solving mechanical equations not only on the surface but also in the bulk. Another group of methods, including boundary-element and spectral methods, focuses exclusively on surface interactions and, in their basic implementation, relies on space-invariant fundamental solutions [Kalker, 1977, Stanley and Kato, 1997, Polonsky and Keer, 1999, Bonnet, 1999, Campaná and Müser, 2006]. Over the past few decades, these methods have been extended to tackle more complex problems involving heterogeneous and nonlinear materials [Putignano et al., 2015, Amuzuga et al., 2016, Frérot et al., 2019, Pérez-Ràfols et al., 2023].

When two conductive solids come into contact, a localized resistance emerges due to the discontinuous nature of the actual contact area. From a mathematical standpoint, at the macro-scale level, this phenomenon is characterized by a discontinuity in potential or temperature at the interface of contact. At smaller scales, as discussed in the pioneering works of Holm [Holm, 1957], the conductivity at contact spots exhibits a continuous change in potential (temperature) for electrical (and thermal) contacts. In a manner analogous to early models of mechanical contact involving rough surfaces, the true contact area can be represented by a collection of discrete circular contact spots distributed across the nominal contact area, and potentially extending even beyond it [Greenwood and Tripp, 1967, Greenwood, 1966]. The latter case is possible when there is a lack of scale separation between contacting shapes and roughness characteristics, see [Greenwood et al., 1984, Yastrebov, 2019].

One of numerous examples that accounts for inelastic deformations in contact is the work of Kogut and Komvopoulos [Kogut and Komvopoulos, 2003], where the authors used a fractal geometry with a simple overlap (cut-off) model to establish a theoretical connection between fractal roughness parameters and constriction resistance, assuming elasto-plastic deformation. However, it is worth recalling here, that overlap models were shown to produce erroneous results and should not be used for quantitative analysis [Pei et al., 2005, Greenwood, 2006, Dapp et al., 2012]. In the elastic regime, Barber [Barber, 2003] rigorously demonstrated the equivalence between electrical/thermal contact resistance and contact normal stiffness. Consequently, within the contact mechanics community, this equivalence is often employed as a justification for disregarding a separate study of thermoelectrical diffusion equation in mechanical contact of rough surfaces.

However, in real-life applications, electrical and thermal resistance are not merely reduced to contact stiffness. First, the contact-induced deformation is often accompanied by inelastic deformations. In addition, in case of thermal resistance, due to the presence of electrically insulating or weakly conducting oxide films and surface contamination, the conductive contact area is reduced [Holm, 1957], [Slade, 2017, Ch.1-4]. Furthermore, in the case of thermal conductivity, additional convective and radiative contributions to heat exchange cause the conductivity to deviate from the strict mathematical equivalence between elastic stiffness and thermal/electric resistance established by Barber [Barber, 2003]. Therefore, in all aforementioned contexts, a study of coupled thermo-/electromechanical problems is worth investigation.

The growing demand for micro-electric devices [Tu, 2003] has spurred an increasing need for *electrical* contact models at the microscale. The constriction resistance model remains consistent with experimental studies [Watanabe, 1986] and can be extended to various contact shapes [Nakamura, 1993, Sano, 1985] as well as to multi-spot contact configura-
tions [Thomas and Probert, 1970, Nakamura and Minowa, 1986]. Contact resistance is influenced by the hardness and resistivity of the contacting solids, but is also affected by oxidation at the contact surface [Cuthrell and Tipping, 1973]. However, the study of simple conductive effects reaches its limits at small scales, where electrical resistance must incorporate the ballistic description of electrons' motion, also known as Sharvin's resistance [Mikrajuddin et al., 1999, Jensen et al., 2005]. Similar research has been conducted for *thermal contact resistance*, which has numerous applications in aerospace, automotive and electronic domains, including conductivity of bolted joints and of thermal cooling devices [Kumano et al., 1994]. This issue is experimentally represented by the thermal contact resistance (TCR), which necessitates careful attention to the set-up precision [Kempers et al., 2009] in both steady-state and transient studies [Fieberg and Kneer, 2008, Burghold et al., 2015]. Experimental advancement was accompanied by theoretical [Cooper et al., 1969, Lambert and Fletcher, 1997, Zou et al., 2008, Jackson et al., 2008] and numerical investigations [Sadowski and Stupkiewicz, 2010, Murashov and Panin, 2015] of thermal resistance between rough surfaces in contact.

The primary motivation for this study stems from the observation of contact spot geometry formed between model rough surfaces (self-affine random geometry). Figure 1(a) illustrates how individual contact spots evolve under increasing load [Yastrebov et al., 2015], while in (b) three separate contact spots obtained from similar simulations demonstrate high complexity both in terms of connectedness and boundary shape. Even under relatively small loads, very complex contact spots can be formed if the spectral content of roughness [Nayak, 1971] is sufficiently rich. However, basic models of contact resistance assume that individual contact spots distributed over the nominal contact area possess simple shapes: elliptical or circular.

In contrast, the realistic contact shapes (Fig. 1) can be non-simply connected (having holes) and exhibit highly complex boundaries. This complexity can be characterized by the ratio of the square root of the area to its perimeter, also known as compactness. In this study, we investigate the impact of connectedness and compactness of individual contact spots on their conductive properties. Instead of studying contact spots obtained in direct numerical simulations or rough contact, as shown in Fig. 1, we construct relatively simple models that capture the primary features of such spots: (1) connectedness, and (2) low compactness. The first effect is represented by an annulus shape with a varying ratio of inner to outer radius. Compactness, as a first approximation, is modeled through a flower-shaped spot with varying numbers of petals and their lengths. Additionally, we explore "multi-petal" configurations to generalize our findings. The complexity of contact spot geometry is further addressed by modeling self-affine contact spots, akin to multiscale petals, paraphrasing Archard with "petals on top of petals on top of petals." The *primary objectives* of this study are twofold: firstly, to understand the subtle relationship between the geometry of such complex shapes and their thermal and electrical conductivity; and secondly, to construct a simplified model capable of predicting this conductivity based on a set of basic geometrical characteristics.

<figure><img src="image_3.png" /><figcaption>Figure 1: Simulation results of the true contact-area between self-affine random rough surfaces taken from [Yastrebov et al., 2015]: (a) evolution of separate contact spots under increasing load, (b) snapshots of separate clusters highlighting the complexity of the shape.</figcaption></figure>
This paper is structured as follows: Section 2 briefly introduces the numerical methods employed in this study. The subsequent three sections study the effects of the topology and the shape of contact spots on simple examples. Specifically, Section 3 explores the influence of non-simple connectedness using the flux through an annular spot as a case study. Section 4 examines the effect of compactness on conductivity of a simplified model of a flower-shaped spot, characterized by equally spaced, uniformly sized petals, and similar formations. Section 5 extends this examination to self-affine spots, which could be seen as flower-shaped spots with polydispersed petals following a self-affine distribution. The paper concludes with Section 6.

## 2 Methods

Electrostatics and stationary thermal conductivity are described by elliptic partial differential equation, also known as Poisson's equation $\Delta U = S$, where $U$ takes a role of electrostatic potential or temperature, $S$ is the source/sink term and the flux is given by $j = -k\nabla U$ with $k$ being the electric or thermal conductivity. Throughout the paper we assume linear isotropic, homogeneous and constant conductivity $k$, i.e. $\nabla k = 0$, $\partial k / \partial U = 0$.

In this study, to solve Poisson's equation (elliptic PDE) on an isotropic, linear and homogeneous half-space we used an in-house implementation of the boundary element method (BEM) [Bonnet, 1999] in its fast-BEM version based on hierarchical matrices [Grasedyck, 2005, Grasedyck and Hackbusch, 2003, Hackbusch, 2015, Chaillat et al., 2017]. The solution of the Poisson's equation can make use of a particular solution, formulated with a Green function $G(x, y)$, resulting in an integral equation involving only the contact area $A$ where a constant potential $U(x) = U_0$ is prescribed:

$$U(x) = \int_A \frac{j_n(y)}{k} G(|x - y|) dS_y \quad (1)$$

The Green's function depends only on the distance between the "source" point $y$ and the "observation" point $x$ $G(|x - y|)$ stands for the elementary solution of the Poisson's equation.

In this study we focus on a homogeneous Poisson's problem $\Delta U = 0$ with a constant potential $U_0$ prescribed on the contact area $A$ and zero flux $j_n = 0$ outside, i.e. on $\bar{A}$, see Fig. 2(a). This conductivity problem is equivalent to a mechanical contact problem of elastic indenting a half-space by a flat indenter with the identical section $A$ [Barber, 2003]. The only difference is that the flux $j_n$ is replaced by the contact pressure $p$ and the potential $U$ is replaced by the normal displacement $u_z$, which could be readily derived from the Boussinesq's solution [Boussinesq, 1885] as

$$u_z(x) = \int_A \frac{2(1 - v^2)p(y)}{E} G(|x - y|) dS_y \quad (2)$$

Here, we assume that the half-space is made of an isotropic material linearly elastic material with $E$ being the Young's modulus and $v$ being the Poisson's ratio. To go beyond linear contact problems, the analogy still persists between the normal contact stiffness, the derivative of the contact force $F$ to normal displacement $u_z$, and the electrical/thermal resistance $R$ as soon as the contact area is identical to the conducting area $A$:

$$R = \frac{E}{2k(1 - v^2)} \left[ \frac{\partial F}{\partial u_z} \right]^{-1} \quad (3)$$

For a further discussion and a rigorous derivative of this equivalence the reader is referred to [Barber, 2003].

The integral equation (Eq. 1) can be turned into a linear system to solve, involving the construction of a fully-populated matrix. Hierarchical matrices allow to overcome the drawbacks of the classical BEM, i.e. excessive storage and a tedious resolution of linear system of equations with a full matrix. The open-source Python code is shared under BSD3 license in [Beguin and Yastrebov, ] along with the details of its implementation and a short documentation. Constant interpolation triangular elements
<figure><img src="image_4.png" /><figcaption>Figure 2: Conductivity through a contact spot on the half-space: (a) – the problem set-up for an annulus with prescribed potential $U_0$ in the zone A and zero flux $j_n = 0$ outside, i.e. on $\bar{A}$; (b-c) – an illustration of the Richardson extrapolation technique: the same problem is solved on two meshes possessing different reference sizes $h_1, h_2$, the resulting fluxes $Q_1, Q_2$ are extrapolated to a limit value $h \to 0$ assuming a linear convergence with the mesh size $h$ using Eq. (4).</figcaption></figure>

were used. In addition, some problems were also solved using our in-house finite element suite Z-set [Besson and Foerch, 1997, Z-set, 2023] on a cylindrical approximating of the half-space with the height and radius of the cylinder set to be much greater than the size of the conducting spot.

To achieve accurate results using the Boundary Element Method (BEM), the problem is addressed using two distinct meshes, each characterized by a relatively fine granularity and differing in the number of nodes. Results are then extrapolated employing the Richardson extrapolation technique [Richardson, 1911], as illustrated in Fig. 2. Two meshes are built, approximating the same contact area, but with different reference mesh size $h_1 > h_2$ (here, we take $h_1 = 2h_2$). Since the flux converges linearly with the mesh size, a simple extrapolation would give [Zienkiewicz and Taylor, 2000]:

$$Q^* \approx 2Q(h_2) - Q(2h_2).$$

In a general case of two different reference mesh-sizes, we get

$$Q^* \approx \frac{h_1 Q(h_2) - h_2 Q(h_1)}{h_1 - h_2}. \qquad (4)$$

More details and numerical experiments could be found in [Beguin, 2024]. Since we use linear elements, the accuracy of the geometry approximation is also dependent on the mesh size, therefore the Richardson extrapolation is not exact. Nevertheless, for the geometrical error, the ratio of the difference between approximate and the true area to the true area scales as $\approx h^2/(8R^2)$ where $R$ is the curvature radius, therefore we ignore this error.

For the FEM, in addition to mesh-extrapolation we need to introduce a correction related to the finite size geometry, namely, the dependence on the ratio of the average spot radius to the size of the simulated domain, which can be seen as a geometrically related small parameter $\epsilon_g = b/R$ (see Fig. 4). The primary, mesh-related, small parameter is given by $\epsilon_h = h/b$. Then, if we assume a linear dependence on $\epsilon_g$, we can write first terms of the expansion around the exact solution $Q^*$ as:

$$Q(\epsilon_h, \epsilon_g) = Q^* + c_h \epsilon_h + c_g \epsilon_g + o(\epsilon_h, \epsilon_g).$$

In this case, three simulations are needed to determine the value of $Q^*$ and coefficients $c_h$ and $c_g$ in which two values of $\epsilon_h$ and of $\epsilon_g$ should be used in any combination. More details on the mesh and geometry corrections could be found in [Beguin, 2024].
## 3 Conductivity of a not simply connected spot

The total flux of a single circular spot of radius $a$ between two half-spaces made of the material with the same thermal/electric conductivity $k$ was given in [Carslaw and Jaeger, 1959]:

$$Q_0 = 4kaU_0, \qquad (5)$$

where $U_0$ is the difference between electric potential or temperature between the spot and an infinitely remote boundary.

It assumes stationary conductivity and perfect insulation outside the circular spot, i.e. we ignore convective and radiative heat exchange for the thermal problem and ignore tunnel effect or electrical breakdown for the electric problem. Notably, the total flux is proportional to the radius of the spot. The flux distribution within the spot is axially symmetric and thus could be expressed in polar coordinates as a function of radius $r$:

$$j_n^\circ(r) = j \cdot n = \frac{2U_0 k}{\pi \sqrt{a^2 - r^2}}, \qquad (6)$$

where $n$ is the outer normal, see [Carslaw and Jaeger, 1959]. The flux diverges as $1/\sqrt{\xi}$ near the boundary where $\xi$ is the distance to this boundary. This solution is equivalent to the contact pressure of a circular stamp pressed in an elastic half space [Sneddon, 1995] and to the normal stress distribution of an external circular crack [Tada et al., 2000, p. 377]. An analytical solution for conductivity of an elliptic spot on a half-space was later obtained in [Holm, 1957].

In this section, we study the total flux through an annulus spot of different internal radius. The question is to which extent the flux is altered by the presence of small internal holes in not simply connected spots. The internal radius is $r = b$, the external one $r = a$, their ratio is denoted by $\xi = b/a \in [0, 1)$. The boundary conditions remain the same as for the circular spot: the potential is set constant $U = U_0$ at the annulus ($b \le r \le a$) the zero flux is set elsewhere ($r > a$ and $r < b$), at infinity the potential is set to $U = 0$. The geometry is shown in Fig. 3 for three different values of $\xi$. We are mainly interested in the asymptotic evolution of the total flux $Q(\xi)$ in the limit of $\xi \ll 1$. The flux diverges at the two borders of the annulus, however the total flux $Q_a$ should be continuous and decreasing with respect to the relative hole radius $\xi$.

As $\xi \to 0$, the total flux should tend to $Q_0$, whereas for $\xi \to 1$, the flux should vanish $Q_a \to 0$. We could also conjecture, from a physical nature of the phenomenon, that the flux should be a concave function of $\xi$ with zero derivative at $\xi = 0$, meaning that small holes should not affect considerably the flux. We could also conjecture that the local singularity on the inner edge should be less pronounced than on the outer edge being regularized to some extent by the interaction with the whole inner border especially for small size holes.

Smythe [Smythe, 1951] was the first to solve this problem using a superposition method and provided methods to approximately evaluate the flux along with few tabulated results for the annulus conductivity. Other authors [Cooke, 1963, Collins, 1963, Fabrikant, 1993] expanded Smythe's work by reformulating the problem as a triple integral or as a Fredholm integral equation of second kind. However, those solutions do not provide closed form formulas for the flux; the influence of the effect of the hole cannot be easily deduced. Fabrikant [Fabrikant, 1993] proposed an iterative method for the resolution of the integral system. An alternative solution was obtained by Love [Love, 1976] based on constructing upper and lower limits series provided the first terms as:

$$\frac{Q_{\text{Love}}}{Q_0} = 1 - \frac{4}{3\pi^2} \xi^3 - \frac{8}{15\pi^2} \xi^5 - \frac{16}{27\pi^4} \xi^6 - \frac{92}{315\pi^2} \xi^7 - \frac{416}{675\pi^4} \xi^8 + o(\xi^8) \qquad (7)$$

We used the finite element method (FEM) to solve this problem numerically as an axisymmetric problem. The mesh has to be refined near the annulus edges to capture singularities. To do so, two semicircular insertions of radius $r_e$ are constructed near the singularity points (see Fig. 4). The element size at the annulus's edge is set to $h_{\min}$ such that $h_{\min}/b \ll 1$. Far from the annulus, the mesh size is set
<figure><img src="image_5.png" /><figcaption>Figure 3: Geometry of the annular contact spot for various ratios of the internal to external radii $\xi = \{0.2, 0.5, 0.8\}$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Mesh definition and boundary conditions for modeling of an annulus contact spot, using the parameters $\xi = 0.5$, $h_{min}/r_0 = 0.005$ and $R/a = 5$.</figcaption></figure>

to be coarser $h_{max}$, as $h_{max} \gg h_{min}$. Supplementing the mesh refinement near the annulus’ edge, the mesh size is set to $h_e$, at the edge of the semicircular object setting $h_e = 10h_{min}$.

The geometry of the axially symmetric problem is restricted to a square plane of length and width $R \gg b$. The mesh size on the outer edge must be much smaller compared to $b$, i.e. $h \ll b$ and near the internal edge $h \ll a$. However, far from singularities the mesh could be coarser.

The finite element results for the flux distribution for several values of $\xi$ are shown in Fig. 5 and compared with the analytical solution for the circular spot (6). The density of result points allows the reader to judge on the density of mesh used for this solution (the finite element mesh is provided in Supplementary material [Beguin and Yastrebov, 2024]). As expected, the normal flux diverges at both edges and the singularity at the hole (internal edge) is less pronounced than on the outer edge: it decreases faster than the one on the outer edge. Nakamura [Nakamura, 1993] also studied this problem, by both FEM and BEM to assess the BEM accuracy with respect to the FEM resolution, however, the author ignored the existing analytical solution.

Simulation results for the total flux are shown in Fig. 6 with the flux normalized by the one of a circular spot $Q_\circ$ (see Eq. (5)). To identify numerically an asymptotic solution for small holes, we consider a contribution from the hole with a power-law of $\xi$:

$$Q_{\text{fit}} = Q_{\circ} (1 - \alpha k^{\beta}), \qquad (8)$$
whose parameters $\alpha, \beta$ were identified by least squares method as $\alpha \approx 0.1435$ and $\beta \approx 3.028$ in the interval $\xi < \xi_{lim} = 0.2$, which are very close to analytical results of [Love, 1976, Fabrikant, 1993], $\alpha = 4/(3\pi^2) \approx 0.1351$ and $\beta = 3$ (see Eq. (7)). Nakamura's [Nakamura, 1993] results calculated by the BEM are also displayed for comparison, however, because of the lack of Richardson type extrapolation and convergence study, they underestimate the flux value. The Love's solution slightly overestimates the flux for higher values of $\xi$, but this could be readily improved by including a larger number of terms in his series. In fact, all polynomial coefficients of $\xi$ for the series expansion of the flux function are negative and adding new terms will slightly reduce the flux. Nevertheless, the first terms in Eq. (7) are in very good agreement with the numerical results at least for small values of $\xi$.

Expectedly, we conclude that the total flux is very weakly dependent on the presence of small holes in annulus spots because the corrective term is of order $\xi^3$ with a small factor $\sim 0.1$. Therefore, we could conclude that not simply connected spots, at least for holes located far from the outer boundary, conduct almost as well as simply connected spots with the same outer boundary. In addition to this axisymmetric study, one could conduct a similar study but with a hole placed with some eccentricity with respect to the center. Such a study would provide an even stronger argument on the effect of such holes in conductivity problems, however, this study is not included in the scope of the current paper.

<figure><img src="image_7.png" /><figcaption>Figure 5: Finite element results of the normal flux distribution for annulus spot for different values of $k$, dashed line represents the normal flux of a circular spot.</figcaption></figure>

# 4 Conductivity of flower-shaped spots

To mimic complex shapes formed by contact between random rough surface, we first consider a simple geometrical model which we call *flower-shaped* spot whose boundary is described by the following equation in polar coordinates

$$r(\theta) = r_0 + r_1 \cos(n\theta) = r_0 (1 + \xi \cos(n\theta)), \quad r_1 < r_0 \qquad (9)$$
<figure><img src="image_8.png" /><figcaption>Figure 6: Total flux of an annulus as a function of internal to external radius ratio $\xi$: finite element results (circles), the least squares fit of an offset power-law for the asymptotic solution $Q'_{\text{fit}}$ for $\xi < 0.2$ and an approximate analytical solution by [Love, 1976] $Q_{\text{Love}}/Q_{\circ}$. (a) – normalized total flux $Q'(k) = Q(k)/Q_{\circ}$; (b) – normalized flux difference compared to the circular flux $(Q_{\circ} - Q(k))/Q_{\circ}$ in log-log scale highlighting the power-law flux evolution.</figcaption></figure>

where $r_0$ is the mean radius, $r_1$ is the half-length of petals and $n$ is their number. So the two positive dimensionless parameters describing the shape are $\xi = r_1/r_0 < 1$ and $n$. Different flower-shaped spots are presented in the Fig. 7. Note that the average radius does not change with the number of petals nor with their length and is equal to $\langle r \rangle = r_0$. Circles of radius $r_0$ and $r_0(1 + \xi) = r_0 + r_1$ are also shown in the figure; corresponding to minimal and maximal limits for the resulting flux $Q_{\text{min}}(r_0) = Q_{\circ} = 4kr_0U_0$, $Q_{\text{up}} = Q_{\circ}(r_0 + r_1) = 4kr_0(1 + \xi)U_0$,

$$1 \le \frac{Q(r_0, \xi, n)}{Q_{\text{min}}} \le 1 + \xi$$

<figure><img src="image_9.png" /><figcaption>Figure 7: Examples of flower-shaped spots with $\xi = \{0.1, 0.2\}$ and the number of petals $n = \{4, 7\}$.</figcaption></figure>

The particularity of these flower-shaped spots is that the perimeter increases with the number of
petals whereas the area remains constant. They are given by the following equations

$$A = \pi r_0^2 \left(1 + \frac{\xi^2}{2}\right), \qquad (10a)$$

$$P = r_0 \int_0^{2\pi} \sqrt{1 + (n\xi \sin(\theta))^2} d\theta, \qquad (10b)$$

We can notice that the integrand depends on the dimensionless parameter $n' = n\xi \in \mathbb{R}$ which can take arbitrary real values contrary to the integer number of petals $n \in \mathbb{N}$. To characterize the shape of the flower-shaped spot, we could also use compactness $C$ being the ratio of the square root of the area to the perimeter:

$$C(n') = \frac{\sqrt{A}}{P} = \sqrt{\pi} \frac{\sqrt{1 + \xi^2/2}}{4E(in')}, \qquad (11)$$

where $E(x) = \int_0^{\pi/2} \sqrt{1 - x^2 \sin^2(\theta)} d\theta$ is the complete elliptic integral of the second kind and $i$ is imaginary unit. For small values of $\xi$, we could assume that the compactness depends on $n'$ only, for $n' \to 0$, $E(in') \to \pi/2$ and the compactness tends to the maximal value, i.e. the compactness of a circle $C \to C_0 = 1/(2\sqrt{\pi}) \approx 0.282$.

## 4.1 Conductivity results

This study was conducted using both FEM and BEM. Despite the fact that the BEM is more appropriate for half-space approximation, the comparison of two methods with an extrapolation in terms of the size of simulated domain employed in the FEM is also used to assess the validity of the implemented BEM. The flower shape presents a dihedral symmetry $D_n$ which allows us to use only one mesh section (half-petal) with symmetric boundary conditions imposed on the lateral sides out of $2n$ sections needed to construct the full spot for the FEM. The same size reduction can be performed with the BEM using instead the repeatability of the solution. With an increasing number of petals, the angle of the central element near the axis of symmetry sharpens. Consequently, to prevent any deterioration in the solution quality, the smallest simulated angular sector was configured to be approximately $\pi/6$. This setup necessitates simulating a larger segment than what is required by symmetry considerations only.

Simulation results of the normal flux are presented in Fig. 8 for $n = \{4, 7, 10\}$ petals and for $\xi = 0.1$. For the first two cases, the symmetry is fully exploited whereas for $n = 10$ the sector's angle is set to $\pi/5$ to preserve good mesh quality in the center. Meshes are refined near the outer edges where the flux is singular. For $n = 4$ the total flux is $Q \approx 1.0084Q_{min}$, for $n = 7$ the total flux is $Q \approx 1.0150Q_{min}$, and for $n = 10$ the total flux is $Q \approx 1.0209Q_{min}$. So there is a trend to increase the total flux with the increasing number of petals. Visually we can also observe that the singularity in the trough (petal's root) is weaker than the one near the crest (petal's extremity). The more petals we have, the weaker the flux intensity in the trough because of the increasing interaction with the neighboring petals; a similar trend was observed for the annular spot for small internal radii.

In total, 97 simulations<fnref n="1" /> were carried out for $\xi = 0.1$ and $n \in (1, 256)$ as well as for $\xi = 0.2$ and $n \in (1, 100)$. The resulting total flux with respect to the normalized number of petals $n'$ is presented in Fig. 9. Those results are computed with Richardson extrapolation obtained by two meshes of different density for FEM and BEM and using a similar extrapolation for domain's dimensions in case of FEM simulations.

The total flux, offset by the minimal flux $Q - Q_0$ could be normalized by the difference between the maximal and the minimal theoretical fluxes, corresponding to circular spots of radii $r_0(1 + \xi)$ and $r_0$, respectively, i.e. we get in the denominator $Q_{up} - Q_0 = 4kr_0\xi U_0 = \xi Q_0$. This normalization results in a

<footnote n="1">¹Without counting simulations used to employ Richardson extrapolation.</footnote>
<figure><img src="image_10.png" /><figcaption>Figure 8: Simulation results for the normal flux for flower-shaped spots with $n = 4, 7, 10$ and $\xi = 0.1$.</figcaption></figure>

universal curve for the total normalized flux for any $r_0$ and $\xi$:

$$Q' = \frac{Q - Q_{\circ}}{Q_{\text{up}} - Q_{\circ}} = \frac{Q - Q_{\circ}}{\xi Q_{\circ}} \qquad (12)$$

The normalized flux evolution seems to be logarithmic, but from the physical point of view, the flux cannot overpass $Q_{\text{up}}$, therefore we shall require that the bounds $0 \le Q' \le 1$ are respected for all $n'$. The resulting flux is well fitted by a two-parameter function which, however, was found empirically:

$$Q'_{\text{fit}}(n') = a \left( 1 - \frac{1}{bn' + 1} \right), \quad 0 < a \le 1. \qquad (13)$$

This fit function is plotted in Fig. 9 along with FEM and BEM simulation results. The coefficients determined by least squares fit are presented in Table 1. The slope at $n' = 0+$ is equal to the product ab. Even if the coefficients a and b are slightly different for different sets, this slope remains close for all three independent fits, and roughly equals to 0.3. Combining Eqs. (12) and (13), we obtain the following phenomenological equation for the total flow of a flower-shaped spot:

$$Q = Q_{\circ} \left( 1 + a\xi \left( 1 - \frac{1}{bn\xi + 1} \right) \right) \approx Q_{\circ} \left( 1 + 0.923\xi \left( 1 - \frac{1}{0.326n\xi + 1} \right) \right), \qquad (14)$$

where $Q_{\circ} = 4kr_0U_0$. For the infinite number of petals n of finite half-length factor $\xi$, there is a limit flux given by this fit, $\lim_{n' \to \infty} (Q'_{\text{fit}}) = a \approx 0.923$ and this limit is independent of $\xi$. However, the validity of the suggested fit beyond the studied interval of $n'$ cannot be taken for granted. An argument in favor of such a limit $a < 1$, i.e. that the total flux for the infinite number of petals remains below the flux of a circular shape of radius $r_0(1 + \xi)$, could me made based on the area of the conductive spot. Indeed, the area of the full circular spot is considerably bigger than this of the flower-shaped spot $\pi r_0^2(1 + \xi)^2 > \pi r_0^2(1 + \xi^2/2)$. On the other hand, the small features (infinitely thin petals) should not strongly affect the conductivity of the spot, thus suggesting that possibly the flux should simply tend to 1... However, the current fit function could not be properly approximate the data if one sets $a = 1$. The question of a rigorous definition of the limit flux for the infinite number of petals *remains open*. But our initial guess, for the fitted parameters is given by the following limit flux:

$$\lim_{n\xi \to \infty} Q \approx Q_{\circ}(1 + 0.923\xi) < Q_{\text{up}} = Q_{\circ}(1 + \xi). \qquad (15)$$

Even though such a flower-shaped geometry, especially in the limit of infinite number of petals, is not very relevant to contact problems between isotropic surfaces, which was at the origin for this study, this limit value presents an interesting by-product of this study. Among other results one can deduce a
relation between $n$ and $\xi$ which ensures $x$ fraction conductivity in the interval $Q_\circ$ and $Q_{up} = (1 + \xi)Q_\circ$, with $Q = Q_\circ$ for $x = 0$ and $Q = Q_\circ(1 + \xi/2)$ for $x = 0.5$:

$$x = a \left( 1 - \frac{1}{bn\xi + 1} \right) \Leftrightarrow n' = n\xi = \left[ \frac{x}{b(a-x)} \right] \Rightarrow n\xi \approx \left[ \frac{x}{0.326(0.923-x)} \right], \quad (16)$$

Therefore, to reach the mean flux between two limits, i.e. for $x = 0.5$, one would need a spot with $n\xi \approx 3.61$, i.e. for $\xi = 0.2$ one would need approximately 18 petals and for $\xi = 0.1$ a double of that. However, to reach 75% ($x = 0.75$), for $\xi = 0.1$, one would need a spot with approximately 130 petals.

<table><caption>Table 1: Least squares fit for coefficients of Eq. (13) for the sets of results of flower-shaped spot obtained using FEM and BEM simulations.</caption><thead><tr><th rowspan="2">Simulation</th><th colspan="2">Parameters</th><th colspan="3">Coefficients</th></tr><tr><th>ξ</th><th>n ∈</th><th>a</th><th>b</th><th>ab</th></tr></thead><tbody><tr><td>FEM</td><td>0.1</td><td>[1, 150]</td><td>0.928</td><td>0.327</td><td>0.304</td></tr><tr><td>FEM</td><td>0.2</td><td>[1, 100]</td><td>0.923</td><td>0.326</td><td>0.301</td></tr><tr><td>BEM</td><td>0.1</td><td>[1, 256]</td><td>0.923</td><td>0.326</td><td>0.301</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>Figure 9: Simulation results for the normalized total conductivity for the flower-shaped spot as a function of normalized petal's parameter $n'$ for different $\xi$, three independent least squares fit of function (13) are also plotted; the three corresponding tangents at the origin are also plotted.</figcaption></figure>
## 4.2 Alternative "multi-petal" spots

The same conductivity study could be conducted on other simple forms possessing a single-scale petal-like structures with the same symmetry properties. Specifically we identified the following shapes: star-shaped and gear-shaped spots shown in Fig. 10 and Fig. 11, respectively. For "stars", each petal is made up by straight lines connecting the roots and extremities of "petals", i.e. points with radial coordinates $r_0(1 - \xi)$, $r_0(1 + \xi)$. The number of "petals" (or "rays") as previously is denoted by $n$, and the half-petal length $r_1$ is again determined by the ratio $r_1 = \xi r_0$. The gear-shaped spots are made of circular arcs with constant $r = r_0(1 - \xi)$ and $r = r_0(1 + \xi)$ over equal angular segments. Contrary to $C^\infty$ flower-shaped spots, star-shaped ones are only of class $C^0$ with respect to $\theta$ whereas gear-shaped spots represent multivalued mapping, so they are not even injective even for a single "petal" or "tooth".

<figure><img src="image_12.png" /><figcaption>Figure 10: Examples of star-shaped spots with $n \in \{10, 20\}$ petals, and half-petal length defined by $\xi \in \{0.1, 0.2\}$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11: Examples of gear-shaped spots with $n \in \{10, 20\}$ petals, and half-petal length defined by $\xi \in \{0.1, 0.2\}$.</figcaption></figure>

The area of these shapes can be readily expressed by an elementary sum of triangles and a regular polygon for the stars, or circular sectors for the gears. Similarly, the perimeters are easy to find, and the resulting compactness can be also readily computed:

$$ 
\begin{aligned}
P_{\text{star}} &= 2\sqrt{2}nr_0\sqrt{1 + \xi^2 - (1 - \xi^2)\cos(\pi/n)}, \\
A_{\text{star}} &= nr_0^2(1 - \xi^2)\sin(\pi/n), \\
C_{\text{star}} &= \frac{\sqrt{(1 - \xi^2)\sin(\pi/n)}}{2\sqrt{2n}\sqrt{1 + \xi^2 - (1 - \xi^2)\cos(\pi/n)}}.
\end{aligned} 
\quad (17) 
$$
$$ 
\begin{aligned}
P_{\text{gear}} &= 2\pi r_0 (1 + 2n'/\pi), \\
A_{\text{gear}} &= \pi r_0^2 (1 + \xi^2), \\
C_{\text{gear}} &= \frac{\sqrt{(1 + \xi^2)}}{2\sqrt{\pi}(1 + 2n'/\pi)}.
\end{aligned} 
\qquad (18) $$

All those geometric features are summarized in appendix A and in Fig. 25. In summary, the star-shaped area converges to

$$ A_{\text{star}} \xrightarrow{n \to \infty} \pi r_0^2 (1 - \xi^2), $$

but for all n it is always smaller than the area of gear- and flower-shaped spots:

$$ A_{\text{star}} < A_{\text{flower}} < A_{\text{gear}}. $$

The gear-shaped spots also have a bigger perimeter than the one of flower and star for a given number of petals:

$$ P_{\text{flower}} \approx P_{\text{star}} < P_{\text{gear}}. $$

Finally, the gear-shaped spots appear to be the least compact shape, while the flower shape is as compact as the star

$$ C_{\text{gear}} < C_{\text{flower}} \approx C_{\text{star}}. $$

The conductive simulations for these shapes were conducted using the fast-BEM. The star- and gear-shaped spots present the same dihedral symmetry $D_n$ as for the flower-shaped, which allows us again to reduce the problem size to some extent. Example of simulation results are presented in Fig. 12. In total, 48 simulations were performed for both gear- and star-shaped spots, for $\xi = 0.1$, and $n \in (4, 256)$. As previously, the conductive property is assessed by computing the overall flux and helped by the Richardson extrapolation. Those are again normalized according to (12), and finally presented in Fig. 13 complemented with previous results for the flower-shaped spots. The least square fitted parameters for the same Eq. (13) are shown in Table 2.

Qualitatively, all three types of shapes show the same trend in the total flux evolution with the number of "petals": an initial steep increase and further saturation to a more or less constant value. The thermal conductivity of the gears is higher than that of the flowers, while the star-shaped configuration displays a lower conductivity. We could attribute this ordering to the only basic geometric parameter which significantly differs for all three of them, namely the area. Another consideration is the amount of area located closer to the outward boundary, which of course is higher for the gear like geometry than those of flower and star. In the limit of the infinite number of petals the following result (initial guess obtained by extrapolation) is obtained (see parameter a in Table 2):

$$ Q_{\text{star}}^{\lim} \approx 0.903 < Q_{\text{flower}}^{\lim} = 0.923 < Q_{\text{gear}}^{\lim} = 0.978. $$

However, we would like to highlight once again that these values must be seen as a first guess, and a more rigorous assessment (e.g., using an accurate asymptotic analysis) is needed.

<table><caption>Table 2: Fit parameters for Eq. (13) for the total flux of different multi-petal shapes.</caption><thead><tr><th rowspan="2">Simulation</th><th colspan="2">Parameters</th><th colspan="3">Coefficients</th></tr><tr><th>ξ</th><th>n ∈</th><th>a</th><th>b</th><th>ab</th></tr></thead><tbody><tr><td>Gear</td><td>0.1</td><td>[4, 256]</td><td>0.978</td><td>0.894</td><td>0.875</td></tr><tr><td>Flower</td><td>0.1</td><td>[1, 256]</td><td>0.923</td><td>0.326</td><td>0.301</td></tr><tr><td>Star</td><td>0.1</td><td>[4, 256]</td><td>0.903</td><td>0.220</td><td>0.199</td></tr></tbody></table>
<figure><img src="image_14.png" /><figcaption>Figure 12: BEM result of the flux through a gear- and star-shaped spots of 20 petals and $\xi = 0.1$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 13: Results of the normalized flux for star-, gear- and flower-shaped spots with respect to the normalized "petal" parameter $n' = \xi n$. Least squares fit of Eq. (13) is also shown.</figcaption></figure>
# 5 Conductivity of self-affine spots

## 5.1 Geometry of spots

Being inspired by shapes of contact clusters occurring in contact of random rough surfaces (see Fig. 1), in this section we study contact spots of model complex shapes presenting some randomness. The shapes under study (see Fig. 15) could also recall coffee or ink stains. To take up the Archard’s image of ”protuberances on protuberances on protuberances” [Archard, 1957], we constructed contact-spots in a self-affine fashion by summing-up multiple harmonics. The goal is to imitate to some extent realistic contact spots occurring for surfaces with a rich spectral content and to expand the results obtained for single-harmonic flower-shaped spots to more complex forms.

The first step to generate a spot with self-affine boundary is to introduce a periodic perturbation function $h(\theta)$ as a superposition of cosines which individually could be seen as flower-shaped spots:

$$h(\theta) = \sum_{k=k_l}^{k_s} \xi_k \cos(k\theta + \theta_k^0), \quad h(\theta) = h(\theta + 2\pi), \quad \langle h \rangle = 0 \quad (19a)$$

$$\xi_k = \xi \left( \frac{k}{k_l} \right)^{-(0.5+H)} \quad (19b)$$

The summed up harmonics include all integer modes from a fixed interval $k \in (k_l, k_s)$ with amplitudes $\xi_k$ which decay as a strict power-law of the mode number with an exponent involving the Hurst exponent $H \in (0, 1)$ ensuring a self-affinity of the boundary. The randomness is provided exclusively by the phases $\theta_k^0$ which follow a uniform distribution on $\theta_k^0 \in [-\pi, \pi)$. The perturbation $h(\theta)$ thus constructed follows a Gaussian distribution. The power spectral density (PSD) decays as a power law of the wavenumber with the exponent $-(1+2H)$. This richness of the spectrum could be defined by the ”magnification” parameter $\zeta$ presenting the ratio of the highest to lowest wavenumbers $\zeta = k_s/k_l$ [Persson, 2001].

The radius of a contact spot in polar coordinates $r(\theta)$ can be readily defined with the perturbation $h(\theta)$ as:

$$r(\theta) = r_0(1 + h(\theta)), \quad (20)$$

naturally $\langle r \rangle = r_0$. Nevertheless, in this construction, even imposing $\xi < 1$ does not guarantee positivity, as the factor $1+h(\theta)$ may become negative, consequently leading to negative values for $r(\theta)$ as well. To overcome this problem, some smarter transformation should be used, for example:

$$r(\theta) = r_0 \exp(h(\theta)), \quad (21)$$

whose two first terms of Taylor expansion is equivalent to (20), but the transformation (21) keeps the final shape well defined without self-intersections: even for $h \to -\infty$, $r \to 0$. However, this transformation does not preserve the mean radius at $r_0$; the radius will change with $H$, $\xi$, $k_l$ and $\zeta$; this deviation will be characterized by a dimensionless quantity $\bar{r} = \langle r(\theta) \rangle / r_0$. Within this formulation the parameter $\xi$ plays a similar role as in the study of flower-shaped spots: here, to the first order it presents the ratio between the amplitude of the first mode to the nominal radius. The two transformations (20),(21) are illustrated in Fig. 14 highlighting a situation, where for transformation (20) the radius becomes negative. Therefore, for this study we adopt the exponential transformation (21). Several examples of complex shapes generated using the presented algorithm are shown in Fig. 15 for $\xi = 0.1$ and different values of $k_l, k_s$ and $H$ but with the same set of random phases $\theta_k^0$. As the Hurst exponent increases, the spot naturally becomes smoother.

## 5.2 Geometrical characteristics

The initial height perturbation $h$ follows normal distribution with zero mean and standard deviation $\sigma_h$, but it is not preserved by the exponential transformation (21). The obtained radius follows a log-normal
<figure><img src="image_16.png" /><figcaption>Figure 14: Example of self-affine spot construction for $\xi = 0.25$, $k_l = 4$, $k_s = 32$ and $H = 0.5$: (a) height perturbation $h(\theta)$ and (b) resulting spot radius for linear (20) and exponential (21) transformations.</figcaption></figure>

distribution

$$P(r) = \frac{1}{r\sigma_h\sqrt{2\pi}} \exp\left(-\frac{\ln^2(r)}{2\sigma_h^2}\right). \qquad (22)$$

Histograms shown in Fig. 16 present the probability density of the spot radius constructed for $H = 0.25$, $k_l = 8$, $k_s = 16$ and two values of $\xi = \{0.05, 0.1\}$ and computed over 1000 generated spots, the least squares fitted normal and log-normal densities are also plotted. For small values of $\xi$, the distribution is very close to the normal one, whereas for a higher value, it clearly follows the log-normal one.

The standard deviation $\sigma_h$ or the variance $\sigma_h^2$ of the height perturbation $h(\theta)$ can be computed taking into account the orthogonality of cosine functions:

$$\sigma_h^2 = \frac{\xi^2}{2} \sum_{k=l}^{k_s} \left(\frac{k}{k_l}\right)^{-(1+2H)} \qquad (23)$$

After the transformation, the mean radius $\langle r \rangle$ and the variance of the height distribution $\sigma_r^2$ can be found as

$$\langle r \rangle = r_0 \exp(\sigma_h^2/2), \qquad (24a)$$

$$\sigma_r^2 = [\exp(\sigma_h^2) - 1] \exp(\sigma_h^2), \qquad (24b)$$

The mean value of the radius is no longer equal to $r_0$, but it tends to $\langle r \rangle \to r_0$ as $\sigma_h \to 0$. Note that the variance could be expressed through the first terms of the Taylor expansion as $\sigma_r^2 \approx \sigma_h^2 + 1.5\sigma_h^4 + 7/6\sigma_h^6 + O(\sigma_h^8)$, which demonstrates that for small values of $\sigma_h$, $\sigma_r \approx \sigma_h$ with a high accuracy. The comparison of the analytical expression of the variance Eq. (24b) with numerically evaluated standard deviation is presented in Appendix B.

The standard deviation of radius is an important geometric characteristic of the spot. However, the standard deviation of the gradient and Laplacian of the radius as well as Nayak parameter [Nayak, 1971] could also have an effect on the conductivity of the spot. These geometrical characteristics are related to
<figure><img src="image_17.png" /><figcaption>Figure 15: Examples of self-affine spots and their geometrical characteristics for $\xi = 0.1$.</figcaption></figure>

In the upper panel (a,b,c) $k_l = 2$, $k_s = \{8, 16, 32\}$, $H = 0.25$: (a) $\bar{r} = 1.0017$, $\sigma = 0.057$, (b) $\bar{r} = 1.0020$, $\sigma = 0.063$, (c) $\bar{r} = 1.0023$, $\sigma = 0.067$;

in the middle panel (d,e,f) $k_l = 4$, $k_s = \{8, 16, 32\}$, $H = 0.25$: (d) $\bar{r} = 1.0019$, $\sigma = 0.062$, (e) $\bar{r} = 1.0029$, $\sigma = 0.075$, (f) $\bar{r} = 1.0036$, $\sigma = 0.084$;

in the lower panel (g,h,i) $k_l = 4$, $k_s = 128$, $H = \{0.25, 0.5, 0.75\}$: (g) $\bar{r} = 1.0044$, $\sigma = 0.094$, (h) $\bar{r} = 1.0028$, $\sigma = 0.074$, (i) $\bar{r} = 1.0020$, $\sigma = 0.063$.

spectral moments $m_0$, $m_2$ and $m_4$ as follows:

$$m_0 = \sigma^2, \qquad (25a)$$

$$m_2 = \langle |\nabla r|^2 \rangle = \frac{1}{2\pi} \int_0^{2\pi} \left( \frac{1}{r(\theta)} \frac{\partial r}{\partial \theta} \right)^2 d\theta, \qquad (25b)$$

$$m_4 = \frac{1}{2\pi} \langle |\Delta r|^2 \rangle = \int_0^{2\pi} \left( \frac{1}{r^2(\theta)} \frac{\partial^2 r}{\partial \theta^2} \right)^2 d\theta, \qquad (25c)$$
<figure><img src="image_18.png" /><figcaption>Figure 16: Distribution of the radius for constructed self-affine spots for $H = 0.25$, $k_l = 8$, $k_s = 16$ and (a) $\xi = 0.05$ and (b) $\xi = 0.1$. Least squares fit of the normal and log-normal distributions is also presented.</figcaption></figure>

In the limit of infinitesimal perturbations $\xi \ll 1$, we can use the following analytical equations for the 2nd and 4th spectral moments (the 0th moment is nothing but the variance of radius computed in (24b)):

$$m_2 = \frac{\xi^2}{2} \sum_{k_l} k^2 \left(\frac{k}{k_l}\right)^{-(1+2H)}, \quad m_4 = \frac{\xi^4}{2r_0^2} \sum_{k_l} k^4 \left(\frac{k}{k_l}\right)^{-(1+2H)} \quad (26)$$

The comparison of Eq. (26) with numerically evaluated moments is presented in Appendix B. For flower-shaped spots, the spectral moments simplify to the following forms:

$$m_0^f = \frac{r_0^2}{2} = \frac{\xi^2 r_0^2}{2}, \quad m_2^f = \frac{r_1^2 n^2}{2r_0^2} = \frac{\xi^2 n^2}{2} = \frac{n'^2}{2}, \quad m_4^f = \frac{r_1^2 n^4}{2r_0^4} = \frac{\xi^2 n^4}{2r_0^2} = \frac{n'^2 n^2}{2r_0^2} \quad (27)$$

Since the normalized conductivity for flower-shaped spots was shown to depend exclusively on $n' = \xi n$ (see Eq. (13)), we could suggest that for a similar normalization, the main characteristic affecting the conductivity of self-affine spots will be the standard deviation of the gradient $\sqrt{\langle |\nabla r|^2 \rangle} = \sqrt{m_2}$. In addition, it could be shown that the area of a spot is an affine function of $m_0$: $A \approx \pi r_0^2 (1 + am_0)$, and its perimeter is an affine function of $\sqrt{m_2}: P \approx 2\pi r(1 + b\sqrt{m_2})$, where a, b are positive constants.

## 5.3 Normalization and result

In contrast to the study of flower-shaped spots, this study is no longer deterministic and requires taking into account the randomness of the studied geometries. The total flux for a given geometry $Q_i$, defined as an independent event, is assumed to have a consistent average value $\mu(Q)$ and a standard deviation $\sigma(Q)$ across the same set of spot-generative parameters. This study aims to understand the average behavior based on these parameters. To achieve this, we compute average values for numerous realization of spot geometries. However, this yields only an approximate value for $\langle Q \rangle$, which depends on the number of realizations $n$. In practice, the standard deviation of the mean value scales as $\sigma(\langle Q \rangle) = \sigma(Q)/\sqrt{n}$.
The Bienaymé-Chebyshev inequality aids in establishing a confidence interval implying a parameter $\gamma \in (0, 1)$: the probability to find the mean value $\langle Q \rangle$ outside the interval $\pm\sigma(Q)/\sqrt{n\gamma}$ around its theoretical value $\mu(Q)$ is smaller than $\gamma$ whatever the true underlying distribution, i.e.

$$P\left(\left|\mu(Q) - \langle Q \rangle\right| \ge \frac{\sigma(Q)}{\sqrt{n\gamma}}\right) \le \gamma, \qquad (28)$$

Equivalently, the probability is $(1 - \gamma)$ to find the mean value in the confidence interval $\pm\sigma(Q)/\sqrt{n\gamma}$. When there are 11 simulations, the interval of confidence spans approximately $\pm 3\sigma(Q)$ with the probability of 99%, effectively encapsulating the mean value $\langle Q \rangle$. In order to reduce the interval to one standard deviation with the same accuracy, i.e. $\sqrt{n\gamma} = 1$ for $\gamma = 0.01$ one would need to carry out $n = 100$ simulations which is computationally expensive in view of the number of parameters to be studied. For $n = 11$ the probability to find the mean value in one standard deviation interval, i.e. $\sqrt{11\gamma} = 1$ is higher than $\approx 9\%$. While the Bienaymé-Chebyshev inequality provides a rigorous lower bound, the actual accuracy of our results can overpass this conservative limit. To balance the computational efforts and the accuracy, the number of BEM simulation results per combination of parameters was set to $n = 11$, and this dataset was utilized to estimate the mean value and the confidence interval obtained from the measurement of the standard deviation. As seen in practice for this problem the standard deviation is much smaller than the approximate mean value, therefore the confidence interval is very narrow and $n = 11$ seems to be a good compromise.

An example of flux distribution obtained by fast-BEM is displayed in Fig. 17 for $\xi = 0.05$, $k_l = 8$, $k_s = 128$ and $H = 0.25$. The normal flux remains singular at the edge but less pronounced at troughs than at crests as shown previously for the flower-shaped contact spot. Compared to the latter, it appears more difficult to construct a good mesh for self-affine spots efficiently (fine mesh near the border and coarse far from border). The mesh size was prescribed as a function of the edge curvature and as a function of the shortest distance to the border. See zoom in Fig. 17, the finest used mesh reach $N_e = 34340$ elements. As previously, to employ the Richardson extrapolation, two meshes of different density were used to obtain accurate results.

<figure><img src="image_19.png" /><figcaption>Figure 17: Example of simulation results representing the flux distribution for a self-affine spot with $\xi = 0.05$, $k_l = 8$, $k_s = 128$ and $H = 0.25$.</figcaption></figure>

For the global flux analysis, the following geometrical characteristics describe sufficiently well the geometry $\mathcal{A} = \{\langle r \rangle, \sqrt{m_0}, \sqrt{m_2}, \sqrt{m_4}, H\}$: (1) the mean radius $\langle r \rangle$, the standard deviation of (2) radius
$\sqrt{m_0}$, (3) of its gradient $\sqrt{m_2}$, (4) of its Laplacian $\sqrt{m_4}$ and (5) the Hurst exponent. The initial set of independent generative parameters is $I = \{r_0, \xi, k_l, k_s, H\}$ and contains (1') radius $r_0$, (2') amplitude of perturbation $\xi$, (3') lower $k_l$ and (4') upper $k_s$ cutoffs, and (5) the Hurst exponent $H$. The set $\mathcal{A}$ is easy to measure for any spot, and it allows to study the effect of individual parameters on the total flux; the generative set $I$ is complete and dimensionless, but suffers from the fact that it could not be easily derived for an arbitrary geometry. A sensitivity analysis will be carried out on both sets.

The first task is to normalize the total flux $Q$ produced by self-affine spots. There are two main options: it can be normalized by the flux of an equivalent circular contact spot (a) of the same mean radius $Q_0 = 4kU_0\langle r \rangle$ or (b) of the same area. The first option was selected since the mean radius also enters the set of parameters, moreover $Q_0$ is equivalent to the definition of $Q_{min}$ used for the flower-shaped contact spot. The normalized total flux is thus defined by

$$Q' = \frac{Q}{Q_0} = \frac{Q}{4kU_0\langle r \rangle} \qquad (29)$$

Therefore, since the problem does not have an internal length, we can exclude the mean radius $\langle r \rangle$ from the set of parameters defining the total flux $Q'$ and consider two sets of dimensionless parameters:

* Geometrical set of parameters $\mathcal{A}' := \{\sqrt{m_0}/\langle r \rangle, m_2, \langle r \rangle \sqrt{m_4}, H\}$

* Generative set of parameters $I' = \{\xi, k_l, k_s, H\}$

By analogy with the maximal total flux $Q_{up}$ defined for flower-shaped spots, we could define an equivalent upper limit for self-affine spots. It is not a good idea to define it as the maximal radius of the self-affine spot as it could tend to infinity for very high number of modes. However, the variance of the flower-shaped spot $m_0^f$, Eq. (27), provides us with a hint getting back the half-petal length $r_1$ in another way, i.e. using this equation we can express it as $r_1 = (2m_0^f)^{1/2}$; therefore, the limit characteristic radius could be expressed as $\langle r \rangle + \sqrt{2m_0}$. Thus, the difference between two flows used for normalization takes the following form:

$$Q_{up} - Q_0 = 4kU_0\sqrt{2m_0}$$

Following the same renormalization between zero and one as was used for the flower-shaped spots, we could define the renormalized flux as:

$$Q'' = \frac{Q - Q_0}{Q_{up} - Q_0} = \frac{Q}{4kU_0\sqrt{2m_0}} - \frac{\langle r \rangle}{\sqrt{2m_0}} \qquad (30)$$

The total flux normalized according to Eq. (29) is presented in Fig.18 for all simulated self-affine spots. These spots are constructed by changing the lower cutoff $k_l = \{2, 4, 8\}$ and for four values of the upper cutoff $k_s = \zeta k_l$ with the magnification $\zeta = \{4, 8, 16, 32\}$. The Hurst parameter $H$ takes the values $H = \{0.25, 0.3, 0.4, 0.5, 0.6, 0.7, 0.75\}$. The colors are used to distinguish the 3 sets of the results according to different $k_l$. In each color set, the results are distinguished by their marker style according to values of magnification $\zeta$, moreover, the higher the $\zeta$, the darker the color. Along every result-curve the Hurst exponent $H$ changes as shown by the arrow: the smaller the $H$, the higher the flux. The curves are entwined together, but they seem to follow the same trend. Plotting the data with respect to geometrical parameters $\mathcal{A}'$ offers a better representation than the use of the generative set of parameters $I'$. The variation in slope seems to be controlled by parameter $\zeta = k_s/k_l$: increasing $\zeta$ increases the average slope with respect to $\sqrt{m_2}/\langle r \rangle$ and decreases the average slope with respect to $\sqrt{m_2}$ and $\langle r \rangle \sqrt{m_4}$.

The results of the renormalized flux Eq.(30) with respect to $\sqrt{m_2}$ and $\langle r \rangle \sqrt{m_4}$ are presented in Fig. 19. However, the role of the fourth moment in the form $\langle r \rangle \sqrt{m_4}$ seems to be strongly correlated with $\sqrt{m_2}$ and does not bring much additional information. The lack of simple dependency of the normalized flux with respect to geometrical characteristics pushes us to suggest an alternative normalization. The shown results seem to depend strongly on $k_l$ parameter. By exploring a wider spreading of $\sqrt{m_2}$ for different $k_l$, we hope to easier identify their influence, which is the objective of the following sections.
<figure><img src="image_20.png" /><figcaption>Figure 18: Normalized flux (29) of self-affine spots: for the following generative parameters $I'$: $\xi = 0.05$, $k_l = \{2,4,8\}$, $\zeta = \{4,8,16,32\}$ and $H = \{0.25,0.3,0.4,0.5,0.6,0.7,0.75\}$, note that $k_s = \zeta k_l$. The mean normalized flux is plotted with a marker and shaded around according to the confidence interval half-width of $\pm 1.35\sigma$ defined for $\gamma = 0.05$ (namely with a rate of confidence of 95%).</figcaption></figure>

To provide a visual geometrical interpretation in the flux variation (see Fig. 20) related to geometrical characteristics of self-affine spots, we present particular shapes along with the values of corresponding geometrical characteristics and of the total flux in Table 3, the location of these particular spots is also highlighted in Fig. 19. The pairs $\{S_1, S_2\}$, $\{S_3, S_4\}$ and $\{S_5, S_6\}$ have close values of $\sqrt{m_2}$. The difference between the flux of spots $S_5$ and $S_6$ highlights the fact that the flux does not depend only on the second moment. Nevertheless, this could be seen as a second order effect compared to that of $\sqrt{m_2}$. Spots $S_7$ and $S_8$ are among the "roughest" spots and possess the highest flux.

<table><caption>Table 3: Parameters for spots shown in Fig. 20 and the resulting total flux.</caption><thead><tr><th rowspan="2">Spot #</th><th rowspan="2">Q''</th><th colspan="3">Parameters <i>A</i>'</th><th colspan="4">Parameters <i>I</i>'</th></tr><tr><th>√<i>m</i><sub>0</sub>/⟨<i>r</i>⟩</th><th>√<i>m</i><sub>2</sub></th><th>√<i>m</i><sub>4</sub>⟨<i>r</i>⟩</th><th>ξ</th><th><i>k</i><sub>l</sub></th><th><i>k</i><sub>s</sub></th><th><i>H</i></th></tr></thead><tbody><tr><td>S<sub>1</sub></td><td>0.151</td><td>0.0664</td><td>0.658</td><td>14</td><td>0.05</td><td>2</td><td>32</td><td>0.25</td></tr><tr><td>S<sub>2</sub></td><td>0.172</td><td>0.0664</td><td>0.666</td><td>12.7</td><td>0.05</td><td>4</td><td>32</td><td>0.6</td></tr><tr><td>S<sub>3</sub></td><td>0.173</td><td>0.0657</td><td>0.958</td><td>39.8</td><td>0.05</td><td>2</td><td>64</td><td>0.3</td></tr><tr><td>S<sub>4</sub></td><td>0.184</td><td>0.0645</td><td>0.972</td><td>63.8</td><td>0.05</td><td>4</td><td>128</td><td>0.7</td></tr><tr><td>S<sub>5</sub></td><td>0.317</td><td>0.0800</td><td>2.06</td><td>163</td><td>0.05</td><td>4</td><td>128</td><td>0.4</td></tr><tr><td>S<sub>6</sub></td><td>0.380</td><td>0.0852</td><td>1.97</td><td>130</td><td>0.05</td><td>8</td><td>128</td><td>0.75</td></tr><tr><td>S<sub>7</sub></td><td>0.644</td><td>0.126</td><td>5.12</td><td>446</td><td>0.05</td><td>8</td><td>128</td><td>0.25</td></tr><tr><td>S<sub>8</sub></td><td>0.749</td><td>0.182</td><td>8.77</td><td>1493</td><td>0.05</td><td>8</td><td>256</td><td>0.25</td></tr></tbody></table>
<figure><img src="image_21.png" /><figcaption>Figure 19: Renormalized flux, Eq. (30) of self-affine spots: for $\xi = 0.05$, $k_l = \{2, 4, 8\}$, $\zeta = \{4, 8, 16, 32\}$, $H = \{0.25, 0.3, 0.4, 0.5, 0.6, 0.7, 0.75\}$. The mean normalized flux is plotted with a marker and shaded around according to the confidence interval half-width of $\pm 1.35\sigma$ and of 95% of accuracy; vertical lines and numbers correspond to spots shown in Fig. 20.</figcaption></figure>

## 5.4 Results with renormalized standard deviation

As presented in Figs. 18 and 19, the results are clustered with respect to $k_l$. To have more control on geometrical characteristics, we renormalize the generative function $h(\theta)$ in order to prescribe its dimensionless standard deviation $\sigma_h = \sqrt{m_{0,h}}$ (see Eq. (23)):

$$h(\theta) = \sqrt{\frac{2m_{0,h}}{s}} \sum_{k=k_l}^{k_s} \xi_k \cos(k\theta + \theta_k^0), \quad s = \sum_{k=k_l}^{k_s} \xi_k^2 \qquad (31)$$

with $\xi_k$ defined by Eq. (19b). The exponential transformation from $h(\theta)$ to radius $r(\theta)$ (21) remains intact. Then, for the normalized generative parameters we have: $I'' = \{\sqrt{m_{0,h}}, k_l, k_s, H\}$ and for the geometrical ones we still have $\mathcal{A}'$. The main goal for such a choice is to decorrelate $\sqrt{m_0}/\langle r \rangle$ and $\sqrt{m_2}$ and thus to level down $m_0$ for $k_l = 8$, and to level it up for $k_l = 2$. Note also that $\sqrt{m_0}/\langle r \rangle \approx \sqrt{m_{0,h}}$ for small values of the latter, see (24b) and its Taylor expansion.

The results of this set of simulations for the renormalized flux (30) are presented in Fig. 21 with respect to the standard deviation of the radius gradient $\sqrt{m_2}$; lower cutoffs $k_l = 2$ (orange circles and red squares) and $k_l = 8$ (green crosses and cyan triangles) were used. In addition, different colors correspond to different magnifications: $\zeta = 4$ for red squares and cyan triangles, $\zeta = 8$ for orange circles and green crosses. The Hurst exponent takes three values $H = \{0.25, 0.50, 0.75\}$. For the same value of $\sqrt{m_{0,h}} \approx \sqrt{m_0}/\langle r \rangle$, thanks to the variation of the Hurst exponent, the value of $\sqrt{m_2}$ varies within a certain interval, such data points are connected by a line. Every set of such lines (of the same shade) align along their master curve. Such results demonstrate that even though the standard deviation $\sqrt{m_0}/\langle r \rangle$ controls the thermal flux to a large extent, the standard deviation of the gradient $\sqrt{m_2}$ also influences the result. This conclusion is possible since $\sqrt{m_2}$ does not enter the flux normalization (30). The flux increases with respect to both $\sqrt{m_0}/\langle r \rangle$ and $\sqrt{m_2}$ as well as with respect to $\langle r \rangle \sqrt{m_4}$ because of the strong correlation of the latter with the second moment. For an equivalent $\sqrt{m_2}$, the flux is higher for spots with a lower
<figure><img src="image_22.png" /><figcaption>Figure 20: Examples of self-affine spots with corresponding generative characteristics; the corresponding data points are highlighted in Fig. 19, the resulting flux and geometrical characteristics can be found in Table 3.</figcaption></figure>

magnification $\zeta$.

To provide a geometrical meaning to these results, coupled pairs of self-affine spots {$S_1, S_2$}, {$S_3, S_4$}, {$S_5, S_6$} with different spectral content but similar value of $m_2$ are displayed in Fig. 22 and are highlighted in Fig. 21 and in Table 4. Remarkably, the three foremost right lines (d) in the figure seem to continue each other. The spots {$S_5, S_6$} well illustrate this link: they in fact have the same number of modes, but different $H$ and $m_{0,h}$.

<table><caption>Table 4: Parameters for spots shown in Fig. 22 and the resulting total flux. The spots are ordered according to the increasing value of $\sqrt{m_0}$.</caption><thead><tr><th rowspan="2">Spot #</th><th rowspan="2">Q''</th><th colspan="3">Parameters A'</th><th colspan="4">Parameters I''</th></tr><tr><th>√m<sub>0</sub>/⟨r⟩</th><th>√m<sub>2</sub></th><th>√m<sub>4</sub>/⟨r⟩</th><th>√m<sub>0,h</sub></th><th>k<sub>l</sub></th><th>k<sub>s</sub></th><th>H</th></tr></thead><tbody><tr><td>S<sub>1</sub></td><td>0.189</td><td>0.100</td><td>0.507</td><td>5.10</td><td>0.100</td><td>2</td><td>16</td><td>0.5</td></tr><tr><td>S<sub>2</sub></td><td>0.213</td><td>0.143</td><td>0.515</td><td>2.84</td><td>0.141</td><td>2</td><td>8</td><td>0.5</td></tr><tr><td>S<sub>3</sub></td><td>0.510</td><td>0.0995</td><td>1.56</td><td>33.6</td><td>0.100</td><td>8</td><td>32</td><td>0.5</td></tr><tr><td>S<sub>4</sub></td><td>0.476</td><td>0.0705</td><td>1.56</td><td>61.9</td><td>0.005</td><td>8</td><td>64</td><td>0.5</td></tr><tr><td>S<sub>5</sub></td><td>0.665</td><td>0.0997</td><td>2.62</td><td>114</td><td>0.100</td><td>8</td><td>64</td><td>0.25</td></tr><tr><td>S<sub>6</sub></td><td>0.666</td><td>0.141</td><td>2.61</td><td>94.1</td><td>0.141</td><td>8</td><td>64</td><td>0.75</td></tr></tbody></table>

Finally, we would like to point out that (1) the influence of $\sqrt{m_0}$ is ultimately handled by normalization; (2) the results are rather well clustered along a simple trend line in terms of $m_2$, (3) but clearly there is a dependence on $\zeta$. At the same time, the shape of the trend is very similar to what was observed
<figure><img src="image_23.png" /><figcaption>Figure 21: Results for the renormalized total flux Eq. (30) of self-affine spots with respect to standard deviation of the radius gradient $\sqrt{m_2}$ obtained with controlled standard deviation $\sqrt{m_{0,h}} \approx \sqrt{m_0}/\langle r \rangle$.</figcaption></figure>

for multi-petal spots and we recall that $\sqrt{m_2}$ is analogous parameter to $n'$ used there. Qualitatively the slope of the normalized total flux with respect to $\sqrt{m_2}$ decreases suggesting an ultimate saturation as in multi-petal shapes (see Figs. 9 and 13). For the extra generative parameter $\zeta \in I''$ to which some dependence is observed, it should be expressed through spectral moments which could be easily measured for arbitrary shape, it will be handled in the following subsection.

## 5.5 Conductivity model

This study aims to quantify the flux transmitted through a spot of complex shape. While the numerical results encompass a broad parametric space, they are not readily comprehensible in their full scope and pose challenges for generalization. We thus make an attempt to construct a general phenomenological model relying on geometrical characteristics and inspired from the model used for flower-shaped spots.

### 5.5.1 Covariance matrix

The first simple step is the construction of a covariance $C_{ij}$ matrix of the normalized flux and all available normalized parameters $\tilde{x}_i$:

$$C_{ij} = \langle \tilde{x}_i \tilde{x}_j \rangle, \quad \tilde{x}_i = \frac{x_i - \langle x_i \rangle}{\sigma(x_i)}, \qquad (32)$$

where, as previously, $\langle x_i \rangle$ denotes the average value, and $\sigma(x_i)$ denotes its standard deviation. The covariance matrix constructed based on all available simulation data is provided in Table 5. There is a strong correlation between $\tilde{Q}', \tilde{Q}''$ and parameters $\tilde{k}_l$, $\sqrt{\tilde{m}_0}/\langle r \rangle$, $\sqrt{\tilde{m}_2}$ and $\sqrt{\tilde{m}_4}\langle r \rangle$. However, because of the strong correlation, the effect of the moment $m_4$ is hard to isolate from the effect of $m_2$. Very small correlation of the flux is found with $\tilde{\xi}, \tilde{H}$ and $\tilde{\alpha}$; slightly more correlation exists with $\tilde{\zeta}$. According
<figure><img src="image_24.png" /><figcaption>Figure 22: Comparison of pairs of self-affine spots for the same $\sqrt{m_2}$ and for different values of $\sqrt{m_0}/\langle r \rangle$. To simplify the reading of the parameters, they are equivalently displayed in Table 4.</figcaption></figure>

to the covariance matrix, the Hurst exponent seems to be negligible, which is surprising in the light of our previous results. In conclusion, we could confirm that the covariance matrix and eventual Principal Component Analysis, which access only first order correlations, present too coarse tools to determine subtle non-linear correlations. Finally, since generative parameters *I* are strongly linked to the method of spot generation, in constructing our model we will focus exclusively on geometrical parameters *A* which could be measured for an arbitrary shape.

### 5.5.2 Phenomenological model

Analyzing the obtained results, we noticed a weak logarithmic dependence of the total normalized flux on the magnification parameter $\zeta = k_s/k_l$. A relatively simple phenomenological model including this parameter could be constructed, but since this parameter is generative, it is of no help for a general case. Nevertheless, it is clear that the magnification $\zeta$ is intimately related to another geometrical parameter, known as Nayak parameter [Nayak, 1971, Yastrebov et al., 2017] $\alpha = m_0 m_4 / m_2^2$ (see covariance matrix in Table 5). Remark that from Eqs. (27) for the flower-shaped spot, the Nayak parameter is simply 1 so, consistently it does not enter the phenomenological equation for the conductivity of such simple forms Eq. (13). A rigorous link between the generative parameter $\zeta$ and the geometrical parameter $\alpha$ could be provided, see Appendix B. The concrete form of this link was not used, but we could notice that another geometrical characteristic, namely the Hurst exponent *H* is involved, regardless the results of the
<table><caption>Table 5: Covariance matrix of normalized flux and all normalized parameters according to Eq. (32).</caption><thead><tr><th></th><th>Q&#771;</th><th>Q''</th><th>&xi;&#771;</th><th>k&#771;</th><th>&zeta;</th><th>H&#771;</th><th>&radic;<span style="text-decoration: overline;">m<sub>0</sub></span>/&lang;r&rang;</th><th>&radic;<span style="text-decoration: overline;">m<sub>2</sub></span></th><th>&radic;<span style="text-decoration: overline;">m<sub>4</sub></span>/&lang;r&rang;</th><th>&alpha;&#771;</th></tr></thead><tbody><tr><th>Q&#771;</th><td>1</td><td>0.9</td><td>0.1</td><td>0.7</td><td>0.3</td><td>-0.2</td><td>0.7</td><td>0.9</td><td>0.7</td><td>0.06</td></tr><tr><th>Q''</th><td>-</td><td>1</td><td>-0.04</td><td>0.8</td><td>0.4</td><td>-0.2</td><td>0.7</td><td>0.9</td><td>0.7</td><td>0.1</td></tr><tr><th>&xi;&#771;</th><td>-</td><td>-</td><td>1</td><td>-0.2</td><td>-0.2</td><td>0.1</td><td>0.7</td><td>-0.06</td><td>-0.07</td><td>-0.08</td></tr><tr><th>k&#771;</th><td>-</td><td>-</td><td>-</td><td>1</td><td>0.06</td><td>0.0</td><td>0.4</td><td>0.7</td><td>0.4</td><td>-0.03</td></tr><tr><th>&zeta;</th><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td><td>0.0</td><td>0.06</td><td>0.5</td><td>0.5</td><td>0.9</td></tr><tr><th>H&#771;</th><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td><td>-0.2</td><td>-0.2</td><td>0.3</td><td>0.2</td></tr><tr><th>&radic;<span style="text-decoration: overline;">m<sub>0</sub></span>/&lang;r&rang;</th><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td><td>0.6</td><td>0.4</td><td>-0.05</td></tr><tr><th>&radic;<span style="text-decoration: overline;">m<sub>2</sub></span></th><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td><td>0.9</td><td>0.2</td></tr><tr><th>&radic;<span style="text-decoration: overline;">m<sub>4</sub></span>/&lang;r&rang;</th><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td><td>0.3</td></tr><tr><th>&alpha;&#771;</th><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td></tr></tbody></table>

covariance analysis. So, the ultimate set of geometrical dimensionless parameters is chosen to be:

$$ \mathcal{A}^f = \{\sqrt{m_0}/\langle r \rangle, \sqrt{m_2}, H, \alpha\} \equiv \left\{ \sigma/\langle r \rangle, \sqrt{\langle (\nabla r)^2 \rangle}, H, \alpha \right\} \quad (33) $$

The fourth moment $m_4$ does not enter explicitly in the set of parameters, only through the Nayak parameter $\alpha$ similarly to models of rough contact [Greenwood, 2006, Carbone and Bottiglione, 2008, Yastrebov et al., 2017]. Finally, we suggest the following form for the phenomenological model:

$$ Q'' = a \left[ 1 - \frac{1}{b\sqrt{m_2} + 1} \right] (1 + cH) \left\{ 1 + \frac{d}{e\alpha^f + 1} \right\} \quad (34) $$

with the core term in square brackets which is equivalent to the phenomenological law obtained for multi-petal spots, see Eq. (13). In addition, the effects of $H$ and $\alpha$ enter the equation through the product of linear and non-linear functions, respectively (normal and curly brackets). The former is the increasing function of $H$ and the latter is a decreasing function of $\alpha$. Both terms provide a slight factor adjustment: in the interval $(1, 1+c)$ for $H \in (0, 1)$, and in the interval $(1+d/(1+e), 1)$ for $\alpha \in (1, \infty)$. Due to a weak dependence on the Nayak parameter, we made an attempt to integrate it through a logarithmic dependence, like in [Yastrebov et al., 2017], but the constructed model could not fulfil the physical consistency, i.e. to always ensure positive normalized flux $Q''$ which increases monotonically for increasing $\alpha$ and $\zeta$ (see Appendix C). This physical consistency could be formulated as an inequality for the exponent parameter $f$:

$$ f \geq \frac{1-H}{2H}. $$

The issue with this bound is that it diverges for $H \to 0$. Therefore, we deliberately fixed the minimal value of the Hurst exponent that we took into consideration $H \geq 0.25$, providing the following condition for the exponent $f \geq 1.5$. A further study should be carried out to formulate a physically consistent phenomenological model for the flux for spots with lower values of the Hurst exponent.

Combining Eqs. (34) and (30), the final equation for the flux is obtained as:

$$ Q = Q_0 \left( 1 + a \frac{\sqrt{2m_0}}{\langle r \rangle} \left[ 1 - \frac{1}{b\sqrt{m_2} + 1} \right] (1 + cH) \left\{ 1 + \frac{d}{e\alpha^f + 1} \right\} \right) \quad (35) $$

The coefficients are found by the least square fit of all simulation results, see Table 6. Results of the fitting law are shown in Fig. 23 separately for two sets of simulation data: in Fig. 23(a) for the set of contact spots parametrized by Eq. (19a), and in Fig. 23(b) for those defined by Eq. (31). A relative error could be defined as:

$$ E = \frac{1}{N} \sum_{i=1}^{N} \frac{|Q_i - Q_i^{\text{fit}}|}{Q_i}, \quad (36) $$
<table><caption>Table 6: Parameters of the phenomenological model (34),(35) optimized through least square fit and resulting in relative error $E = 4.3\%$ and Pearson's correlation factor $\rho = 0.9976$.</caption><thead><tr><th>Parameter</th><th>a</th><th>b</th><th>c</th><th>d</th><th>e</th><th>f</th></tr></thead><tbody><tr><td>Value</td><td>0.968</td><td>0.255</td><td>0.0867</td><td>4.38</td><td>5.49</td><td>1.50</td></tr></tbody></table>

and for the fit coefficients the error reduces to $E = 4.3\%$. The Pearson's correlation factor for the set of identified parameters is equal to $\rho = 0.9976$.

The high concentration of data for low $\sqrt{m_2}$ might have introduced biases during the fitting process. Nonetheless, the obtained model captures well all the trends observed in our simulation results, notably, it represents well the flux of the roughest contact spots with highest values of $\sqrt{m_2}$ and $\alpha$. In summary, the obtained model could be seen as a generalization of the initial model formulated for multi-petal shapes (flower-, star- and gear-like). The ultimate model integrates the combined effects not only of standard deviation $\sigma = \sqrt{\langle\langle(r - \langle r \rangle)^2\rangle} = \sqrt{m_0}$ and $\sqrt{\langle\langle(\nabla r)^2\rangle} = \sqrt{m_2}$ but also of more subtle shape parameters such as the Nayak parameter $\alpha$ and the Hurst exponent $H$, which are related to bandwidth length and fractal dimension, respectively.

As a by product, the form of the phenomenological model (35) permits us to access the fractal limit of the self-affine spots, when the magnification $\zeta = k_s/k_l \to \infty$, then a very simple form for the limit flux could be obtained, depending only on the standard deviation of the spot and its Hurst exponent:

$$ \lim_{\zeta \to \infty} (Q) = Q_\circ \left( 1 + a \frac{\sqrt{2m_0}}{\langle r \rangle} (1 + cH) \right) = 4kU_0 \left( \langle r \rangle + a \sqrt{2m_0} (1 + cH) \right). \quad (37) $$

As a first order approximation, one could use the following value $4kU_0(\langle r \rangle + \sqrt{2m_0})$ which remains relatively accurate due to the factor $a(1 + cH)$ having minimal variation, remaining within the range $(0.968, 1.052)$. In general, this fractal limit remains speculative and could be seen as our conjecture as for the case of multi-petal shapes.

# 6 Conclusion

In establishing the bounds on the conductivity of rough contacts [Barber, 2003], Barber argued that "its greatest potential probably lies in establishing the maximum effect of neglected microscales of roughness in a solution of the contact problem for bodies with multiscale or fractal roughness." In our contribution, we focus on these "microscales" and make an attempt to assess their quantitative effect on the conductivity. If we repeat after Samuel Karlin that "the purpose of models is not to fit the data but to sharpen the questions", this study indeed permitted to sharpen few of them.

## 6.1 Flower-shaped spots and other simple forms

For simple multi-petal shapes: flower-, star- and gear-like conductive spots we could obtain the following results. In the limit of the infinite number of petals, rays, and teeth, the conductivity seems to converge to different finite limits. The bigger the area, the higher its limit, therefore the gear-like shapes have the highest and star-like shapes the lowest conductivity. We determined these limits by an extrapolation of a constructed phenomenological model, and these results should be interpreted as a first guess. Hence, the first question is whether a conductivity of such spots could be determined analytically in the limit of infinite number of petals, rays or teeth? Expectedly, these limits are bounded between the conductivity of a circle with the average radius $r_0$ and a circle with the radius equal to the maximal extent of these spots $r_0(1 + \xi)$. On the other hand, in this limit, the boundary of the conducting spot could be seen as fuzzy, with the same geometrical bounds but different "fuzziness" types, which surprisingly significantly affects the limit. The physical and mathematical limits could be different here because of radiative and
<figure><img src="image_25.png" /><figcaption>Figure 23: Simulation results for the normalized flux through self-affine spots (color markers and interpolation lines) and phenomenological prediction (smaller black markers of the same type). Upper row: initial set of data (1), lower raw: data with controlled standard deviation (2). The letters in the second series refer to those defined Fig. 21.</figcaption></figure>

eventually convective heat exchanges or because of tunneling effects. The physical conductivity should
probably hit the upper limit $Q_{up}$ defined by the conductivity of a circular spot of radius $r_0(1 + \xi)$.
## 6.2 Conductivity of self-affine spots

In terms of conductivity of self-affine random spots, based on numerous simulation results and being inspired by the phenomenological model constructed for a flower-shaped spot, we suggested a phenomenological model including four parameters: (1) mean radius, (2) its standard deviation (or the square root of the zeroth spectral moment), (3) its gradient's standard deviation (or the square root of the second spectral moment), (4) its Hurst exponent and (5) its Nayak parameter. The model is applicable in a relatively large interval of parameters and properly describes the change in flux with these geometrical parameters. It is worth noting that the model shows an interplay between the second spectral moment and a specific combination of the Hurst exponent and the Nayak parameter. The conductivity increases with the former and decreases with the latter. In the generative model employed in this study, under increasing magnification, the second spectral moment and the Nayak parameter increase in such a way that the flux is always a monotonically increasing function (by construction). Nonetheless, it is conceivable to design shapes where these two parameters are independently controlled. Consequently, the *second question* arises: could an increase in the Nayak parameter actually lead to a reduction in flux in practical scenarios? An affirmative response would intriguingly imply the existence of an optimal Nayak parameter (linked to an ideal shape) that maximizes conductivity. However, such a scenario seems rather unlikely.

Similar to our analysis of simple multi-petal shapes, the phenomenological model enabled us to determine the ultimate fractal limit for the conductivity of self-affine shapes as the magnification $\zeta$ approaches infinity. This limit depends solely on the mean radius, standard deviation, and only weakly on the Hurst exponent. However, as in our earlier findings, this identified limit should be regarded as a preliminary estimate. The mathematical question of the conductivity of self-affine shapes in the fractal limit remains open for further exploration. From a physical perspective, similarly to observations with flower-shaped spots, the diffusive nature of the boundary could provide a more practical approach to determining this limit.

## 6.3 Contact spots between rough surfaces

Concerning the conductivity of contact spots formed between randomly rough surfaces in contact, we can highlight several pertinent findings. The non-simple connectedness of these spots, characterized by non-contact areas surrounded by contact ones, does not appear to significantly affect overall conductivity. However, the complexity of their shapes undoubtedly influences this conductivity. Drawing from our analyses of relatively simpler cases, a set of parameters proves effective for estimating conductivity using the developed phenomenological models (Eqs. (34) and (35)). These parameters include (1) average radius, (2) standard deviation, (3) second spectral moment, (4) Hurst exponent, and (5) Nayak's parameter of the outer contour. But this model should be applied to realistic contact spots with caution. In most cases, such spots cannot be parametrized through a function in polar coordinates $r(\theta)$ and, in general, polar coordinates do not make much sense for complex spots (see Fig. 1). Instead, a more general parametrization using convective coordinates defined along the outer boundary is needed. In this study, by limiting ourselves to relatively simple geometrical models, we left more realistic contact shapes for future research. As a preliminary simulation for this research, we present a conductivity analysis of Koch snowflakes [Koch, 1904], with conductivity results for several initial iterations shown in Fig. 24. Qualitatively, these results align with the observed conductivity saturation at the fractal limit, a natural outcome for the physical problem of conduction.

## Acknowledgement

This work is a part of PhD project of PB which was funded by MINES Paris - PSL as "Thèse-Open". The authors are grateful to Samuel Forest and Cristian Ovalle-Rodas for thoughtful discussions and for their useful remarks during the whole PhD project of PB. The authors also would like to thank Stéphanie Chaillat for her guidance with hierarchical matrices and Basile Marchand for his help with programming
<figure><img src="image_26.png" /><figcaption>Figure 24: Normalized flux at Koch snowflake spots at several first fractal iterations and the evolution of the total flux.</figcaption></figure>

issues. In this work, a Large Language Model assisted the authors with polishing the text and, to a little extent, with the code improvement.

## A Geometrical characteristics of flower-, star- and gear-like shapes

Geometrical characteristics (perimeter, area and compactness) are summarized in Table 7 and Fig. 25 for flower-, star- and gear-like spots. We remind that the mean radius is denoted by $r_0$, and half-length of petals (stars or teeth) is equal to $r_0\xi$.

## B Geometrical characteristics of self-affine spots

In Fig. 26, the analytical form for the standard deviation of self-affine spots $\sigma_s$ Eq. (24) is compared with the one evaluated over 1000 generated spots for each combination of generative parameters: $\xi = 0.05$, $k_l \in \{2, 4, 8, 16, 32, 64, 128\}$, $\zeta \in \{4, 8, 16, 32\}$ and $H \in [0.2, 0.8]$. These results are quite sensitive to $k_l$, but the maximal relative error is 0.05% for $k_l = 32$ and $\zeta = 8$.

The second $m_2$ and the fourth $m_4$ moments have been computed for the same set of generative parameters over the same 1000 spots. These moments could be computed in three different ways. First, the discretized contour geometry could be used to evaluate these moments $m_p^D$ using Eq. (25). The discretization consists of splitting the contour in $N = \max\{10000, 100k_s\}$ straight segments with $d\theta = 2\pi/N$ and
<table><caption>Table 7: Geometrical characteristics (area, perimeter and compactness) for multi-petal shapes: flower, star and gear; $E(in') = E(inξ) = \int_{0}^{\pi/2} \sqrt{1 + (nξ\sin(\theta))^2} d\theta$ is the complete elliptic integral of the second kind and $i$ is the imaginary unit.</caption><thead><tr><th>Characteristics</th><th>Flower</th><th>Star</th><th>Gear</th></tr></thead><tbody><tr><td>Area A</td><td>πr<sub>0</sub><sup>2</sup>(1 + ξ<sup>2</sup>/2)</td><td>nr<sub>0</sub><sup>2</sup>(1 - ξ<sup>2</sup>) sin(π/n)</td><td>πr<sub>0</sub><sup>2</sup>(1 + ξ<sup>2</sup>)</td></tr><tr><td>Perim. P</td><td>r<sub>0</sub>E(in')</td><td>2√2nr<sub>0</sub>√(1 + ξ<sup>2</sup> - (1 - ξ<sup>2</sup>)cos(π/n))</td><td>2πr<sub>0</sub>(1 + 2n'/π)</td></tr><tr><td>Comp. C = √A/P</td><td>√π√(1 + ξ<sup>2</sup>/2) / (4E(in'))</td><td>√(1 - ξ<sup>2</sup>) sin(π/n) / (2√2n√(1 + ξ<sup>2</sup> - (1 - ξ<sup>2</sup>)cos(π/n)))</td><td>√(1 + ξ<sup>2</sup>) / (2√π(1 + 2n'/π))</td></tr></tbody></table>

<figure><img src="image_27.png" /><figcaption>Figure 25: Perimeter, area and compactness for the different shapes of contact spot</figcaption></figure>
evaluation gradient and Laplacian as

$$ \nabla r_i = \frac{2(r_{i+1} - r_i)}{(r_{i+1} + r_i)d\theta}, \quad \Delta r_i = \frac{4(r_{i+1} - 2r_i + r_{i-1})}{(r_{i+1} + 2r_i + r_{i-1})^2 d\theta^2}, \qquad (38) $$

where $r_i = r(id\theta)$, $i = 1, N$. This method was used throughout the paper. Second, the moments could be approximated by discrete sums of all mode contributions as in Eq. (26):

$$ m_p^S = \frac{(r_0 \xi)^2}{2} \sum_{k_l}^{k_s} k^p \left( \frac{k}{k_l} \right)^{-(1+2H)} . \qquad (39) $$

This method is however valid only for relatively small values of $\xi$ because it ignores the exponential transformation (21). Third, for sufficiently large values of $k_l$, these discrete sums could be turned into integrals with wavenumber $k$ becoming a continuum variable of integration:

$$ m_p^C = \frac{1}{2} \int_{k_l}^{\zeta k_l} (r(k))^2 k^p dk, \qquad (40) $$

where $r(k) = \xi r_0(k/k_l)^{-H-0.5}$. Analytical formulas derive from the development for the moments $m_0^C$, $m_2^C$ and $m_4^C$, as follows.

$$ m_0^C = -\frac{(r_0 \xi)^2 k_l}{4H} (\zeta^{-2H} - 1), \quad m_2^C = -\frac{(r_0 \xi)^2 k_l^3}{2(2-2H)} (\zeta^{2-2H} - 1), \quad m_4^C = -\frac{(r_0 \xi)^2 k_l^5}{2(4-2H)} (\zeta^{4-2H} - 1). \qquad (41) $$

For high values of $k_l$ the discrete spectrum is closer to a continuous one, and thus the spectral moments can be deduced from Eq. (26) as detailed by Nayak [Nayak, 1971]. These analytical values are compared with the numerically evaluated ones in Fig.27 and Fig.28. The maximum deviation is of only 0.7 % for the $m_2$ for parameters $k_l = 2, \zeta = 4$ and $H = 0.2$. However, for $m_4$ an average discrepancy of 10 % is observed and could raise to as much as 27 % in certain instances. Nevertheless, in all results presented in the paper only actual values of the moments and of their combinations were used.

The Nayak's parameter is determined using the moments $m_0, m_2$ and $m_4$, $\alpha = m_0 m_4 / m_2^2$. The three models mentioned above could be used to compute the Nayak parameter as $\alpha^D, \alpha^S$ and $\alpha^C$, respectively. The average values computed over a set of 1000 spots are compared in Fig. 29. This analysis demonstrates that a continuum model could be successfully used in practical applications. In the limit of high magnification $\zeta$, the second moment and the Nayak parameter scale as $m_2^C \sim \zeta^{2-2H}$ and $\alpha^C \sim \zeta^{2H}$.

## C Physical consistency of the phenomenological model

The phenomenological model of flux Eqs. (34), (35) exhibits an increasing behavior with respect of $m_2$, but decreases with $\alpha$. Nevertheless, from general physical considerations, we conjecture that the flux should be a monotonically non-decreasing function of the magnification $\zeta$. So, we require that the derivative of flux $Q$ with respect to magnification $\zeta$ remains non-negative:

$$ \frac{\partial Q}{\partial \zeta} = \frac{\partial Q}{\partial m_2} \frac{\partial m_2}{\partial \zeta} + \frac{\partial Q}{\partial \alpha} \frac{\partial \alpha}{\partial \zeta} \geq 0 \qquad (42) $$

The terms $\partial Q / \partial m_2$ and $\partial Q / \partial \alpha$ could be readily derived from Eq. (35).

$$ \frac{\partial Q}{\partial m_2} = a \frac{b}{2\sqrt{m_2}(b\sqrt{m_2}+1)^2} (1+cH) \left(1+\frac{d}{e\alpha f+1}\right) \qquad (43a) $$

$$ \frac{\partial Q}{\partial \alpha} = a \left\{ 1 - \frac{b}{b\sqrt{m_2} + 1} \right\} (1 + cH) \frac{-def\alpha^{f-1}}{(e\alpha f + 1)^2} \qquad (43b) $$
<figure><img src="image_28.png" /><figcaption>Figure 26: Results of standard deviation by spectral and sample analysis, with $k_l = \{2, 4, 8, 16, 32, 64, 128\}$ in both figures: (a) $\zeta = \{4, 8, 16, 32\}$ and $H = 0.2$; (b) $H = \{0.2, 0.4, 0.6, 0.8\}$ and $\zeta = 4$.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Figure 27: Results of mean square gradient by spectral and sample analysis, with $k_l = \{2, 4, 8, 16, 32, 64, 128\}$ in both figures: (a) $\zeta = \{4, 8, 16, 32\}$ and $H = 0.2$; (b) $H = \{0.2, 0.4, 0.6, 0.8\}$ and $\zeta = 4$.</figcaption></figure>

The derivatives of $m_2$ and $\alpha$ with respect to $\zeta$ could be found from Eq.(41), resulting in the following asymptotic forms:

$$ \frac{\partial Q}{\partial m_2} \frac{\partial m_2}{\partial \zeta} \sim \frac{1}{\sqrt{m_2(\sqrt{m_2}+1)^2}} \frac{\partial m_2}{\partial \zeta} \sim \zeta^{-2+H} \qquad (44a) $$

$$ \frac{\partial Q}{\partial \alpha} \frac{\partial \alpha}{\partial \zeta} \sim \frac{\alpha^{f-1}}{(\alpha^{f}+1)^2} \frac{\partial \alpha}{\partial \zeta} \sim \zeta^{-2fH-1} \qquad (44b) $$

These expressions enable us to define a constraint criterion to ensure the derivative of the flux law with respect to $\zeta$ remains non-negative for all $\zeta$. The exponent of $\zeta$ in Eq. (44a) must be lower than the one in Eq.(44a), resulting in the following inequality that the exponent $f$ should satisfy:

$$ f \geq \frac{1-H}{2H} \qquad (45) $$
<figure><img src="image_30.png" /><figcaption>Figure 28: Results of mean square Laplacian, with $k_l = \{2, 4, 8, 16, 32, 64, 128\}$ in both figures: (a) $\zeta = \{4, 8, 16, 32\}$ and $H = 0.2$; (b) $H = \{0.2, 0.4, 0.6, 0.8\}$ and $\zeta = 4$.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Figure 29: Comparison of different models which could be used to evaluate spectral moments and the Nayak parameter: (1) discretized geometrical evaluation, (2) discrete sum for the generative function, (3) continuous version of this discrete sum.</figcaption></figure>

The problem with this constraint is that it results in too high values of $f$ for small $H$ and, ultimately, it diverges for $H \to 0$. In the current study we set the minimal value of the Hurst exponent to $H = 0.25$ thus resulting in $f \ge 1.5$. The results for derivatives using the continuous expressions for $m_2$ and $\alpha$ are presented Fig. 30 for $H = 0.25$, $k_l = 8$, and $\xi = 0.05$, and the fitting parameters shown in Table 6. The
two derivative terms are distinguished: one positive, as given by Eq. (43a), and the other negative, as given by Eq. (43a). The full derivative remains positive, however, thus keeping the required assumptions true, even for value of ζ significantly far from the initial set of parameter. The absolute values for these derivatives are also depicted in inset in log-log scale, showing similar power-laws of the two competing derivatives (in dots) for high values of ζ.

<figure><img src="image_32.png" /><figcaption>Figure 30: Derivatives of the flux phenomenological model in function of the magnification ζ for $k_l = 8$, $H = 0.25$, and $\xi = 0.05$ for parameters from Table 6.</figcaption></figure>

## References

[Afferrante et al., 2012] Afferrante, L., Carbone, G., and Demelio, G. (2012). Interacting and coalescing hertzian asperities: A new multasperity contact model. *Wear*, 278:28–33.

[Amuzuga et al., 2016] Amuzuga, K. V., Chaise, T., Duval, A., and Nelias, D. (2016). Fully coupled resolution of heterogeneous elastic-plastic contact problem. *Journal of Tribology*, 138(2):021403.

[Archard, 1953] Archard, J. F. (1953). Contact and rubbing of flat surfaces. *Journal of Applied Physics*, 24(8):981–988.

[Archard, 1957] Archard, J. F. (1957). Elastic deformation and the laws of friction. *Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences*, 243(1233):190–205.

[Barber, 2003] Barber, J. R. (2003). Bounds on the electrical resistance between contacting elastic rough bodies. *Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences*, 459(2029):53–66.

[Baydoun et al., 2020] Baydoun, S., Arnaud, P., and Fouvry, S. (2020). Modelling adhesive wear extension in fretting interfaces: An advection-dispersion-reaction contact oxygenation approach. *Tribology International*, 151:106490.

[Baydoun et al., 2022] Baydoun, S., Arnaud, P., and Fouvry, S. (2022). Explicit formulations of adhesive wear extension in fretting interfaces applying the contact oxygenation concept. *Wear*, 488:204147.
[Beguin, 2024] Beguin, P. (2024). *Conductivity of Contact Interfaces: A Multi-Scale Study*. PhD thesis, MINES Paris - PSL, Paris, France.

[Beguin and Yastrebov, 2024] Beguin, P. and Yastrebov, V. A. Fast-bem code for poisson's equation on a half-space. https://github.com/vyastreb/HBEM, version 1.0.

[Beguin and Yastrebov, 2024] Beguin, P. and Yastrebov, V. A. (2024). Supplementary material to paper. Zenodo. doi:10.5281/zenodo.10200997.

[Besson and Foerch, 1997] Besson, J. and Foerch, R. (1997). Large scale object-oriented finite element code design. *Computer Methods in Applied Mechanics and Engineering*, 142(1-2):165-187.

[Blok, 1963] Blok, H. (1963). The flash temperature concept. *Wear*, 6(6):483-494.

[Bonnet, 1999] Bonnet, M. (1999). *Boundary Integral Equation Methods for Solids and Fluids*. Wiley.

[Boussinesq, 1885] Boussinesq, J. (1885). *Application des potentiels à l'étude de l'équilibre et du mouvement des solides élastiques*. Gauthier-Villars.

[Bowden et al., 2001] Bowden, F. P., Bowden, F. P., and Tabor, D. (2001). *The Friction and Lubrication of Solids*. Oxford University Press.

[Burghold et al., 2015] Burghold, E. M., Frekers, Y., and Kneer, R. (2015). Determination of time-dependent thermal contact conductance through ir-thermography. *International Journal of Thermal Sciences*, 98:148-155.

[Bush et al., 1975] Bush, A. W., Gibson, R. D., and Thomas, T. R. (1975). The elastic contact of a rough surface. *Wear*, 35(1):87-111.

[Campaná and Müser, 2006] Campaná, C. and Müser, M. H. (2006). Practical green's function approach to the simulation of elastic semi-infinite solids. *Physical Review B*, 74(7):075420.

[Carbone and Bottiglione, 2008] Carbone, G. and Bottiglione, F. (2008). Asperity contact theories: Do they predict linearity between contact area and load? *Journal of the Mechanics and Physics of Solids*, 56(8):2555-2572.

[Carslaw and Jaeger, 1959] Carslaw, H. S. and Jaeger, J. C. (1959). *Conduction of Heat in Solids*. Oxford University Press.

[Chaillat et al., 2017] Chaillat, S., Desiderio, L., and Ciarlet, P. (2017). Theory and implementation of h-matrix based iterative and direct solvers for helmholtz and elastodynamic oscillatory kernels. *Journal of Computational Physics*, 351:165-186.

[Collins, 1963] Collins, W. D. (1963). On the solution of some axisymmetric boundary value problems by means of integral equations: VIII. potential problems for a circular annulus. *Proceedings of the Edinburgh Mathematical Society*, 13(3):235-246.

[Cooke, 1963] Cooke, J. C. (1963). Some further triple integral equation solutions. *Proceedings of the Edinburgh Mathematical Society*, 13(4):303-316.

[Cooper et al., 1969] Cooper, M. G., Mikic, B. B., and Yovanovich, M. M. (1969). Thermal contact conductance. *International Journal of heat and mass transfer*, 12(3):279-300.

[Cuthrell and Tipping, 1973] Cuthrell, R. E. and Tipping, D. W. (1973). Electric contacts. ii. mechanics of closure for gold contacts. *Journal of Applied Physics*, 44(10):4360-4365.

[Dapp et al., 2012] Dapp, W. B., Lücke, A., Persson, B. N. J., and Müser, M. H. (2012). Self-affine elastic contacts: percolation and leakage. *Physical Review Letters*, 108(24):244301.
[Dodds and Robson, 1973] Dodds, C. J. and Robson, J. D. (1973). The description of road surface roughness. *Journal of Sound and Vibration*, 31(2):175-183.

[Fabrikant, 1993] Fabrikant, V. I. (1993). Dirichlet problem for an annular disk. *Zeitschrift für angewandte Mathematik und Physik ZAMP*, 44(2):333–347.

[Fieberg and Kneer, 2008] Fieberg, C. and Kneer, R. (2008). Determination of thermal contact resistance from transient temperature measurements. *International Journal of Heat and Mass Transfer*, 51(5-6):1017–1023.

[Frérot et al., 2019] Frérot, L., Bonnet, M., Molinari, J.-F., and Anciaux, G. (2019). A fourier-accelerated volume integral method for elastoplastic contact. *Computer Methods in Applied Mechanics and Engineering*, 351:951–976.

[Gao and Bower, 2006] Gao, Y.-F. and Bower, A. F. (2006). Elastic-plastic contact of a rough surface with weierstrass profile. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 462(2065):319–348.

[Goldsby and Tullis, 2011] Goldsby, D. L. and Tullis, T. E. (2011). Flash heating leads to low frictional strength of crustal rocks at earthquake slip rates. *Science*, 334(6053):216–218.

[Grasedyck, 2005] Grasedyck, L. (2005). Adaptive recompression of matrices for bem. *Computing*, 74(3):205–223.

[Grasedyck and Hackbusch, 2003] Grasedyck, L. and Hackbusch, W. (2003). Construction and arithmetics of h-matrices. *Computing*, 70(4):295–334.

[Greenwood, 1966] Greenwood, J. A. (1966). Constriction resistance and the real area of contact. *British Journal of Applied Physics*, 17(12):1621.

[Greenwood, 2006] Greenwood, J. A. (2006). A simplified elliptic model of rough surface contact. *Wear*, 261(2):191–200.

[Greenwood et al., 1984] Greenwood, J. A., Johnson, K. L., and Matsubara, E. (1984). A surface roughness parameter in hertz contact. *Wear*, 100(1-3):47–57.

[Greenwood and Tripp, 1967] Greenwood, J. A. and Tripp, J. H. (1967). The Elastic Contact of Rough Spheres. *Journal of Applied Mechanics*, 34(1):153-159.

[Greenwood and Williamson, 1966] Greenwood, J. A. and Williamson, J. B. P. (1966). Contact of nominally flat surfaces. *Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences*, 295(1442):300–319.

[Hackbusch, 2015] Hackbusch, W. (2015). *Hierarchical Matrices: Algorithms and Analysis*. Springer.

[Holm, 1957] Holm, R. (1957). *Electric Contacts: Theory and Application*. Springer.

[Hyun et al., 2004] Hyun, S., Pei, L., Molinari, J.-F., and Robbins, M. O. (2004). Finite-element analysis of contact between elastic self-affine surfaces. *Physical Review E*, 70(2):026117.

[Jackson et al., 2008] Jackson, R. L., Bhavnani, S. H., and Ferguson, T. P. (2008). A multiscale model of thermal contact resistance between rough surfaces. *Journal of Heat Transfer*, 130(8):081301.

[Jacobson and Hogmark, 2009] Jacobson, S. and Hogmark, S. (2009). Surface modifications in tribo-logical contacts. *Wear*, 266(3-4):370–378.
[Jensen et al., 2005] Jensen, B. D., Huang, K., Chow, L. L.-W., and Kurabayashi, K. (2005). Low-force contact heating and softening using micromechanical switches in diffusive-ballistic electron-transport transition. *Applied Physics Letters*, 86(2):023507.

[Kalker, 1977] Kalker, J. J. (1977). Variational principles of contact elastostatics. *IMA Journal of Applied Mathematics*, 20(2):199-219.

[Kempers et al., 2009] Kempers, R., Kolodner, P., Lyons, A., and Robinson, A. J. (2009). A high-precision apparatus for the characterization of thermal interface materials. *Review of Scientific Instruments*, 80(9):095111.

[Koch, 1904] Koch, H. V. (1904). Sur une courbe continue sans tangente, obtenue par une construction géométrique élémentaire. *Arkiv for Matematik, Astronomi och Fysik*, 1:681-704.

[Kogut and Komvopoulos, 2003] Kogut, L. and Komvopoulos, K. (2003). Electrical contact resistance theory for conductive rough surfaces. *Journal of Applied Physics*, 94(5):3153-3162.

[Kumano et al., 1994] Kumano, H., Sawa, T., and Hirose, T. (1994). Mechanical behavior of bolted joints under steady heat conduction. *Journal of Pressure Vessel Technology*, 116(1):42-48.

[Lambert and Fletcher, 1997] Lambert, M. A. and Fletcher, L. S. (1997). Thermal contact conductance of spherical rough metals. *Journal of Heat Transfer*, 119(4):684-690.

[Love, 1976] Love, E. R. (1976). Inequalities for the capacity of an electrified conducting annular disc. *Proceedings of the Royal Society of Edinburgh Section A: Mathematics*, 74:257-270.

[Majumdar and Bhushan, 1991] Majumdar, A. and Bhushan, B. (1991). Fractal model of elastic-plastic contact between rough surfaces. *Journal of Tribology*, 113(1):1-11.

[Majumdar and Tien, 1990] Majumdar, A. and Tien, C. L. (1990). Fractal characterization and simulation of rough surfaces. *Wear*, 136(2):313-327.

[Manners and Greenwood, 2006] Manners, W. and Greenwood, J. A. (2006). Some observations on persson's diffusion theory of elastic contact. *Wear*, 261(5-6):600-610.

[Mikrajuddin et al., 1999] Mikrajuddin, A., Shi, F. G., Kim, H. K., and Okuyama, K. (1999). Size-dependent electrical constriction resistance for contacts of arbitrary size: from sharvin to holm limits. *Materials Science in Semiconductor Processing*, 2(4):321-327.

[Murashov and Panin, 2015] Murashov, M. V. and Panin, S. D. (2015). Numerical modelling of contact heat transfer problem with work hardened rough surfaces. *International Journal of Heat and Mass Transfer*, 90:72-80.

[Nakamura, 1993] Nakamura, M. (1993). Constriction resistance of conducting spots by the boundary element method. *IEEE transactions on components, Hybrids, and Manufacturing technology*, 16(3):339-343.

[Nakamura and Minowa, 1986] Nakamura, M. and Minowa, I. (1986). Computer simulation for the conductance of a contact interface. *IEEE transactions on components, Hybrids, and Manufacturing technology*, 9(2):150-155.

[Nayak, 1971] Nayak, P. R. (1971). Random process model of rough surfaces. *ASME Journal of Lubrication Technology*, 93:398-407.

[Pei et al., 2005] Pei, L., Hyun, S., Molinari, J. F., and Robbins, M. O. (2005). Finite element modeling of elasto-plastic contact between rough surfaces. *Journal of the Mechanics and Physics of Solids*, 53(11):2385-2409.
[Pérez-Ràfols et al., 2023] Pérez-Ràfols, F., Van Dokkum, J. S., and Nicola, L. (2023). On the interplay between roughness and viscoelasticity in adhesive hysteresis. *Journal of the Mechanics and Physics of Solids*, 170:105079.

[Persson, 2001] Persson, B. N. J. (2001). Theory of rubber friction and contact mechanics. *The Journal of Chemical Physics*, 115(8):3840–3861.

[Polonsky and Keer, 1999] Polonsky, I. A. and Keer, L. M. (1999). A numerical method for solving rough contact problems based on the multi-level multi-summation and conjugate gradient techniques. *Wear*, 231(2):206–219.

[Putignano et al., 2015] Putignano, C., Carbone, G., and Dini, D. (2015). Mechanics of rough contacts in elastic and viscoelastic thin layers. *International Journal of Solids and Structures*, 69:507–517.

[Rice, 2006] Rice, J. R. (2006). Heating and weakening of faults during earthquake slip. *Journal of Geophysical Research: Solid Earth*, 111:B05311.

[Richardson, 1911] Richardson, L. F. (1911). The approximate arithmetical solution by finite differences with an application to stresses in a masonry dam. *Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character*, 210:307–357.

[Sadowski and Stupkiewicz, 2010] Sadowski, P. and Stupkiewicz, S. (2010). A model of thermal contact conductance at high real contact area fractions. *Wear*, 268(1-2):77–85.

[Sano, 1985] Sano, Y. (1985). Effect of space angle on constriction resistance and contact resistance for a point contact. *Journal of Applied Physics*, 58(7):2651–2654.

[Slade, 2017] Slade, P. G. (2017). *Electrical Contacts: Principles and Applications*. CRC press.

[Smythe, 1951] Smythe, W. R. (1951). The capacitance of a circular annulus. *Journal of Applied Physics*, 22(12):1499–1501.

[Sneddon, 1995] Sneddon, I. N. (1995). *Fourier Transforms*. Courier Corporation.

[Stanley and Kato, 1997] Stanley, H. M. and Kato, T. (1997). An fft-based method for rough surface contact. *Journal of tribology*, 119(3):481–485.

[Tada et al., 2000] Tada, H., Paris, P. C., and Irwin, G. R. (2000). *The Stress Analysis of Cracks Handbook*. ASME, 3rd edition.

[Thomas and Probert, 1970] Thomas, T. R. and Probert, S. D. (1970). Establishment of contact parameters from surface profiles. *Journal of Physics D: Applied Physics*, 3(3):277.

[Tu, 2003] Tu, K.-N. (2003). Recent advances on electromigration in very-large-scale-integration of interconnects. *Journal of Applied Physics*, 94(9):5451–5473.

[Watanabe, 1986] Watanabe, Y. (1986). New instrument for measuring contact resistance developed for studying electrical contact phenomena. *Wear*, 112(1):1–15.

[Yamashita et al., 2015] Yamashita, F., Fukuyama, E., Mizoguchi, K., Takizawa, S., Xu, S., and Kawakata, H. (2015). Scale dependence of rock friction at high work rate. *Nature*, 528(7581):254–257.

[Yastrebov, 2019] Yastrebov, V. A. (2019). The elastic contact of rough spheres investigated using a deterministic multi-asperity model. *Journal of Multiscale Modelling*, 10(01):1841002.
[Yastrebov et al., 2015] Yastrebov, V. A., Anciaux, G., and cois Molinari, J.-F. (2015). From infinitesimal to full contact between rough surfaces: evolution of the contact area. *International Journal of Solids and Structures*, 52:83–102.

[Yastrebov et al., 2017] Yastrebov, V. A., Anciaux, G., and Molinari, J.-F. (2017). The role of the roughness spectral breadth in elastic contact of rough surfaces. *Journal of the Mechanics and Physics of Solids*, 107:469–493.

[Yastrebov et al., 2011] Yastrebov, V. A., Durand, J., Proudhon, H., and Cailletaud, G. (2011). Rough surface contact analysis by means of the finite element method and of a new reduced model. *Comptes Rendus Mécanique*, 339(7-8):473–490.

[Z-set, 2023] Z-set (2023). Z-set: Non-linear material and structure analysis suite. http://www.zset-software.com.

[Zienkiewicz and Taylor, 2000] Zienkiewicz, O. C. and Taylor, R. L. (2000). *The Finite Element Method: The Basis*, volume 1. Butterworth-Heinemann, Oxford, 5th edition.

[Zou et al., 2008] Zou, M., Yu, B., Cai, J., and Xu, P. (2008). Fractal model for thermal contact conductance. *Journal of Heat Transfer*, 130(10):101301.