# 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 multi-scale 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.

###### Contents

  1. 1 Introduction
  2. 2 Methods
  3. 3 Conductivity of a not simply connected spot
  4. 4 Conductivity of flower-shaped spots
    1. 4.1 Conductivity results
    2. 4.2 Alternative ”multi-petal” spots
  5. 5 Conductivity of self-affine spots
    1. 5.1 Geometry of spots
    2. 5.2 Geometrical characteristics
    3. 5.3 Normalization and result
    4. 5.4 Results with renormalized standard deviation
    5. 5.5 Conductivity model
      1. 5.5.1 Covariance matrix
      2. 5.5.2 Phenomenological model
  6. 6 Conclusion
    1. 6.1 Flower-shaped spots and other simple forms
    2. 6.2 Conductivity of self-affine spots
    3. 6.3 Contact spots between rough surfaces
  7. A Geometrical characteristics of flower-, star- and gear-like shapes
  8. B Geometrical characteristics of self-affine spots
  9. C  Physical consistency of the phenomenological model

##  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/transformations,
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 multi-scale 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 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 thermo-electrical 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 with 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-/electro-mechanical problems is worth
investigation.

The growing demand for micro-electric devices [Tu, 2003] has spurred an
increasing need for electrical contact models at the micro-scale. 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 configurations [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 with 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
multi-scale 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 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.

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 $\bm{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]. Hierarchical matrices allow to
overcome main drawbacks of the classical BEM requiring storage and resolution
of linear system of equations with full matrices. The open-source Python code
is shared under BSD3 license in [Beguin and Yastrebov, 2023a] along with the
details of its implementation and a short documentation. Constant
interpolation triangular elements 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]. For the FEM, in addition to mesh-extrapolation, the
problem was solved on two domains of various size and the extrapolated result
obtained with the finest mesh was also extrapolated to obtain an accurate
result.

##  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_{\circ}=4kaU_{0},$ (1)

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^{\circ}_{n}(r)=\bm{j}\cdot\bm{n}=\frac{2U_{0}k}{\pi\sqrt{a^{2}-r^{2}}},$
(2)

where $\bm{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\leq r\leq 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. 2 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_{\circ}$, 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{\tiny
Love}}}{Q_{\circ}}=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})$
(3)

$\displaystyle\xi=0.2$$\displaystyle\xi=0.5$$\displaystyle\xi=0.8$

Figure 2:  Geometry of the annular contact spot for various ratios of the
internal to external radii $\xi=\\{0.2,0.5,0.8\\}$.

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. The geometry of the axially symmetric problem is
restricted to a square plane of length and width $R\gg a$. The mesh size on
the outer edge must be much smaller compared to $a$, i.e. $h\ll a$ and near
the internal edge $h\ll b$. 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. 3 and compared with the analytical solution for the
circular spot (2). 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, 2023b]). 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. Hakamura [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. 4 with the flux
normalized by the one of a circular spot $Q_{\circ}$ (see Eq. (1)). 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}\left(1-\alpha k^{\beta}\right),$ (4)

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. (3)). 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 slighly
reduce the flux. Nevertheless, the first terms in Eq. (3) 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 extentricity with respect to the
centre. 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.

$\displaystyle{0.0}$$\displaystyle{0.2}$$\displaystyle{0.4}$$\displaystyle{0.6}$$\displaystyle{0.8}$$\displaystyle{1.0}$$\displaystyle
r^{\prime}$$\displaystyle{0}$$\displaystyle{1}$$\displaystyle{2}$$\displaystyle{3}$$\displaystyle{4}$$\displaystyle{5}$$\displaystyle
j_{n}$$\displaystyle
j^{\circ}_{n}(r^{\prime})$$\displaystyle\xi$$\displaystyle 0.2$$\displaystyle
0.4$$\displaystyle 0.6$$\displaystyle 0.8$

Figure 3:  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.

$\displaystyle{0.0}$$\displaystyle{0.2}$$\displaystyle{0.4}$$\displaystyle{0.6}$$\displaystyle{0.8}$$\displaystyle{1.0}$$\displaystyle\xi$$\displaystyle
1$$\displaystyle 0.95$$\displaystyle 0.9$$\displaystyle 0.85$$\displaystyle
0.8$$\displaystyle Q^{\prime}=Q(\xi)/Q_{\circ}$$\displaystyle\xi_{lim}=0.2$FEM
resultsNakamura’s results$\displaystyle Q^{\prime}_{\text{Love}}$
eq.(3)$\displaystyle Q^{\prime}_{fit}$
eq.(4)$\displaystyle{10^{-2}}$$\displaystyle{10^{-1}}$$\displaystyle{10^{0}}$$\displaystyle\xi$$\displaystyle{10^{-6}}$$\displaystyle{10^{-5}}$$\displaystyle{10^{-4}}$$\displaystyle{10^{-3}}$$\displaystyle{10^{-2}}$$\displaystyle{10^{-1}}$$\displaystyle
1-Q^{\prime}=1-Q(\xi)/Q_{\circ}$$\displaystyle\xi_{lim}$

Figure 4: Total flux of an annulus as a function of internal to external
radius ratio $\xi$: finite element results (circles), least squares fit of an
offsetted power-law for the asymptotic solution $Q^{\prime}_{\text{fit}}$ for
$\xi<0.2$ and an approximate analytical solutions by [Love, 1976]
$Q_{\text{\tiny Love}}/Q_{\circ}$. (a) – normalized total flux
$Q^{\prime}(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.

##  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}\left(1+\xi\cos(n\theta)\right),\quad
r_{1}<r_{0}$ (5)

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. 5. 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_{\min}(r_{0})=Q_{\circ}=4kr_{0}U_{0}$, $Q_{\text{\tiny
up}}=Q_{\circ}(r_{0}+r_{1})=4kr_{0}(1+\xi)U_{0}$,

$1\leq\frac{Q(r_{0},\xi,n)}{Q_{\min}}\leq 1+\xi$

$\displaystyle n=4,\xi=0.1$$\displaystyle r_{0}$$\displaystyle
n=4,\xi=0.2$$\displaystyle r_{0}(1+\xi)$$\displaystyle
n=7,\xi=0.1$$\displaystyle n=7,\xi=0.2$

Figure 5: Examples of flower-shaped spots with $\xi=\\{0.1,0.2\\}$ and the
number of petals $n=\\{4,7\\}$.

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

$\displaystyle A$ $\displaystyle=\pi
r_{0}^{2}\left(1+\frac{\xi^{2}}{2}\right),$ (6a) $\displaystyle P$
$\displaystyle=r_{0}\int\limits_{0}^{2\pi}\sqrt{1+\left(n\xi\sin(\theta)\right)^{2}}d\theta,$
(6b)

We can notice that the integrand depends on the dimensionless parameter
$n^{\prime}=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^{\prime})=\frac{\sqrt{A}}{P}=\sqrt{\pi}\frac{\sqrt{1+\xi^{2}/2}}{4E(in^{\prime})},$
(7)

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^{\prime}$ only, for $n^{\prime}\to 0$, $E(in^{\prime})\to\pi/2$ and the
compactness tends to the maximal value, i.e. the compactness of a circle $C\to
C_{\circ}=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. 6 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 centre. 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 neighbouring petals; a similar trend was
observed for the annular spot for small internal radii.

Figure 6: Simulation results for the normal flux for flower-shaped spots with
$n={4,7,10}$ and $\xi=0.1$.

In total, $97$ simulations111Without counting simulations used to employ
Richardson extrapolation. 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^{\prime}$ is presented in Fig. 7. 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, offsetted by the minimal flux $Q-Q_{\circ}$ 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_{\text{\tiny
up}}-Q_{\circ}=4kr_{0}\xi U_{0}=\xi Q_{\circ}$. This normalization results in
a universal curve for the total normalized flux for any $r_{0}$ and $\xi$:

$Q^{\prime}=\frac{Q-Q_{\circ}}{Q_{\text{\tiny
up}}-Q_{\circ}}=\frac{Q-Q_{\circ}}{\xi Q_{\circ}}$ (8)

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

$Q^{\prime}_{\mathrm{fit}}(n^{\prime})=a\left(1-\frac{1}{bn^{\prime}+1}\right),\quad
0<a\leq 1.$ (9)

This fit function is plotted in Fig. 7 along with FEM and BEM simulation
results. The coefficients determined by least squares fit are presented in
Table 1. The slope at $n^{\prime}=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. (8) and (9), 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),$ (10)

where $Q_{\circ}=4kr_{0}U_{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^{\prime}\to\infty}(Q^{\prime}_{\mathrm{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^{\prime}$ 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
folllowing limit flux:

$\lim\limits_{n\xi\to\infty}Q\approx Q_{\circ}(1+0.923\xi)<Q_{\text{\tiny
up}}=Q_{\circ}(1+\xi).$ (11)

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_{\text{\tiny 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)\quad\Leftrightarrow\quad
n^{\prime}=n\xi=\left[\frac{x}{b(a-x)}\right]\quad\Rightarrow\quad
n\xi\approx\left[\frac{x}{0.326(0.923-x)}\right],$ (12)

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.

| Parameters | Coefficients
---|---|---
Simulation | $\xi$ | $n\in$ | $a$ | $b$ | $ab$
FEM | $0.1$ | $[1,150]$ | 0.928 | 0.327 | 0.304
FEM | $0.2$ | $[1,100]$ | 0.923 | 0.326 | 0.301
BEM | $0.1$ | $[1,256]$ | 0.923 | 0.326 | 0.301
Table 1: Least squares fit for coefficients of Eq. (9) for the sets of results
of flower-shaped spot obtained using FEM and BEM simulations.

$\displaystyle 0$$\displaystyle 5$$\displaystyle 10$$\displaystyle
15$$\displaystyle 20$$\displaystyle 25$$\displaystyle
n^{\prime}$$\displaystyle 0$$\displaystyle 0.2$$\displaystyle
0.4$$\displaystyle 0.6$$\displaystyle 0.8$$\displaystyle 1$$\displaystyle
Q^{\prime}=(Q-Q_{\circ})/(Q_{up}-Q_{\circ})$FEM, $\displaystyle\xi=0.1$FEM,
$\displaystyle\xi=0.2$BEM, $\displaystyle\xi=0.1$$\displaystyle
0$$\displaystyle 0.5$$\displaystyle 1$$\displaystyle 1.5$$\displaystyle
2$$\displaystyle 2.5$$\displaystyle 0$$\displaystyle 0.1$$\displaystyle
0.2$$\displaystyle 0.3$$\displaystyle 0.4$$\displaystyle 0.5$

Figure 7: Simulation results for the normalized total conductivity for the
flower-shaped spot as a function of normalized petal’s parameter $n^{\prime}$
for different $\xi$, three independent least squares fit of function (9) are
also plotted; the three corresponding tangents at the origin are also plotted.

###  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. 8 and Fig. 9, 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”.

$\displaystyle n=10,\xi=0.1$$\displaystyle r_{0}$$\displaystyle
n=20,\xi=0.1$$\displaystyle r_{0}(1+\xi)$$\displaystyle
n=10,\xi=0.2$$\displaystyle n=20,\xi=0.2$

Figure 8:  Examples of star-shaped spots with $n\in\\{10,20\\}$ petals, and
half-petal length defined by $\xi\in\\{0.1,0.2\\}$.

$\displaystyle n=10,\xi=0.1$$\displaystyle r_{0}$$\displaystyle
n=20,\xi=0.1$$\displaystyle r_{0}(1+\xi)$$\displaystyle
n=10,\xi=0.2$$\displaystyle n=20,\xi=0.2$

Figure 9:  Examples of gear-shaped spots with $n\in\\{10,20\\}$ petals, and
half-petal lenght defined by $\xi\in\\{0.1,0.2\\}$.

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:

$\displaystyle P_{\text{\tiny star}}$
$\displaystyle=2\sqrt{2}nr_{0}\sqrt{1+\xi^{2}-(1-\xi^{2})\cos(\pi/n)},$ (13)
$\displaystyle A_{\text{\tiny star}}$
$\displaystyle=nr_{0}^{2}(1-\xi^{2})\sin(\pi/n),$ $\displaystyle
C_{\text{\tiny star}}$
$\displaystyle=\frac{\sqrt{(1-\xi^{2})\sin(\pi/n)}}{2\sqrt{2n}\sqrt{1+\xi^{2}-(1-\xi^{2})\cos(\pi/n)}}.$

$\displaystyle P_{\text{\tiny gear}}$ $\displaystyle=2\pi
r_{0}(1+2n^{\prime}/\pi),$ (14) $\displaystyle A_{\text{\tiny gear}}$
$\displaystyle=\pi r_{0}^{2}(1+\xi^{2}),$ $\displaystyle C_{\text{\tiny
gear}}$
$\displaystyle=\frac{\sqrt{(1+\xi^{2})}}{2\sqrt{\pi}(1+2n^{\prime}/\pi)}.$

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

$A_{\text{\tiny 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{\tiny star}}<A_{\text{\tiny flower}}<A_{\text{\tiny 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{\tiny flower}}\approx P_{\text{\tiny star}}<P_{\text{\tiny 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{\tiny gear}}<C_{\text{\tiny flower}}\approx C_{\text{\tiny 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. 10. 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 (8), and finally
presented in Fig. 11 complemented with previous results for the flower-shaped
spots. The least square fitted parameters for the same Eq. (9) 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^{\lim}_{\text{\tiny star}}\approx 0.903\quad<\quad Q^{\lim}_{\text{\tiny
flower}}=0.923\quad<\quad Q^{\lim}_{\text{\tiny gear}}=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.

Figure 10:  BEM result of the flux through a gear- and star-shaped spots of
$20$ petals and $\xi=0.1$.

$\displaystyle{0}$$\displaystyle{4}$$\displaystyle{8}$$\displaystyle{12}$$\displaystyle{16}$$\displaystyle{20}$$\displaystyle{24}$$\displaystyle
n^{\prime}$$\displaystyle{0.0}$$\displaystyle{0.2}$$\displaystyle{0.4}$$\displaystyle{0.6}$$\displaystyle{0.8}$$\displaystyle{1.0}$$\displaystyle
Q^{\prime}=(Q-Q_{\circ})/(Q_{up}-Q_{\circ})$Flower BEM,
$\displaystyle\xi=0.1$Star BEM, $\displaystyle\xi=0.1$Gear BEM,
$\displaystyle\xi=0.1$

Figure 11:  Results of the normalised flux for star-, gear- and flower-shaped spots with respect to the normalized ”petal” parameter $n^{\prime}=\xi n$. Least squares fit of Eq. (9) is also shown.  | Parameters | Coefficients
---|---|---
Simulation | $\xi$ | $n\in$ | $a$ | $b$ | $ab$
Gear | $0.1$ | $[4,256]$ | 0.978 | 0.894 | 0.875
Flower | $0.1$ | $[1,256]$ | 0.923 | 0.326 | 0.301
Star | $0.1$ | $[4,256]$ | 0.903 | 0.220 | 0.199
Table 2:  Fit parameters for Eq. (9) for the total flux of different multi-
petal shapes.

##  5 Conductivity of self-affine spots

###  5.1 Geometry of spots

Being inspired by shapes of contact clusters occuring 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.
13) 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
occuring 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:

$\displaystyle h(\theta)$
$\displaystyle=\sum_{k=k_{l}}^{k_{s}}\xi_{k}\cos{(k\theta+\theta^{0}_{k})},\quad
h(\theta)=h(\theta+2\pi),\quad\langle h\rangle=0$ (15a) $\displaystyle\xi_{k}$
$\displaystyle=\xi\left(\frac{k}{k_{l}}\right)^{-(0.5+H)}.$ (15b)

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^{0}_{k}$ which follow a uniform
distribution on $\theta^{0}_{k}\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)),$ (16)

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)),$ (17)

whose two first terms of Taylor expansion is equivalent to (16), but the
transformation (17) 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
(16),(17) are illustrated in Fig. 12 highlighting a situation, where for
transformation (16) the radius becomes negative. Therefore, for this study we
adopt the exponential transformation (17). Several examples of complex shapes
generated using the presented algorithm are shown in Fig. 13 for $\xi=0.1$ and
different values of $k_{l}$, $k_{s}$ and $H$ but with the same set of random
phases $\theta^{0}_{k}$. As the Hurst exponent increases, the spot naturally
becomes smoother.

$\displaystyle
h$(a)$\displaystyle\theta$$\displaystyle\pi$$\displaystyle-\pi$$\displaystyle
h=-1$$\displaystyle r_{0}$(b)$\displaystyle r_{0}\exp(h(\theta))$ Eq.
(17)$\displaystyle r_{0}(1+h(\theta))$ Eq. (16)

Figure 12: 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 (16) and exponential (17) transformations.

(a) $\displaystyle k_{l}=2,k_{s}=8$(b) $\displaystyle k_{l}=2,k_{s}=16$(c)
$\displaystyle k_{l}=2,k_{s}=32$(d) $\displaystyle k_{l}=4,k_{s}=8$(e)
$\displaystyle k_{l}=4,k_{s}=16$(f) $\displaystyle k_{l}=4,k_{s}=32$(g)
$\displaystyle H=0.25$(h) $\displaystyle H=0.5$(i) $\displaystyle H=0.75$

Figure 13: Examples of self-affine spots and their geometrical characteristics
for $\xi=0.1$.
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$.

###  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 (17). The obtained radius follows a log-normal
distribution

$P(r)=\frac{1}{r\sigma_{h}\sqrt{2\pi}}\exp{\left(-\frac{\ln^{2}(r)}{2\sigma_{h}^{2}}\right)}.$
(18)

Histograms shown in Fig. 14 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.

$\displaystyle{0.5}$$\displaystyle{1.0}$$\displaystyle{1.5}$$\displaystyle{0}$$\displaystyle{1}$$\displaystyle{2}$$\displaystyle{3}$(a)$\displaystyle{0.5}$$\displaystyle{1.0}$$\displaystyle{1.5}$(b)$\displaystyle
r_{0}\exp(h(\theta))$ Eq. (17)Log-Normal Eq. (18)NormalNormalized radius,
$\displaystyle r/r_{0}$Probability density, $\displaystyle P(r/r_{0})$

Figure 14:  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.

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=k_{l}}^{k_{s}}\left(\frac{k}{k_{l}}\right)^{-(1+2H)}$
(19)

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

$\displaystyle\langle r\rangle$
$\displaystyle=r_{0}\exp\left(\sigma_{h}^{2}/2\right),$ (20a)
$\displaystyle\sigma_{r}^{2}$
$\displaystyle=\left[\exp\left(\sigma_{h}^{2}\right)-1\right]\exp(\sigma_{h}^{2}),$
(20b)

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. (20b) 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
spectral moments $m_{0}$, $m_{2}$ and $m_{4}$ as follows:

$\displaystyle m_{0}$ $\displaystyle=\sigma^{2},$ (21a) $\displaystyle m_{2}$
$\displaystyle=\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,$ (21b) $\displaystyle m_{4}$
$\displaystyle=\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,$
(21c)

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 (20b)):

$\quad
m_{2}=\frac{\xi^{2}}{2}\sum_{k_{l}}^{k_{s}}k^{2}{\left(\frac{k}{k_{l}}\right)}^{-(1+2H)},\quad
m_{4}=\frac{\xi^{4}}{2r_{0}^{2}}\sum_{k_{l}}^{k_{s}}k^{4}{\left(\frac{k}{k_{l}}\right)}^{-(1+2H)}$
(22)

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

$m^{\text{\tiny f}}_{0}=\frac{r_{1}^{2}}{2}=\frac{\xi^{2}r^{2}_{0}}{2},\quad
m^{\text{\tiny
f}}_{2}=\frac{r_{1}^{2}n^{2}}{2r_{0}^{2}}=\frac{\xi^{2}n^{2}}{2}=\frac{n^{\prime
2}}{2},\quad m^{\text{\tiny
f}}_{4}=\frac{r_{1}^{2}n^{4}}{2r_{0}^{4}}=\frac{\xi^{2}n^{4}}{2r_{0}^{2}}=\frac{n^{\prime
2}n^{2}}{2r_{0}^{2}}$ (23)

Since the normalized conductivity for flower-shaped spots was shown to depend
exclusively on $n^{\prime}=\xi n$ (see Eq. (9)), 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 to take 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(|\mu(Q)-\langle
Q\rangle|\geq\frac{\sigma(Q)}{\sqrt{n\gamma}}\right)\leq\gamma,$ (24)

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. 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.

An example of flux distribution obtained by fast-BEM is displayed in Fig. 15
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. 15, the finest used mesh
reach $N_{e}=34\,340$ elements. As previously, to employ the Richardson
extrapolation, two meshes of different density were used to obtain accurate
results.

Figure 15:  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$.

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
$\mathcal{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 $\mathcal{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_{\circ}=4kU_{0}\langle r\rangle$ or (b) of the same area. The first option
was selected since the mean radius also enters in the set of parameters,
moreover $Q_{\circ}$ is equivalent to the definition of $Q_{\min}$ used for
the flower-shaped contact spot. The normalized total flux is thus defined by

$Q^{\prime}=\frac{Q}{Q_{\circ}}=\frac{Q}{4kU_{0}\langle r\rangle}$ (25)

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^{\prime}$ and consider two sets of dimensionless parameters:

  * •

Geometrical set of parameters $\mathcal{A}^{\prime}:=\\{\sqrt{m_{0}}/\langle
r\rangle,m_{2},\langle r\rangle\sqrt{m_{4}},H\\}$

  * •

Generative set of parameters $\mathcal{I}^{\prime}=\\{\xi,k_{l},k_{s},H\\}$

By analogy with the maximal total flux $Q_{\text{\tiny 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^{\text{\tiny f}}_{0}$, Eq.
(23), 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}^{\text{\tiny 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_{\text{\tiny up}}-Q_{\circ}=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^{\prime\prime}=\frac{Q-Q_{\circ}}{Q_{\text{\tiny
up}}-Q_{\circ}}=\frac{Q}{4kU_{0}\sqrt{2m_{0}}}-\frac{\langle
r\rangle}{\sqrt{2m_{0}}}$ (26)

Figure 16:  Normalized flux (25) of self-affine spots: for the following
generative parameters $\mathcal{I}^{\prime}$: $\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\%$).

The total flux normalized according to Eq. (25) is presented in Fig.16 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}^{\prime}$ offers a better representation than the use of the
generative set of parameters $\mathcal{I}^{\prime}$. 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}}$.

Figure 17:  Renormalized flux, Eq. (26) 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. 18.

The results of the renormalized flux Eq.(26) with respect to $\sqrt{m_{2}}$
and $\langle r\rangle\sqrt{m_{4}}$ are presented in Fig.17. 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.

To provide a visual geometrical interpretation in the flux variation (see Fig.
18) 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. 17. 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.

Figure 18:  Examples of self-affine spots with corresponding generative characteristics; the corresponding data points are highlighted in Fig. 17, the resulting flux and geometrical characteristics can be found in Table 3. |  | Parameters $\mathcal{A}^{\prime}$ | Parameters $\mathcal{I}^{\prime}$
---|---|---|---
Spot # | $Q^{\prime\prime}$ | $\sqrt{m_{0}}/\langle r\rangle$ | $\sqrt{m_{2}}$ | $\sqrt{m_{4}}\langle r\rangle$ | $\xi$ | $k_{l}$ | $k_{s}$ | $H$
$S_{1}$ | $0.151$ | $0.0664$ | $0.658$ | $14$ | $0.05$ | $2$ | $32$ | $0.25$
$S_{2}$ | $0.172$ | $0.0664$ | $0.666$ | $12.7$ | $0.05$ | $4$ | $32$ | $0.6$
$S_{3}$ | $0.173$ | $0.0657$ | $0.958$ | $39.8$ | $0.05$ | $2$ | $64$ | $0.3$
$S_{4}$ | $0.184$ | $0.0645$ | $0.972$ | $63.8$ | $0.05$ | $4$ | $128$ | $0.7$
$S_{5}$ | $0.317$ | $0.0800$ | $2.06$ | $163$ | $0.05$ | $4$ | $128$ | $0.4$
$S_{6}$ | $0.380$ | $0.0852$ | $1.97$ | $130$ | $0.05$ | $8$ | $128$ | $0.75$
$S_{7}$ | $0.644$ | $0.126$ | $5.12$ | $446$ | $0.05$ | $8$ | $128$ | $0.25$
$S_{8}$ | $0.749$ | $0.182$ | $8.77$ | $1493$ | $0.05$ | $8$ | $256$ | $0.25$
Table 3: Parameters for spots shown in Fig. 18 and the resulting total flux.

###  5.4 Results with renormalized standard deviation

As presented in Figs. 16 and 17, 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. (19)):

$h(\theta)=\sqrt{\frac{2m_{0,h}}{s}}\sum_{k=k_{l}}^{k_{s}}\xi_{k}\cos{(k\theta+\theta^{0}_{k})},\quad
s=\sum_{k=k_{l}}^{k_{s}}\xi_{k}^{2}$ (27)

with $\xi_{k}$ defined by Eq. (15b). The exponential transformation from
$h(\theta)$ to radius $r(\theta)$ (17) remains intact. Then, for the
normalized generative parameters we have:
$\mathcal{I}^{\prime\prime}=\\{\sqrt{m_{0,h}},k_{l},k_{s},H\\}$ and for the
geometrical ones we still have $\mathcal{A}^{\prime}$. 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 (20b) and its Taylor expansion.

The results of this set of simulations for the renormalised flux (26) are
presented in Fig. 19 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 in the flux normalization (26).
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
magnification $\zeta$.

To provide a geometrical meaning to these results, a 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.
20 and are highlighted in Fig. 19 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}$.

Figure 19:  Results for the renormalized total flux Eq. (26) 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$.

Figure 20:  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.  |  | Parameters $\mathcal{A}^{\prime}$ | Parameters $\mathcal{I}^{\prime\prime}$
---|---|---|---
Spot # | $Q^{\prime\prime}$ | $\sqrt{m_{0}}/\langle r\rangle$ | $\sqrt{m_{2}}$ | $\sqrt{m_{4}}\langle r\rangle$ | $\sqrt{m_{0,h}}$ | $k_{l}$ | $k_{s}$ | $H$
$S_{1}$ | $0.189$ | $0.100$ | $0.507$ | $5.10$ | $0.100$ | $2$ | $16$ | $0.5$
$S_{2}$ | $0.213$ | $0.143$ | $0.515$ | $2.84$ | $0.141$ | $2$ | $8$ | $0.5$
$S_{3}$ | $0.510$ | $0.0995$ | $1.56$ | $33.6$ | $0.100$ | $8$ | $32$ | $0.5$
$S_{4}$ | $0.476$ | $0.0705$ | $1.56$ | $61.9$ | $0.005$ | $8$ | $64$ | $0.5$
$S_{5}$ | $0.665$ | $0.0997$ | $2.62$ | $114$ | $0.100$ | $8$ | $64$ | $0.25$
$S_{6}$ | $0.666$ | $0.141$ | $2.61$ | $94.1$ | $0.141$ | $8$ | $64$ | $0.75$
Table 4: Parameters for spots shown in Fig. 20 and the resulting total flux.
The spots are ordered according to increasing $\sqrt{m_{0}}$.

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 for multi-petal spots and we recall that
$\sqrt{m_{2}}$ is analogous parameter to $n^{\prime}$ 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. 7 and 11). For the extra generative parameter
$\zeta\in\mathcal{I}^{\prime\prime}$ 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})},$ (28)

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}^{\prime}$,
$\tilde{Q}^{\prime\prime}$ and parameters $\tilde{k}_{l}$,
$\tilde{\sqrt{m_{0}}/\langle r\rangle}$, $\tilde{\sqrt{m_{2}}}$ and
$\tilde{\sqrt{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 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 $\mathcal{I}$ are strongly
linked to the method of spot generation, in constructing our model we will
focus exclusively on geometrical parameters $\mathcal{A}$ which could be
measured for an arbitrary shape.

| $\tilde{Q}^{\prime}$ | $\tilde{Q}^{\prime\prime}$ | $\tilde{\xi}$ | $\tilde{k}_{l}$ | $\tilde{\zeta}$ | $\tilde{H}$ | $\tilde{\sqrt{m_{0}}/\langle r\rangle}$ | $\tilde{\sqrt{m_{2}}}$ | $\tilde{\sqrt{m_{4}}\langle r\rangle}$ | $\tilde{\alpha}$
---|---|---|---|---|---|---|---|---|---|---
$\tilde{Q}^{\prime}$ | 1 | 0.9 | 0.1 | 0.7 | 0.3 | -0.2 | 0.7 | 0.9 | 0.7 | 0.06
$\tilde{Q}^{\prime\prime}$ | - | 1 | -0.04 | 0.8 | 0.4 | -0.2 | 0.7 | 0.9 | 0.7 | 0.1
$\tilde{\xi}$ | - | - | 1 | -0.2 | -0.2 | 0.1 | 0.7 | -0.06 | -0.07 | -0.08
$\tilde{k}_{l}$ | - | - | - | 1 | 0.06 | 0.0 | 0.4 | 0.7 | 0.4 | -0.03
$\tilde{\zeta}$ | - | - | - | - | 1 | 0.0 | 0.06 | 0.5 | 0.5 | 0.9
$\tilde{H}$ | - | - | - | - | - | 1 | -0.2 | -0.2 | 0.3 | 0.2
$\tilde{\sqrt{m_{0}}/\langle r\rangle}$ | - | - | - | - | - | - | 1 | 0.6 | 0.4 | -0.05
$\tilde{\sqrt{m_{2}}}$ | - | - | - | - | - | - | - | 1 | 0.9 | 0.2
$\tilde{\sqrt{m_{4}}\langle r\rangle}$ | - | - | - | - | - | - | - | - | 1 | 0.3
$\tilde{\alpha}$ | - | - | - | - | - | - | - | - | - | 1
Table 5: Covariance matrix of normalized flux and all normalized parameters
according to Eq. (28).

####  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. (23) for the flower-
shaped spot, the Nayak parameter is simply $1$ so, consistenly it does not
enter the phenomenological equation for the conductivity of such simple forms
Eq. (9). 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 covariance analysis. So, the ultimate set of geometrical
dimensionless parameters is chosen to be:

$\mathcal{A}^{f}=\left\\{\sqrt{m_{0}}/\langle
r\rangle,\sqrt{m_{2}},H,\alpha\right\\}\equiv\left\\{\sigma/\langle
r\rangle,\sqrt{\langle(\nabla r)^{2}\rangle},H,\alpha\right\\}$ (29)

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^{\prime\prime}=a\left[1-\frac{1}{b\sqrt{m_{2}}+1}\right]\left(1+cH\right)\left\\{1+\frac{d}{e\alpha^{f}+1}\right\\}$
(30)

with the core term in square brackets which is equivalent to the
phenomenological law obtained for multi-petal spots, see Eq. (9). In addition,
the effects of $H$ and $\alpha$ enter in 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 decreaseing
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. ensure always posivite normalized flux
$Q^{\prime\prime}$ 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.(30) and (26), the final equation for the flux is obtained as:

$Q=Q_{\circ}\left(1+a\frac{\sqrt{2m_{0}}}{\langle
r\rangle}\left[1-\frac{1}{b\sqrt{m_{2}}+1}\right]\left(1+cH\right)\left\\{1+\frac{d}{e\alpha^{f}+1}\right\\}\right)$
(31)

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. 21 separately for
two sets of simulation data: in Fig. 21(a) for the set of contact spots
parametrized by Eq.(15a), and in Fig. 21(b) for those defined by Eq. (27). A
relative error could be defined as:

$E=\frac{1}{N}\sum\limits_{i=1}^{N}\frac{|Q_{i}-Q_{i}^{\text{\tiny
fit}}|}{Q_{i}},$ (32)

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$.

Parameter | $a$ | $b$ | $c$ | $d$ | $e$ | $f$
---|---|---|---|---|---|---
Value | $0.968$ | $0.255$ | $0.0867$ | $4.38$ | $5.49$ | $1.50$
Table 6:  Parameters of the phenomenological model (30),(31) optimized through
least square fit and resulting in relative error $E=4.3$ % and Pearson’s
correlation factor $\rho=0.9976$.

Figure 21:  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. 19.

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 higest 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(r-\langle
r\rangle)^{2}\rangle}=\sqrt{m_{0}}$ and $\sqrt{\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 (31) 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\limits_{\zeta\to\infty}(Q)=Q_{\circ}\left(1+a\frac{\sqrt{2m_{0}}}{\langle
r\rangle}\left(1+cH\right)\right)=4kU_{0}\left(\langle
r\rangle+a\sqrt{2m_{0}}\left(1+cH\right)\right).$ (33)

As a first order approximation, one could use the following value
$4kU_{0}\left(\langle r\rangle+\sqrt{2m_{0}}\right)$ 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 conjucture 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 quantiative 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 eventually convective heat exchanges
or because of tunneling effects. The physical conductivity should probably hit
the upper limit $Q_{\text{\tiny 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 zero-th 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. (30) and (31)). 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. 22. Qualitatively, these results align with the
observed conductivity saturation at the fractal limit, a natural outcome for
the physical problem of conduction.

Figure 22:  Normalized flux at Koch snowflake spots at several first fractal
iterations and the evolution of the total flux.

## 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 issues. In this work, a Large Language Model
assisted the authors with polishing the text and, to a little extent, with the
code improvement.

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

Geometrical characteristics (perimeter, area and compactness) are summarized
in Table 7 and Fig. 23 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$.

Characteristics | Flower | Star | Gear
---|---|---|---
Area $A$ | $\pi r_{0}^{2}\left(1+\frac{\xi^{2}}{2}\right)$ | $nr_{0}^{2}(1-\xi^{2})\sin(\pi/n)$ | $\pi r_{0}^{2}(1+\xi^{2})$
Perim. $P$ | $r_{0}E(in^{\prime})$ | $2\sqrt{2}nr_{0}\sqrt{1+\xi^{2}-(1-\xi^{2})\cos(\pi/n)}$ | $2\pi r_{0}(1+2n^{\prime}/\pi)$
Comp. $C=\frac{\sqrt{A}}{P}$ | $\displaystyle\sqrt{\pi}\frac{\sqrt{1+\xi^{2}/2}}{4E(in^{\prime})}$ | $\displaystyle\frac{\sqrt{(1-\xi^{2})\sin(\pi/n)}}{2\sqrt{2n}\sqrt{1+\xi^{2}-(1-\xi^{2})\cos(\pi/n)}}$ | $\displaystyle\frac{\sqrt{(1+\xi^{2})}}{2\sqrt{\pi}(1+2n^{\prime}/\pi)}$
Table 7: Geometrical characteristics (area, perimeter and compactness) for
multi-petal shapes: flower, star and gear;
$E(in^{\prime})=E(in\xi)=\int\limits_{0}^{\pi/2}\sqrt{1+\left(n\xi\sin(\theta)\right)^{2}}d\theta$
is the complete elliptic integral of the second kind and $i$ is the imaginary
unit. Figure 23:  Perimeter, area and compactness for the different shapes of
contact spot

##  Appendix B Geometrical characteristics of self-affine spots

In Fig. 24, the analytical form for the standard deviation of self-affine
spots $\sigma_{s}$ Eq.(20) is compared with the one evaluated over $1\,000$
generated spots for each combination of generative parameters: $\xi=0.05$,
$k_{l}\in\\{2,4,8,16,32,64,126\\}$, $\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 forth $m_{4}$ moments have been computed for the
same set of generative parameters over the same $1\,000$ 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. (21).
The discretization consists of splitting the contour in
$N=\max\\{10\,000,100k_{s}\\}$ straight segments with $d\theta=2\pi/N$ and
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}},$
(34)

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. (22):

$m_{p}^{S}=\frac{(r_{0}\xi)^{2}}{2}\sum_{k_{l}}^{k_{s}}k^{p}{\left(\frac{k}{k_{l}}\right)}^{-(1+2H)}.$
(35)

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

$m_{p}^{C}=\frac{1}{2}\int_{k_{l}}^{\zeta k_{l}}\left(r(k)\right)^{2}k^{p}dk,$
(36)

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}\left(\zeta^{-2H}-1\right),\quad
m_{2}^{C}=-\frac{(r_{0}\xi)^{2}k_{l}^{3}}{2(2-2H)}\left(\zeta^{2-2H}-1\right),\quad
m_{4}^{C}=-\frac{(r_{0}\xi)^{2}k_{l}^{5}}{2(4-2H)}\left(\zeta^{4-2H}-1\right).$
(37)

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. (22) as detailed by
Nayak [Nayak, 1971]. These analytical values are compared with the numerically
evaluated ones in Fig.25 and Fig.26. 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 $1\,000$
spots are compared in Fig 27. 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}$.

Figure 24: 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$.

Figure 25: 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$.

Figure 26: 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$. Figure 27:
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.

##  Appendix C  Physical consistency of the phenomenological model

The phenomenological model of flux Eqs. (30),(31) 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$ (38)

The terms $\partial Q/\partial m_{2}$ and $\partial Q/\partial\alpha$ could be
redily derived from Eq. (31).

$\displaystyle\frac{\partial Q}{\partial m_{2}}$
$\displaystyle=a\frac{b}{2\sqrt{m_{2}}(b\sqrt{m_{2}}+1)^{2}}(1+cH)\left(1+\frac{d}{e\alpha^{f}+1}\right)$
(39a) $\displaystyle\frac{\partial Q}{\partial\alpha}$
$\displaystyle=a\left\\{1-\frac{b}{b\sqrt{m_{2}}+1}\right\\}(1+cH)\frac{-def\alpha^{f-1}}{(e\alpha^{f}+1)^{2}}$
(39b)

The derivatives of $m_{2}$ and $\alpha$ with respect to $\zeta$ could be found
from Eq.(37), resulting in the following asymptotic forms:

$\displaystyle\frac{\partial Q}{\partial m_{2}}\frac{\partial
m_{2}}{\partial\zeta}$
$\displaystyle\sim\frac{1}{\sqrt{m_{2}}{\left(\sqrt{m_{2}}+1\right)}^{2}}\frac{\partial
m_{2}}{\partial\zeta}\sim\zeta^{-2+H}$ (40a) $\displaystyle\frac{\partial
Q}{\partial\alpha}\frac{\partial\alpha}{\partial\zeta}$
$\displaystyle\sim\frac{\alpha^{f-1}}{{\left(\alpha^{f}+1\right)}^{2}}\frac{\partial\alpha}{\partial\zeta}\sim\zeta^{-2fH-1}$
(40b)

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

$f\geq\frac{1-H}{2H}$ (41)

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\geq 1.5$. The results for derivatives using the continuous expressions for
$m_{2}$ and $\alpha$ are presented Fig. 28 for $H=0.25$, $k_{l}=8$, and
$\xi=0.05$, and the fitting parameters shown in Table 6. The two derivative
terms are distingished: one positive, as given by Eq. (39a), and the other
negative, as given by Eq. (39a). The full derivative remains positive,
however, thus keeping the required assumptions true, even for value of $\zeta$
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
$\zeta$.

Figure 28:  Derivatives of the flux phenomenological model in function of the
magnification $\zeta$ for $k_{l}=8$, $H=0.25$, and $\xi=0.05$ for parameters
from Table 6.

## Conflicts of interest

The authors declares no competing financial interest.

## Supplementary data

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.

## References

  * [Afferrante et al., 2012] Afferrante, L., Carbone, G., and Demelio, G. (2012).  Interacting and coalescing hertzian asperities: A new multiasperity 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 and Yastrebov, 2023a] Beguin, P. and Yastrebov, V. (2023a).  Fast-bem code for poisson’s equation on a half-space.  https://github.com/vyastreb/HBEM, version 1.0.
  * [Beguin and Yastrebov, 2023b] Beguin, P. and Yastrebov, V. A. (2023b).  Supplementary material to paper.  Zenodo archive https://doi.org/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.  Meccanica, 34(4):301–302.
  * [Bowden et al., 2001] Bowden, F. P., Bowden, F. P., and Tabor, D. (2001).  The friction and lubrication of solids, volume 1.  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. (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., Mikic, B., and Yovanovich, M. (1969).  Thermal contact conductance.  International Journal of Heat and Mass Transfer, 12(3):279–300.
  * [Cuthrell and Tipping, 1973] Cuthrell, R. and Tipping, D. (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., 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. and Robson, J. (1973).  The description of road surface roughness.  Journal of Sound and Vibration, 31(2):175–183.
  * [Fabrikant, 1993] Fabrikant, V. (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. (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 $\mathcal{H}$-matrices for BEM.  Computing, 74(3):205–223.
  * [Grasedyck and Hackbusch, 2003] Grasedyck, L. and Hackbusch, W. (2003).  Construction and arithmetics of $\mathcal{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. 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, volume 49.  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 tribological 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. (2009).  A high-precision apparatus for the characterization of thermal interface materials.  Review of Scientific Instruments, 80(9):095111.
  * [Koch, 1904] Koch, H. (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. and Fletcher, L. (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. (1990).  Fractal characterization and simulation of rough surfaces.  Wear, 136(2):313–327.
  * [Manners and Greenwood, 2006] Manners, W. and Greenwood, J. (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., 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. (1971).  Random process model of rough surfaces.  Trans. ASME, J. Lub. Tech. 93F, 3:398.
  * [Pei et al., 2005] Pei, L., Hyun, S., Molinari, J., 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 America, 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.
  * [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.