# Topological spin Hall effect in antiferromagnets

Amir N. Zarezad  Józef Barnaś  A. Dyrdał<EMAIL_ADDRESS>Faculty of
Physics, ISQI, Adam Mickiewicz University in Poznań, ul. Uniwersytetu
Poznańskiego 2, 61-614 Poznań, Poland
Alireza Qaiumzadeh  Center for Quantum Spintronics, Department of Physics,
Norwegian University of Science and Technology, NO-7491 Trondheim, Norway

###### Abstract

We investigate topological Hall effects in a metallic antiferromagnetic (AFM)
thin film and/or at the interface of an AFM insulator–normal metal bilayer
with a single skyrmion in the diffusive regime. To determine the spin and
charge Hall currents, we employed a Boltzmann kinetic equation with both spin-
dependent and spin-flip scatterings. The interaction between conduction
electrons and static skyrmions is included in the Boltzmann equation via the
corresponding emergent magnetic field arising from the skyrmion texture. We
compute intrinsic and extrinsic contributions to the topological spin Hall
effect and spin accumulation, induced by an AFM skyrmion. We show that
although the spin Hall current vanishes rapidly outside the skyrmion, the spin
accumulation can be finite at the edges far from the skyrmion, provided the
spin diffusion length is longer than the skyrmion radius. In addition, We show
that in the presence of a spin-dependent relaxation time, the topological
charge Hall effect is finite and we determine the corresponding Hall voltage.
Our results may help to explore antiferromagnetic skyrmions by electrical
means in real materials.

##  I Introduction

Interplay between charge currents and magnetic textures is important from both
fundamental and applied points of view [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
13, 14]. Interactions between current and ferromagnetic skyrmions may lead to
emerging phenomena, including topological Hall and skyrmion Hall effects [15,
16]. The spin-polarized charge current flowing in a ferromagnetic (FM) layer
containing a topological skyrmion drags it along the electric field and
simultaneously deflects it toward one of the film edges, that is, along the
direction perpendicular to the electric field under the influence of
gyrotropic forces. The latter phenomenon is called the skyrmion Hall effect
[17, 18, 19, 20, 21, 22]. This phenomenon was studied theoretically and
observed experimentally; see, e.g. Refs. [20, 21, 23, 24, 25, 26, 27] for an
overview. The motion of rigid topological skyrmions in both longitudinal and
transverse directions in the present of charge currents in FM metals are
conveniently described by Thiele’s equation [19, 28, 29]. In turn, FM
skyrmions affect the flow of spin-polarized charge currents in FM metal by
deflecting their trajectories in the direction perpendicular to the external
electric field [30, 31, 32, 33, 34, 35, 36, 37]. This is a phenomenon similar
to the anomalous Hall effect in FM metals with uniform magnetization direction
[38, 39, 40]. The origin of this skyrmion-induced topological Hall effect is
the emergence of a real-space Berry curvature, induced by skyrmion textures
[31, 41], while in the anomalous Hall effect, the Berry curvature in the
momentum space arises from the spin-orbit couplings [42, 43, 44]. The FM
skyrmion Hall effect and the topological Hall effect of spin-polarized
currents in FM metals are reciprocal phenomena with a similar origin.

On the other hand, AFM skyrmions, unlike their FM counterparts, do not show
any skyrmion Hall effect [45, 46, 47, 48, 49, 50, 51, 52]. This deflection-
free motion of AFM skyrmions is one of the advantages of AFM systems versus
their FM counterparts in practical spintronic devices. However, topological
skyrmion textures may still create a real-space Berry curvature [53] and hence
some type of topological Hall effects are expected to be present in AFM
systems in the presence of skyrmions. Charge currents in compensated AFM
metals are not spin polarized, and thus the two equally populated opposite
spins can be deflected in opposite directions in the presence of the AFM
skyrmion-induced Berry curvature. Therefore, one can expect the topological
_spin_ Hall effect instead of the topological charge Hall effect. Recently,
this effect has been investigated theoretically in a few publications [54, 55,
56, 57], but has not been measured experimentally yet. The absence of a
topological charge Hall effect was recently confirmed in a chiral AFM system
[58].

In Ref. [54], the authors employed an SU(2) semiclassical framework in
combination with the ab-initio description of the electronic structure and
demonstrated the emergence of a sizeable transverse spin current and a
vanishing transverse charge current in a synthetic AFM skyrmion lattice. In
Ref. [55], the authors showed a finite topological spin Hall effect and the
absence of topological charge Hall effect in a compensated AFM skyrmion
crystal on a honeycomb lattice. Akosa et al [56] computed the topological spin
and charge Hall effects in a finite size AFM square lattice using a tight-
binding model in terms of Landauer–Büttiker formalism, implemented in the
Kwant code, in the presence of electrostatic impurities. They found a zero
topological charge Hall effect and a nonzero topological spin Hall effect in
such systems. Finally, in Ref. [57], using Landauer–Büttiker formula, the
authors showed that a vector chirality, formed by the AFM Néel vector, gives
rise to a finite topological spin Hall effect in bulk system of AFM half-
skyrmions or merons but not in skyrmions.

In the present article, we investigate the topological charge and spin Hall
effects in a diffusive AFM layer with a skyrmion texture, by means of the
quasiclassical Boltzmann kinetic equation. We consider both spin-dependent and
spin-flip scatterings in our formalism, and we find analytical expressions for
transverse spin and charge currents as well as spin accumulation in the
presence of an AFM skyrmion. We show that spin-dependent scatterings generate
a nonzero topological charge Hall effect. In addition, we compute both
intrinsic (disorder independent, which is arising from the real-space Berry
curvature) and extrinsic (disorder dependent, which is arising from the spin-
flip scatterings) contributions to the topological spin Hall effect and spin
accumulation.

##  II Model

We consider a two-sublattice compensated AFM square lattice that hosts a
single static skyrmion. The system consists of an electronic subsystem and a
AFM spin subsystem that interact to each other via an sd exchange interaction.
The effective Hamiltonian of the system can be written in the following form
[56, 59, 60, 61]:

$\hat{H}=-t\gamma_{\mathbf{k}}\big{(}\hat{\tau}_{x}\otimes\hat{s}_{0}\big{)}-J\big{(}\hat{\tau}_{z}\otimes\bf{n}\big{)}\cdot\hat{\bf{s}}.$
(1)

Here, ${t>0}$ is the hopping parameter,
${\gamma_{\mathbf{k}}=z^{-1}\sum_{i=1}^{z}\exp{(-i\bm{k}\cdot\bm{\delta}_{i})}}$
is the lattice structure factor with $\bm{\delta}_{i}$ denoting the nearest-
neighbor vectors and $z$ being the coordination number. For a square lattice
($z$=4), we have
$\gamma_{\mathbf{k}}=\big{(}\cos(k_{x}a)+\cos(k_{y}a)\big{)}/2$, where $a$ is
the lattice constant. Furthermore, $J$ is the sd exchange interaction that
parameterizes the coupling of mobile electrons and localized magnetic moments,
$\hat{\bm{\tau}}$ and $\hat{\bm{s}}$ are the vectors of Pauli matrices
representing the AFM sublattice and spin degrees of freedom, respectively,
$\hat{s}_{0}$ is the identity matrix in the spin subspace, and $\mathbf{n}$ is
the unit vector of the staggered Néel order parameter field. In general, the
unit vector of the staggered order parameter field can be decomposed into two
terms, $\mathbf{n}=\mathbf{n}_{0}+\mathbf{n}_{r}$, where $\mathbf{n}_{0}$
represents the background homogeneous field while the second term
$\mathbf{n}_{r}$ represents the noncollinear magnetic texture [62].

The corresponding eigenvalues and eigenstates of the AFM Hamiltonian in the
absence of any texture, $\bf n_{r}=0$, are given by [56, 59];

$\varepsilon_{\eta}(\mathbf{k})=\eta\sqrt{t^{2}\gamma_{\mathbf{k}}^{2}+J^{2}},$
(2)
$\ket{\Psi_{\eta}^{s}}=\dfrac{1}{\sqrt{2}}\Big{(}\sqrt{1+s\eta\mathrm{P}_{\mathbf{k}}}\ket{A}+\eta\sqrt{1-s\eta\mathrm{P}_{\mathbf{k}}}\ket{B}\Big{)}\otimes\ket{\bf\sigma},$
(3)

where, $\eta=+1$ ($-1$) corresponds to the conduction (valence) band, and
$s=+1$ ($-1$) corresponds to the spin up (down) state, $\ket{A(B)}$ refers to
the AFM sublattice A (B) projection, $\ket{\sigma}=\ket{\uparrow(\downarrow)}$
denotes up (down) spin projection, and finally we defined
$\mathrm{P}_{\mathbf{k}}=J/\sqrt{t^{2}\gamma_{\mathbf{k}}^{2}+J^{2}}$.

We consider a compensated AFM system that preserved combined time and
inversion symmetry ($PT$ symmetry) and thus the electronic band dispersion,
Eq. (2), is spin-degenerate. To find analytical results, we assume the Fermi
level near the maximum of the conduction (minimum of the valence) band, and
consider the electronic dispersion around the $\Gamma$ point. Accordingly, the
structure factor is approximated as
$\gamma_{\mathbf{k}}\approx(1-a^{2}k^{2}/4)$, while the energy dispersion
around the $\Gamma$ point becomes
$\varepsilon_{\eta}\approx\eta\big{(}\sqrt{J^{2}+t^{2}}-\hbar^{2}k^{2}/(2m_{\rm
eff})\big{)}$. The first term here is the energy at the $\Gamma$ point and the
second term is the effective kinetic energy of electrons with an effective
mass $m_{\rm eff}=2\hbar^{2}\sqrt{J^{2}+t^{2}}/a^{2}t^{2}$.

The profile of an AFM skyrmion can be modelled in the spherical coordinates
as,

$\mathbf{n}_{r}=\big{(}\cos{\Phi}\sin{\Theta},\sin{\Phi}\sin{\Theta},\cos{\Theta}\big{)},$
(4)

where the polar $\Theta(r)$ and azimuthal $\Phi(r)$ angles are defined as [36,
56]

$\displaystyle\Theta$ $\displaystyle=2\pi-4\,\arctan\Bigl{(}\exp(4r/r_{\rm
sk})\Bigr{)},$ (5a) $\displaystyle\Phi$
$\displaystyle=q\rm{Arg}(x+iy)+c\frac{\pi}{2}.$ (5b)

Here, $r=\sqrt{x^{2}+y^{2}}$ is the distance from the skyrmion center, $r_{\rm
sk}$ denotes the radius of the skyrmion core, $p=\pm 1$ stands for the
skyrmion polarity, $q=\pm 1$ denotes skyrmion vorticity, and $c=\pm 1$ defines
the chirality of the skyrmion.

For a spin texture slowly varying in space, the exchange interaction term –
the second term in the Hamiltonian (1) – can be diagonalized by a unitary
gauge transformation. The result is a uniform spin background in the presence
of an emerging SU(2) gauge field that interacts with itinerant electrons [9,
53, 41, 63, 64, 65]. The emergent magnetic field induced by an AFM texture
depends on the spin, band, and sublattice indices [56, 66],

$B^{\alpha,s}_{\mathrm{em},\eta}(\bm{r})=-s\big{(}1+s\eta\alpha\mathrm{P}_{\mathbf{k}}\big{)}\dfrac{\hbar}{2e}N_{x,y}(\bm{r})\hat{z},\\\
$ (6)

where $\alpha=+(-)$ refers to the $A(B)$ sublattice, $\hbar=h/(2\pi)$ is the
reduced Planck constant, $e$ is the electron charge, and $N_{x,y}(\bm{r})={\bf
n}\cdot(\partial_{x}{\bf n}\times\partial_{y}{\bf n})$ is the topological
charge density.

Without loss of generality, we assume the Fermi level is located in the
conduction band and in the rest of this article, we set $\eta=+1$ and drop
this index. We consider a narrow AFM stripe of width $2w\geq 2r_{\rm sk}$ and
length $2L$, that includes a single skyrmion at its center.

##  III Boltzmann kinetic equation

To describe spin and charge transports in the diffusive regime in the presence
of emerging magnetic field induced by a static AFM skyrmion, we employ a
semiclassical transport theory based on the Boltzmann kinetic equation [67,
68, 69],

$\displaystyle\mathbf{v_{k}}\cdot\frac{\partial
f_{s}}{\partial\mathbf{r}}-\frac{e}{\hbar}\Bigl{(}\mathbf{E}+\mathbf{v_{k}}\times\mathbf{B}_{\mathrm{em}}^{s}\Bigr{)}\cdot\frac{\partial
f_{s}}{\partial\mathbf{k}}$ (7) $\displaystyle=-\frac{f_{s}-\langle
f_{s}\rangle}{\tau_{s}}-\frac{\langle f_{s}\rangle-\langle
f_{-s}\rangle}{\tau_{\rm sf}},$

where $f_{s}=f_{s}(\mathbf{r},\mathbf{k})$ is the nonequilibrium distribution
function for electrons with spin $s=+/-$ (or equivalently
$s=\uparrow/\downarrow$), $\bf v_{k}$ is the electron velocity,
$\mathbf{E}=E_{x}\hat{x}$ is the external electric field applied along the
stripe, $\langle f_{s}\rangle=\int d^{2}\Omega_{\mathbf{k}}f_{s}/\int
d^{2}\Omega_{\mathbf{k}}$ is the angular average over the momentum space, and
$\Omega_{\mathbf{k}}$ is the solid angle in the momentum space. The first term
on the right hand side describes the spin-conserving relaxation processes with
$\tau_{s=\uparrow(\downarrow)}$ being the corresponding spin-dependent
scattering time. The second term describes spin mixing relaxation processes
with $\tau_{\rm sf}$ denoting the corresponding spin-flip relaxation time. The
total spin-dependent emerging magnetic field is the sum of the two sublattices
contributions,

$\mathbf{B}_{\mathrm{em}}^{s}(\bm{r})=\sum_{\alpha=A,B}{\bf
B}^{\alpha,s}_{\mathrm{em}}(\bm{r})=sB_{\rm em}(\bm{r})\hat{z},$ (8)

where for a skyrmion profile, defined in Eqs. (4) and (5), the emerging
magnetic field amplitude reads [18],

$\displaystyle\dfrac{B_{\rm
em}}{B_{0}}=-\frac{8}{r/r_{sk}}\,\sin\left(4\arctan\left(\exp\left(\frac{4r}{r_{\rm
sk}}\right)\right)\right)\frac{\exp(\frac{4r}{r_{\rm
sk}})}{1+\exp(\frac{8r}{r_{\rm sk}})},$ (9)

with $B_{0}=(h/e)(\pi r_{\rm sk}^{2})^{-1}$. The emergent magnetic field is
normal to the AFM layer and has an opposite sign for up and down itinerant
electron spins. Therefore, this magnetic field effectively deflects the
electron trajectory of spin-up and spin-down electrons in opposite directions.
If there is no asymmetry between the spin up and spin down electron subbands,
this leads to a vanishing topological charge Hall effect, while the spin Hall
effect is nonzero. However, in the presence of an asymmetry between spin
subbands, e.g., due to different relaxation times or breaking the time-
reversal symmetry of electronic bands, both charge and spin Hall effects may
occur, as we will show later.

##  IV Analytical solutions of Boltzmann equation

To solve the Boltzmann kinetic equation, Eq. (7), we decompose the
nonequilibrium distribution function into an equilibrium component $f_{s}^{0}$
and small nonequilibrium perturbations [37, 70, 71],

$f_{s}=f_{s}^{0}-\frac{\partial
f_{s}^{0}}{\partial\varepsilon}\Bigl{(}-e\mu_{s}(\mathbf{r})+g_{s}(\mathbf{r},\mathbf{k})\Bigr{)},$
(10)

where $-e\mu_{s}(\mathbf{r})$ and $g_{s}(\mathbf{r},\mathbf{k})$ are the
isotropic (zeroth velocity moment) and anisotropic (first velocity moment)
parts, respectively. The first part is related to spin accumulation, while the
second part represents a shift of the electron sphere in the momentum space.
Both $\mu_{s}(\mathbf{r})$ and $g_{s}(\mathbf{r},\mathbf{k})$ should be
determined from the Boltzmann equation. Inserting Eq. (10) into Eg. (7), and
separating the odd and even velocity moments of the distribution function, up
to linear order in the emerging magnetic field, we find [37, 70, 71],

$\displaystyle g_{s}({\mathbf{k}},{\bf r})=-e\tau_{s}{\bf
v}_{\mathbf{k}}\cdot\left({\bf{E}}-\nabla_{\bf
r}\mu_{s}({\mathbf{r}})-\frac{e\tau_{s}}{m_{\rm eff}}{\bf
E}\times\mathbf{B}_{\rm em}^{s}(\bm{r})\right),$ (11)
$\displaystyle\nabla^{2}\delta\mu(\mathbf{r})-\dfrac{\delta\mu(\mathbf{r})}{\lambda_{\rm
sd}^{2}}=\dfrac{e\tau}{m_{\mathrm{eff}}}(\hat{z}\times{\bf
E})\cdot\nabla\mathbf{\bm{B}}_{\rm em}(\bm{r}).$ (12)

Here, $\delta\mu=(\mu_{\uparrow}-\mu_{\downarrow})/2$ is the net spin
accumulation, $\tau=(\tau_{\uparrow}+\tau_{\downarrow})/2$ is the spin-
averaged relaxation time, $\lambda_{\rm sd}$ is the spin diffusion length,
defined as $\lambda_{\rm sd}^{-2}=(l_{\uparrow}^{-2}+l_{\downarrow}^{-2})/2$,
where $l_{s}^{2}=v_{\rm F}^{2}\tau_{s}\tau_{\rm sf}/2$ and $v_{\rm F}$ is the
electron Fermi velocity. The difference between two spin-dependent scatterings
can be quantified by a spin asymmetry relaxation time parameter
$p_{\tau}=(\tau_{\uparrow}-\tau_{\downarrow})/(\tau_{\uparrow}+\tau_{\downarrow})$.
We can also define spin asymmetry of the spin-dependent conductivity as
$p_{\sigma}=(\sigma_{\uparrow}-\sigma_{\downarrow})/(\sigma_{\uparrow}+\sigma_{\downarrow})$.
In AFM metals with degenerate spin bands $p_{\tau}=p_{\sigma}$ while in spin
nondegenerate bands, such as FM metals, these two parameters can be different
[37].

The spin current density is determined from the formula ${\bf
j}_{s}={-e}(2\pi)^{-2}\int d^{2}\mathbf{k}\,f_{s}(\bf
r,\mathbf{k}){\bf{v}}_{\mathbf{k}}$, which upon using Eq. (11) leads to the
following relation,

${\bf j}_{s}(\bm{r})=\sigma_{s}\Big{(}{\bf E}-\nabla_{\bf
r}\mu_{s}(\bm{r})-\dfrac{e\tau_{s}}{m_{\rm eff}}{\bf E}\times\mathbf{B}_{\rm
em}^{s}(\bm{r})\Big{)},$ (13)

where $\sigma_{s}=(e^{2}/2h)(v_{F}k_{F}\tau_{s})$ is the spin-dependent
conductivity. The total charge and spin current densities are respectively
defined as,

$\displaystyle{\bf j}^{\rm ch}={\bf j}_{\uparrow}+{\bf j}_{\downarrow},$ (14a)
$\displaystyle{\bf j}^{\rm sp}={\bf j}_{\uparrow}-{\bf j}_{\downarrow}.$ (14b)

We are interested in the average of the spin accumulation along the transport
direction, i.e., $x$ direction. Equation (12) can then be rewritten as,

$\dfrac{d^{2}\overline{\delta\mu}(y)}{dy^{2}}-\dfrac{\overline{\delta\mu}(y)}{\lambda_{\rm
sd}^{2}}=\dfrac{e\tau E_{x}}{m_{\rm eff}}\dfrac{d\overline{B}_{\rm
em}(y)}{dy},$ (15)

where we defined $\overline{F}(y)=(2L)^{-1}\int_{-L}^{L}dxF({\bm{r}})$. The
general solution of this differential equation is the sum of homogeneous and
particular solutions. The homogeneous solution includes two unknown constants
that must be determined from the appropriate boundary conditions.

Similarly, using Eqs. (13) and (14), we find the following differential
equations for the transverse charge and spin current densities:

$\displaystyle\overline{j}^{\rm
ch}_{y}(y)=-\sigma_{0}\left(\dfrac{d\overline{\mu}(y)}{dy}+p_{\tau}\dfrac{d\overline{\delta\mu}(y)}{dy}\right)+\dfrac{e\tau\sigma_{0}}{m_{\rm
eff}}p_{\tau}E_{x}\overline{B}_{\rm em}(y),$ (16a)
$\displaystyle\overline{j}^{\rm
sp}_{y}(y)=-\sigma_{0}\left(p_{\tau}\dfrac{d\overline{\mu}(y)}{dy}+\dfrac{d\overline{\delta\mu}(y)}{dy}\right)+\dfrac{e\tau\sigma_{0}}{2m_{\rm
eff}}\big{(}1+p^{2}_{\tau})E_{x}\overline{B}_{\rm em}(y),$ (16b)

where $\sigma_{0}=\sigma_{\uparrow}+\sigma_{\downarrow}$ is the total
longitudinal conductivity of the AFM metal and
$\mu=(\mu_{\uparrow}+\mu_{\downarrow})/2$ is the spin-averaged chemical
potential.

To solve the differential equations (15) and (16), and to find expressions for
the spin accumulation, spin Hall current, and charge Hall voltage in the
assumed stripe geometry, we need now to use the appropriate boundary
conditions. Since the AFM stripe has a finite width with an open boundary
condition, the transverse component of the charge current density must be
zero, $\overline{j}^{\rm ch}_{y}(y)=0$. Consequently, using Eq. (16a) one
finds the average transverse electric field and Hall voltage,

$\displaystyle\overline{E}_{y}(y)$
$\displaystyle=-\dfrac{d{\overline{\mu}}(y)}{dy}=p_{\tau}\left(\dfrac{d\overline{\delta\mu}(y)}{dy}-\dfrac{2e\tau}{m_{\rm
eff}}E_{x}\overline{B}_{\rm em}(y)\right),$ (17a) $\displaystyle V_{H}$
$\displaystyle=\int_{-w}^{w}\overline{E}_{y}(y)dy.$ (17b)

As it is evident from the above expressions, the Hall voltage is finite only
if there is an asymmetry between the relaxation time of up and down spins,
i.e., $p_{\tau}\neq 0$. This emergent Hall voltage consists of two
contributions. The first term on the right-hand side of Eq. (17a) is
proportional to the spatial variation of the spin accumulation, and the second
term on the right-hand side of Eq. (17a) is directly proportional to the
emerging magnetic field.

Figure 1: Spatial variation of the spin accumulation, generated by the
skyrmion-induced real-space Berry curvature, Eq. (18), for different spin-
diffusion lengths. The spin accomulation is normalized to $\mu_{0}=(e\tau
E_{x}B_{0}r_{\rm sk})/m_{\rm eff}$ and we set $w=3r_{\rm sk}$. The inset shows
the profile of the skyrmion-induced magnetic field, Eq. (9). Figure 2: Spatial
variation of the topological spin Hall current density, Eq. (19), generated by
the skyrmion-induced real-space Berry curvature, for several spin diffusion
lengths. The spin current is normalized to
$j_{0}=(1-p_{\tau}^{2})(e\sigma_{0}\tau{E_{x}}B_{0})/m_{\rm eff}$ and we set
$w=3r_{\rm sk}$.

To find the spin current density and the Hall voltage, we need to know the
spatial dependence of the average spin accumulation $\overline{\delta\mu}(y)$.
Now, we should implement other boundary conditions. The spin current should be
zero at both edges, $\overline{j}^{\rm sp}_{y}(\pm w)=0$. Using this boundary
condition, together with Eqs. (15) and (16b), we find the nonequilibrium spin
accumulation profile as,

$\displaystyle\overline{\delta\mu}(y)=$
$\displaystyle\dfrac{e\tau{E_{x}}}{2m_{\mathrm{eff}}}\Bigg{(}\dfrac{\sinh(\frac{y}{\lambda_{\rm{sd}}})\exp(\dfrac{-w}{\lambda_{\rm{sd}}})}{\cosh(\dfrac{w}{\lambda_{\rm{sd}}})}\bigintssss_{-w}^{+w}\overline{B}_{\rm
em}(\tilde{y})\exp(\dfrac{\tilde{y}}{\lambda_{\rm{sd}}})d\tilde{y}$ (18)
$\displaystyle+\bigintssss_{-w}^{+w}\overline{B}_{\rm
em}(\tilde{y})\dfrac{y-\tilde{y}}{|y-\tilde{y}|}\exp(\dfrac{-|y-\tilde{y}|}{\lambda_{\rm{sd}}})d\tilde{y}\Bigg{)},$

The first term on the right-hand side of this expression is the homogeneous
solution of Eq. (15), and the second term is its particular solution.

Having the spin accumulation, Eq. (18), we find the spin Hall current density
from Eq. (16b),

$\displaystyle\overline{j}^{\rm
sp}_{y}(y)=\sigma_{0}(1-p_{\tau}^{2})\dfrac{e\tau E_{x}}{m_{\mathrm{\rm
eff}}}\Bigg{[}\dfrac{1}{2\lambda_{\rm
sd}}\Bigg{(}\bigintssss_{-w}^{+w}\overline{B}_{\rm
em}(\tilde{y})\exp(-\dfrac{|y-\tilde{y}|}{\lambda_{\rm sd}})d\tilde{y}$ (19)
$\displaystyle-\dfrac{\cosh(\frac{y}{\lambda_{\rm
sd}})\exp(\dfrac{-w}{\lambda_{\rm sd}})}{\cosh(\dfrac{w}{\lambda_{\rm
sd}})}\bigintssss_{-w}^{+w}\overline{B}_{\rm
em}(\tilde{y})\exp(\dfrac{\tilde{y}}{\lambda_{\rm
sd}})d\tilde{y}\Bigg{)}+\overline{B}_{\rm em}(y)\Bigg{]}.$

This expression describes the spatial dependence of the spin Hall current
density in the AFM stripe in the presence of an emerging field of a skyrmion.
This transverse spin current has two contributions. There is an extrinsic
contribution arising from the spin accumulation gradient (the first two
integrals) and an intrinsic contribution from the emergent magnetic field (the
last term).

##  V Spin accumulation, spin current density, and Hall voltage

The spin accumulation and spin current density are calculated from Eqs. (18)
and (19), respectively. The integrals cannot be calculated analytically and
thus we integrate them numerically. Figure 1 shows the spatial variation of
the nonequilibrium spin accumulation, Eq. (18), for various spin diffusion
lengths. The inset shows the profile of the emerging magnetic fields, Eq. (9).

The spin accumulation displays a nonmonotonic spatial behaviour, rising within
the first half of the skyrmion radius but declining as it extends towards the
edges. Beyond the skyrmion region, the spin accumulation gradually vanishes
for shorter diffusion lengths, whereas it reaches a constant value for longer
diffusion lengths. Reducing the spin diffusion length results in a decrease in
the spin accumulation amplitude as one may expect.

The associated topological spin Hall current flowing across the stripe is
presented in Fig. 2 for various spin diffusion lengths. The spin Hall current
reduces dramatically outside the skyrmion core, and eventually vanishes
outside the skyrmion. Inside the skyrmion core, the amplitude of the spin
current increases with reducing the spin diffusion length, while it decreases
away from the skyrmion core.

Finally, in Fig. 3 we plot the Hall voltage as a function of the spin
diffusion length for the indicated values of the spin asymmetry parameter
$p_{\tau}$. When the spin asymmetry of the relaxation time increases, the
amplitude of the Hall voltage increases as well. In turn, an increase in the
spin diffusion length leads to a decrease in the Hall voltage.

Figure 3: The Hall voltage as a function of the spin diffusion length for
indicated values of the spin asymmetry of the relaxation time. The Hall
voltage is normalized to $V_{0}={h\tau E_{x}}/({2Lm_{\rm eff}})$.

##  VI Summary and concluding remarks

We computed the topological charge and spin Hall effects as well as the spin
accumulation in a compensated AFM system, arising from a real-space Berry
curvature induced by a skyrmion. We calculated these quantities in the
diffusive regime using the semiclassical Boltzmann formalism. Our model
describes a single skyrmion in either a metallic AFM thin film or at the
interface of an AFM insulator and a metal. We considered both spin-dependent
scattering and spin-flip mechanisms in our calculations. We found both
intrinsic and extrinsic contributions to the spin Hall effect and spin
accumulation, and showed that the spin Hall current vanishes rapidly outside
the skyrmion. On the other hand, the spin accumulation, which is a measurable
quantity, can be finite in systems with the spin diffusion length larger than
the skyrmion size. In addition, we showed that the Hall voltage can be finite
in the presence of a spin-dependent relaxation time in such systems. Direct
detection of skyrmions is a challenge because of the absence of net
magnetization in compensated AFM systems. We argue that the electrical
detection of AFM skyrmions is possible by measuring the spin accumulation
and/or Hall voltage.

## Acknowledgements

This work has been supported by the Norwegian Financial Mechanism 2014- 2021
under the Polish – Norwegian Research Project NCN GRIEG “2Dtronics” no.
2019/34/H/ST3/00515.

## References

  * [1] G. Tatara and H. Fukuyama, Phys. Rev. Lett. 78 (1997) 3773.
  * [2] P. M. Levy and S. Zhang, Phys. Rev. Lett. 79 (1997) 5110.
  * [3] J. B. A. N. van Hoof, K. M. Schep, A. Brataas, G. E. W. Bauer, and P. J. Kelly, Phys. Rev. B 59 (1999) 138.
  * [4] A. Brataas, G. Tatara, and G. E. W. Bauer, Phys. Rev. B 60 (1999) 3406.
  * [5] E. Šimánek, Phys. Rev. B 63 (2001) 224412.
  * [6] M. Tsoi, R. E. Fontana, and S. S. P. Parkin, Appl. Phys. Lett. 83 (2003) 2617-2619.
  * [7] M. Yamanouchi, D. Chiba, F. Matsukura, and H. Ohno, Nature 428 (2004) 539–542.
  * [8] G. Tatara and H. Kohno, Phys. Rev. Lett. 92 (2004) 086601.
  * [9] G. Tatara, H. Kohno, and J. Shibata, Phys. Rep. 468 (2008) 213-301.
  * [10] A. C. Swaving and R. A. Duine, J. Phys. Condens. Matter 24 (2011) 024223.
  * [11] E. G. Tveten, A. Qaiumzadeh, O. A. Tretiakov, and A. Brataas, Phys. Rev. Lett. 110 (2013) 127208.
  * [12] Y. Yamane, J. Ieda, and J. Sinova, Phys. Rev. B 93 (2016) 180408.
  * [13] Y. Araki and K. Nomura, Phys. Rev. B 96 (2017) 165303.
  * [14] J. H. Zheng, A. Brataas, M. Kläui, and A. Qaiumzadeh, Phys. Rev. B 102 (2020) 184413.
  * [15] D. Liang, J. P. DeGrave, M. J. Stolt, Y. Tokura, and S. Jin, Nat. Commun. 6 (2015) 8217 .
  * [16] H. Wang, Y. Dai, G. M. Chow, and J. Chen, Prog. Mater. Sci. 130 (2022) 100971.
  * [17] T. Schulz et al. Nat. Phys. 8 (2011) 136804.
  * [18] N. Nagaosa and Y.Tokura, Nat. Nanotech. 8 (2013) 899–911.
  * [19] J. Iwasaki, M. Mochizuki, and N. Nagaosa, Nat. Commun. 4 (2013) 1463.
  * [20] K. Litzius et al. Nat. Phys. 13 (2017) 170-175.
  * [21] W. Jiang et al. Nat. Phys. 13 (2017) 162-169.
  * [22] I. A. Ado, O. A. Tretiakov, and M. Titov, Phys. Rev. B. 95 (2017) 094401.
  * [23] O. A. Tretiakov, D. Clarke, G.-Wei Chern, Y. B. Bazaliy, and O. Tchernyshyov, Phys. Rev. Lett. 100 (2008) 127204.
  * [24] D. J. Clarke, O. A. Tretiakov, G.-W. Chern, Ya. B. Bazaliy, and O. Tchernyshyov, Phys. Rev. B 78, 134412 (2008).
  * [25] J. Zang, M. Mostovoy, J. Han, and N. Nagaosa, Phys. Rev. Lett. 107 (2011) 136804.
  * [26] H. B. Braun, Adv. Phys. 61 (2012) 1-116.
  * [27] R. Tomasello et al. Sci. Rep. 4 (2014) 6784.
  * [28] A. A. Thiele, Phys. Rev. Lett. 30 (1973) 230.
  * [29] J. C. Martinez and M. B. A. Jalil, New J. Phys. 18 (2016) 033008.
  * [30] J. Ye, Y. B. Kim, A. J. Millis, B. I. Shraiman, P. Majumdar, and Z. Tešanović, Phys. Rev. Lett. 83 (1999) 3737.
  * [31] P. Bruno, V. K. Dugaev, and M. Taillefumier, Phys. Rev. Lett. 93 (2004) 096806.
  * [32] A. Neubauer, C. Pfleiderer, B. Binz, A. Rosch, R. Ritz, P. G. Niklowitz, and P. Böni, Phys. Rev. Lett. 102 (2009) 186602.
  * [33] L. Minhyea, W. Kang, Y. Onose, Y. Tokura, and N. P. Ong, Phys. Rev. Lett. 102 (2009) 186601.
  * [34] Y. Ohuchi, Y. Kozuka, M. Uchida, K. Ueno, A. Tsukazaki, and M. Kawasaki, Phys. Rev. B 91 (2015) 245115.
  * [35] K. S. Denisov, I. V. Rozhansky, N. S. Averkiev, and E. Lähderanta, Phys. Rev. Lett. 117 (2016) 027202.
  * [36] P. B. Ndiaye, C. A. Akosa, and A. Manchon, Phys. Rev. B 95 (2017) 064426.
  * [37] S. S.-L. Zhang and O.-L. Heinonen, Phys. Rev. B. 97 (2018) 134401.
  * [38] K. Robert and J. M. Luttinger,  Phys. Rev. 95 (1954) 1154 .
  * [39] N. Nagaosa, J. Sinova, S. Onoda, A. H. MacDonald, and N. P. Ong,  Rev. Mod. Phys. 82 (2010) 1539.
  * [40] H. Chen, Q. Niu, and A. H. MacDonald,  Phys. Rev. Lett. 112 (2014) 017205.
  * [41] S. E. Barnes and S. Maekawa, Phys. Rev. Lett. 98 (2007) 246601.
  * [42] V. K. Dugaev, A. Crépieux, and P. Bruno,  Phys. Rev. B 64 (2001) 104411.
  * [43] A. Crépieux and P. Bruno, Phys. Rev. B 64 (2001) 014416.
  * [44] D. Culcer, A. MacDonald, and Q. Niu,  Phys. Rev. B 68 (2003) 045327.
  * [45] C. Jin, C. Song, J. Wang, and Q. Liu, Appl. Phys. Lett. 109 (2016) 182404.
  * [46] J. Barker and O. A. Tretiakov, Phys. Rev. Lett. 116 (2016) 147203.
  * [47] X. Zhang, Y. Zhou, and M. Ezawa, Nat. Commun. 7 (2016) 10293.
  * [48] X. Zhang, Y. Zhou, and M. Ezawa, Sci. Rep. 6 (2016) 24795.
  * [49] R. Khoshlahni, A. Qaiumzadeh, A. Bergman, and A. Brataas, Phys. Rev. B 99 (2019) 054423.
  * [50] O. A. Tretiakov, Woodhead Publishing Series in Electronic and Optical Materials (2021) 333-345.
  * [51] B. Göbel, I. Mertig, and O. A. Tretiakov, Phys. Rep. 895 (2021) 1-28.
  * [52] O. J. Amin et al. Nat. Nanotechnol. 18 (2023) 849–853 .
  * [53] R. Cheng and Q. Niu, Phys. Rev. B 86 (2012) 245118.
  * [54] P. M. Buhl, F. Freimuth, S. Blügel, and Y. Mokrousov,  Phys. Status Solidi RRL 11 (2017) 1700007.
  * [55] B. Göbel, A. Mook, J. Henk, and I. Mertig, Phys. Rev. B 96 (2017) 060406.
  * [56] C. A. Akosa, O. A. Tretiakov, G. Tatara, and A. Manchon, Phys. Rev. Lett. 121 (2018) 097204.
  * [57] K. Nakazawa, K. Hoshi, J. J. Nakane, J.-I. Ohe, and H. Kohno,  arXiv:2304.02850 (2023).
  * [58] B. Qiang et al.  Jpn. J. Appl. Phys. 61 (2022) 120901.
  * [59] H. B. M. Saidaoui, X. Waintal, and A. Manchon, Phys. Rev. B 95 (2017) 134424.
  * [60] A. Qaiumzadeh, I. A. Ado, R. A. Duine, M. Titov, and A. Brataas, Phys. Rev. Lett. 120 (2018) 197202.
  * [61] M. M. S. Barbeau, M. Titov, M. I. Katsnelson, and A. Qaiumzadeh, Phys. Rev. Research 5 (2023) L022065.
  * [62] K. S. Denisov, I. V. Rozhansky, N. S. Averkiev, and E. Lähderanta, Phys. Rev. B 98 (20) 195439.
  * [63] G. E. Volovik, J. Phys. C 20 (1987) L83.
  * [64] Y. Tserkovnyak and M. Mecklenburg, Phys. Rev. B. 77 (2008) 134407.
  * [65] S. Zhang and S. S.-L. Zhang, Phys. Rev. Lett. 102 (2009) 086601.
  * [66] C. A. Akosa, _Spin Transport in Ferromagnetic and Antiferromagnetic Textures_ , PhD thesis, King Abdullah University of Science and Technology (KAUST), Thuwal, Saudi Arabia (2016).
  * [67] N. W. Ashcroft and N. D. Mermin, _Solid State Physics_ Harcourt College Publishers, San Diego (1975).
  * [68] Y. Qi and S. Zhang, Phys. Rev. B 67 (2003) 052407.
  * [69] S. S.-L. Zhang and S. Zhang, Phys. Rev. B 82 (2010) 184423.
  * [70] T. Valet and A. Fert, Phys. Rev. B 48 (1993) 7099.
  * [71] A. Fert, T. Valet, and J. Barnaś, Appl. Phys. 75 (1994) 6693.


# Viscoelastic increase of detachment stress of a rigid punch from adhesive
soft viscoelastic layers

Ali Maghami  Michele Ciavarella  Antonio Papangelo 111Corresponding author
Email address<EMAIL_ADDRESS>(Antonio Papangelo).

###### Abstract

The problem of the detachment of a sufficiently large flat indenter from a
plane adhesive viscoelastic strip of thickness “b” is studied. For any given
retraction speed, three different detachment regimes are found: (i) for very
small “b” the detachment stress is constant and equal to the theoretical
strength of the interface, (ii) for intermediate values of “b” the detachment
stress decays approximately as b-1/2, (iii) for thick layers a constant
detachment stress is obtained corresponding to case the punch is detaching
from a halfplane. By using the boundary element method a comprehensive
numerical study is performed which assumes a linear viscoelastic material with
a single relaxation time and a Lennard-Jones force-separation law. Pull-off
stress is found to consistently and monotonically increase with unloading
rate, but to be almost insensitive to the history of the contact. Due to
viscoelasticity, unloading at high enough retraction velocity may allow
punches of macroscopic size to reach the theoretical strength of the
interface. Finally, a corrective term in Greenwood or Persson theories
considering finite size effects is proposed. Theoretical and numerical results
are found in very good agreement.

###### keywords:

Viscoelasticity , Crack propagation , Thin layer , Detachment force.

††journal: …

[inst1]organization=Department of Mechanics, Mathematics and Management,
Politecnico di Bari,addressline=via Orabona 4, city=Bari, postcode=70125,
country=Italy

[inst2]organization=Department of Mechanical Engineering, Hamburg University
of Technology,addressline=Am Schwarzenberg-Campus 1, city=Hamburg,
postcode=21073, country=Germany

##  1 Introduction

Soft materials are of great interest in the scientific community as for their
applicability in many engineering fields ranging from the automotive sector
[1], biomechanics [2], soft robotics [3, 4], manipulators [5], tires grip[1],
sealing of syringes[6], finger-touch-pad interactions [7, 8], soft tissue
adhesion for regenerative medicine [9, 10]and pressure-sensitive adhesives
[11]. With soft polymers macroscopic adhesion due to van der Waals adhesive
interactions remains strong [12], whereas in hard materials it is easily
canceled by the inevitable surface roughness.

In many of these applications the bond strength is a crucial mechanical
property and it is often quantified by measuring the apparent adhesion
strength as given by the maximum pulling force per unit area in a tensile bond
test [13]. Peng et al. [13] have already elucidated how the critical pull-off
force of a flat rigid axisymmetric punch adhered to an elastic film of finite
thickness depends by two dimensionless parameters. The former shows a
transition from uniform detachment (DMT-like) to crack-like propagation (JKR-
like), while the latter is a correction factor due to finite thickness of the
film.

However, in many of the applications mentioned above, soft materials
(polymers, elastomers) are employed, which are known to be viscoelastic, i.e.
they exhibit a frequency-dependent modulus and dissipation [14], and this
complicates their mechanical adhesive behaviour. Numerous experiments in
steady state conditions have shown that the apparent surface energy
$\Delta\gamma$, i.e. the energy per unit area needed to separate two ideally
parallel surfaces, during the crack opening is related to the crack speed $V$
through a power law function [15, 16, 17, 18], commonly referred as the Gent
and Schulz empirical law [19]

$\frac{\Delta\gamma\left(V\right)}{\Delta\gamma_{0}}=1+\left(\frac{V}{V_{ref}}\right)^{n},$
(1)

where $\Delta\gamma_{0}$ is the adiabatic surface energy (or thermodynamic
work of adhesion), $V$ is the crack velocity,
$V_{ref}=\left(ka_{T}^{n}\right)^{-1}$ and $k,n$ are constants with $0<n<1$
and $a_{T}$ is the WLF factor to translate viscoelastic modulus results at
various temperatures $T$ [20]. In its simplest form Gent and Schulz empirical
law Eq. (1) is generally a good phenomenological model for opening cracks,
while for closing cracks a reduced apparent surface energy $\Delta\gamma$ is
observed which generally shows the reciprocal of that law [21, 22, 23, 24, 25,
26]. There is no indication in the law of a limit enhancement, nor it is clear
how far it can be used for transient conditions. Finally, and most
importantly, in this empirical law, there is no indication on possible size
effects, i.e. on how the parameters of the law should be affected by geometry,
and $\frac{\Delta\gamma\left(V\right)}{\Delta\gamma_{0}}>1$ which instead we
shall find is not always true, even for advancing crack.

There are two main approaches which have been attempted to capture more
fundamentally the propagation of viscoelastic crack propagation. One is based
on the description of the processes occurring at the crack tip through a
cohesive zone model [27, 25, 26, 28], and as such is rather general as it can
take into account of initiation of the crack, transient propagation, and
steady state. Also it may show a transition to a cohesive rupture for small
enough cracks (what Peng et al. [13] call uniform DMT-like detachment in their
case, see also [29, 30, 31]), in principle it may be generalized to rate-
dependent cohesive laws, and to non linear materials. In practice, it attempts
to model the real processes occurring at the crack tip. The other, developed
by Persson and coauthors [22, 23, 24], takes an “energy-based” approach and is
restricted to linear materials and to steady state conditions. It is derived
by equating the power input in the system with the power dissipated by
viscoelastic losses and by the rate at which energy is spent to create new
surfaces. The energy-based approach finds different results for finite size
systems where it seems to show non-monotonic
$\Delta\gamma\left(V\right)/\Delta\gamma_{0}$ and also
$\Delta\gamma\left(V\right)/\Delta\gamma_{0}>1$ [23], which is in contrast to
the limit case of cohesive failure where
$\Delta\gamma\left(V\right)/\Delta\gamma_{0}$ can decrease down to zero for
very small crack (see [29]), as we shall discuss in details here with respect
to our geometry. Only for a semi-infinite system and linear material both
approaches yield a very similar monotonically increasing behavior of
$\Delta\gamma\left(V\right)/\Delta\gamma_{0}$ with respect to $V$ up to the
theoretical “high-frequency” limit of
$\Delta\gamma/\Delta\gamma_{0}=E_{\infty}/E_{0}$, where $E_{\infty}$ and
$E_{0}$ respectively represent the glassy (high frequency) and the rubbery
(low frequency) modulus of the viscoelastic material [32]111In their basic
form, both approaches consider $\Delta\gamma_{0}$ to be an intrinsic material
property and therefore rate-independent, as we shall also assume here. .
Furthermore, the cohesive model has been applied modelling of bi-materials
crack (one elastic, one viscoelastic) [33] showing transient dissipation can
be arbitrarily large while loads remain finite and hence dissipation should
not be taken as an indication of true fracture energy [34].

Recently, it has been shown that, depending on the indenter geometry, the
loading history may or may not affect the detachment force. For a Hertzian
indenter, Afferrante and Violano [30, 31] have shown that the loading velocity
and the maximum indentation reached during the loading phase can significantly
influence the pull-off force. One the other hand, Papangelo and Ciavarella
[29] found that when considering an axisymmetric flat punch, the loading
history has a relatively weak effect, and the primary determinant of the pull-
off force is largely the unloading velocity. This is partly related to the
fact that, for a given preload, the contact area achieved at maximum
indentation strongly depends on the loading history only for a Hertzian
indenter, while it remains fixed for a flat punch. This interpretation is
somehow confirmed by the recent work of Muller et al. [35] which further
considers the contact of a flat punch indenter with superimposed small scale
roughness. Their study focuses on the significant hysteresis that is observed
during crack closure and opening which is obtained by the concurrent presence
of viscoelasticity and adhesion. They find that small scale roughness can
indeed leads to loacal jumps -in and -out of contact, which causes the
dependence of the detachment force on the preload.

While previous works have focused on the detachment from viscoelastic semi-
infinite substrates [29, 36, 37], the influence of the layer thickness on the
detachment process has been mostly overlooked. Nevertheless, the latter has a
broad interest in engineering applications where often a thin layer of
viscoelastic material is used, as in Figure 1, as well as to exploit the thin
layer testing geometries, like in Peng et al. [13] which could permit to
extract cohesive properties as well as surface energy properties.

Figure 1: Practical engineering applications where thin polymeric layers are
used: (a) rubber tape, (b) polymeric carpet, (c) gasket, (d) soft gripper for
delicate handling and manipulation (UCSD Jacobs School of Eng., CC BY-NC-SA
4.0).

Hence, in the present study, we focus on the detachment problem of a flat
punch indenting a thin adhesive viscoelastic layer of finite thickness $b$ in
plane strain conditions. The remainder of the paper is organized as follows.
In Section 2 the elastic solution for halfplane contact is recalled. In
Section 3, the case of a thin layer is considered leveraging on the the “thin
strip” assumption by Johnson [38] and applying the Griffith energy balance.
This serves to determine the limiting solutions within which the viscoelastic
results should be confined. In Section 4 the numerical scheme is introduced
which is based on the boundary element method and it assumes a standard linear
model for the viscoelastic layer and a Lennard-Jones force-separation law for
the contact interactions, which is rate-independent. In Section 5 the
numerical results are presented and compared with Greenwood theory for
viscoelastic crack propagation [21] which we extend for finite size effects.
In Section 6 the conclusions are drawn. The study is carried on with
particular emphasis on the effect of (i) the loading history, (ii) the layer
thickness and (iii) the unloading velocity on the pull-off force and on the
effective adhesive energy.

##  2 Detachment from a halfplane

Let us consider the plane contact problem of a flat punch of semi-width $a$
indenting an elastic adhesive frictionless halfplane with Young modulus $E$
and Poisson ratio $\nu$. By applying the Griffith energy balance, the pull-off
force [39, 40] is given by

$P_{po}=L\sqrt{2\pi E^{\ast}\Delta\gamma_{0}a}\;,$ (2)

where $L$ is the layer width, $E^{\ast}=E/\left(1-\nu^{2}\right)$ is the plane
strain elastic modulus. Hence the mean interfacial stress at pull-off is

$\overline{\sigma}_{po}=\sqrt{\frac{\pi E^{\ast}\Delta\gamma_{0}}{2a}}\;,$ (3)

which has the classical Linear Elastic Fracture Mechanics (LEFM) square-root
dependence with respect to the punch semiwidth $a$. Overbar indicates here the
mean value. This implies that smaller punches have a higher pull-off stress,
potentially reaching the theoretical strength (or the cohesive strength) of
the interface, denoted as $\sigma_{0}$. This is typically observed for punches
with a semi-width less than the following typical fracture length

$a_{0}=\frac{\pi}{2}\frac{E^{\ast}\Delta\gamma_{0}}{\sigma_{0}^{2}}\;.$ (4)

In what follows in the paper we shall assume that the punch size is $a>>a_{0}$
as we are interested in the transitions due to the finite size of the layer
rather than the size of the punch. For the latter effect the reader is
referred to Ref. [29].

Hence, in dimensionless form, we have the following relations

$\widehat{\overline{\sigma}}_{po}=\frac{\overline{\sigma}_{po}}{\sigma_{0}}=\frac{1}{\sqrt{a/a_{0}}}\;,$
(5)
$\frac{a_{0}}{h_{0}}=\frac{9\sqrt{3}\pi}{32\Sigma_{0}}\approx\frac{1.53}{\Sigma_{0}}\;,$
(6)

where we have assumed a Lennard-Jones force-separation law, for which
$\Delta\gamma_{0}=\alpha_{LJ}h_{0}\sigma_{0}$, where $h_{0}$ is the range of
interaction and $\alpha_{LJ}=\frac{9\sqrt{3}}{16}$ is a constant,
$\widehat{a}=a/h_{0}$, $\Sigma_{0}=\sigma_{0}/E^{\ast}$ is usually in the
range of $\left[0.1\div 1\right]$ for soft polymers [41, 42, 40], implying
$a_{0}$ to be 1 to 10 times higher than the range of attractive forces. For a
true crystal, this would imply a range of few nanometers, but for soft
materials, the range of adhesive forces may be larger. If a viscoelastic
material with relaxed modulus $E_{0}$ and instantaneous modulus $E_{\infty}$
is considered, then in the limit of very slow and very fast unloading rate we
will have

$\left\\{\begin{array}[c]{cc}\widehat{\overline{\sigma}}_{po}=\sqrt{\frac{9\sqrt{3}\pi}{32\Sigma_{0}\widehat{a}}}=\frac{1}{\sqrt{a/a_{0}}}\;;&\qquad\text{slow
limit }E=E_{0}\;,\\\
\widehat{\overline{\sigma}}_{po}=\sqrt{\frac{9\sqrt{3}\pi}{32\Sigma_{0}k\widehat{a}}}=\frac{1}{\sqrt{ka/a_{0}}}\;;&\qquad\text{fast
limit }E=E_{\infty}\;,\end{array}\right.$ (7)

where $k=E_{0}/E_{\infty}$ and $\widehat{a}=a/h_{0}$.

##  3 Detachment from a thin layer

If the substrate has a finite thickness, it is necessary to consider the
effect of thickness in the analysis. Hence, here we focus on the plane contact
problem of a flat punch with a semi-width of $a$ indenting an adhesive layer
with a thickness of $b$ (Fig. 2). We first focus on the linear elastic
solution, and then we will provide the limiting solutions for the viscoelastic
problem based on the elastic formulation.

###  3.1 Elastic layer

Let us consider the layer in plane strain and supported by a rigid foundation.
In the following the case of frictionless contact between the layer and the
rigid substrate is considered while the correction due to the Poisson effect
for the case of a layer perfectly bonded to the substrate is shown in the
Appendix-I.

Figure 2: On the left is a sketch of a flat punch being loaded on a
viscoelastic adhesive layer of thickness $b$. The numerical implementation
employs the Lennard-Jones force-separation law, while the viscoelastic
material is represented using a standard linear model, as depicted in the
lower-right panel.

Following Johnson [38], we assume that plane sections remain plane upon
loading. Hence for the case of no friction between the layer and the rigid
substrate, the load $P$ and the corresponding elastic strain energy $U_{E}$
stored in the layer are

$\displaystyle P$ $\displaystyle=-2aLE^{\ast}\frac{\delta}{b}\;,$ (8)
$\displaystyle U_{E}$ $\displaystyle=aLE^{\ast}\frac{\delta^{2}}{b}\;,$ (9)

where $\delta$ is the indentation considered positive when the flat punch is
approaching the substrate, consequently, $P$ is positive when tensile. At
unloading, the Griffith energy balance requires the elastic strain energy
released per unit area to be equal to the surface energy. Hence, assuming
detachment occurs immediately, we have the following relations

$\frac{1}{2L}\frac{\partial U_{E}}{\partial
a}=\Delta\gamma\rightarrow\left\\{\begin{array}[c]{c}\delta_{po}=-\sqrt{\frac{2b\Delta\gamma}{E^{\ast}}}\;,\\\
\overline{\sigma}_{po}=\sqrt{\frac{2E^{\ast}\Delta\gamma}{b}}\;,\end{array}\right.$
(10)

where, $\delta_{po}$ and $\overline{\sigma}_{po}$ are the indentation and the
average interfacial stress at pull-off respectively. Notice that the pull-off
stress depends on the layer thickness as $\overline{\sigma}_{po}\propto
b^{-1/2}$, hence it is possible to define a characteristic thickness $b_{0}$
of the substrate where $\overline{\sigma}_{po}$ reaches the theoretical
interfacial strength $\sigma_{0}$, i.e.

$b_{0}=\frac{2E^{\ast}\Delta\gamma}{\sigma_{0}^{2}}=\frac{4}{\pi}a_{0}\approx
1.27a_{0}\;,$ (11)

which is of the same order of magnitude of $a_{0}$. Experiments with PDMS
elastomers in Peng et al. [13], show that this $b_{0}$ is of the order of
$0.1$ mm, where clearly their loading rate corresponds to a certain effective
elastic modulus. On the other hand, in the limit of a very thick layer, we
should obtain the half-plane solution, for which we can utilize Eq. (3) and
Eq. (10) to determine a length scale

$b_{1}=\frac{4a}{\pi}\;,$ (12)

with the meaning that for substrates thicker than $b_{1}$ one should
anticipate the half-plane behaviour. Notice that, while $b_{0}$ is a
characteristic lengthscale that depends on the material and interfacial
properties, $b_{1}$ depends on the punch semi-width. Overall as indicated by
Eq. (11) and Eq. (12), and as illustrated in Fig 3, we identify the following
three regimes

$\left\\{\begin{array}[c]{cr}\overline{\sigma}_{po}=\sigma_{0}\;,&\qquad
b<b_{0}\\\
\overline{\sigma}_{po}=\sqrt{\frac{2E^{\ast}\Delta\gamma}{b}}\;,&\qquad
b_{0}\leq b\leq b_{1}\\\ \overline{\sigma}_{po}=\sqrt{\frac{\pi
E^{\ast}\Delta\gamma}{2a}}\;,&\qquad b>b_{1}\end{array}\right.$ (13)

or, in dimensionless form,

$\left\\{\begin{array}[c]{cr}\widehat{\overline{\sigma}}_{po}=1\;,&\qquad
b/a_{0}<\frac{4}{\pi}\\\
\widehat{\overline{\sigma}}_{po}=\sqrt{\frac{9\sqrt{3}}{8\Sigma_{0}\widehat{b}}}=\sqrt{\frac{4}{\left(b/a_{0}\right)\pi}}\;,&\qquad\frac{4}{\pi}\leq
b/a_{0}\leq\frac{4}{\pi}\frac{a}{a_{0}}\\\
\widehat{\overline{\sigma}}_{po}=\sqrt{\frac{9\sqrt{3}\pi}{32\Sigma_{0}\widehat{a}}}=\frac{1}{\sqrt{a/a_{0}}}\;.&\qquad
b/a_{0}>\frac{4}{\pi}\frac{a}{a_{0}}\end{array}\right.$ (14)

###  3.2 Limiting solutions for a viscoelastic layer

Let us assume that the layer is constituted by a viscoelastic material with
relaxed Young modulus $E_{0}$ and instantaneous Young modulus $E_{\infty}$ so
that $k=E_{0}/E_{\infty}$. In the limit of very slow/very fast unloading rate
the substrate behaves as elastic. Thus, for the case of no friction between
the substrate and the layer, according to Eq. (10), one can anticipate the
following two scenarios

$\left\\{\begin{array}[]{cc}\overline{\sigma}_{po}=\sqrt{\frac{2E_{0}^{\ast}\Delta\gamma}{b}}\;,&\qquad\text{\leavevmode\ltxml@oqmark@open\textquotedblleft\penalty
10000\hskip-0.0002pt\hskip 0.0002ptvery
slow\textquotedblright\ltxml@oqmark@close{}}\\\
\overline{\sigma}_{po}=\sqrt{\frac{2E_{\infty}^{\ast}\Delta\gamma}{b}}\;,&\qquad\text{\leavevmode\ltxml@oqmark@open\textquotedblleft\penalty
10000\hskip-0.0002pt\hskip 0.0002ptvery
fast\textquotedblright\ltxml<EMAIL_ADDRESS>(15)

or in dimensionless form

$\left\\{\begin{array}[]{cc}\widehat{\overline{\sigma}}_{po}=\sqrt{\frac{4}{(b/a_{0})\pi}}\;,&\qquad\text{\leavevmode\ltxml@oqmark@open\textquotedblleft\penalty
10000\hskip-0.0002pt\hskip 0.0002ptvery
slow\textquotedblright\ltxml@oqmark@close{}}\\\
\widehat{\overline{\sigma}}_{po}=\sqrt{\frac{4}{k(b/a_{0})\pi}}\;,&\qquad\text{\leavevmode\ltxml@oqmark@open\textquotedblleft\penalty
10000\hskip-0.0002pt\hskip 0.0002ptvery
fast\textquotedblright\ltxml<EMAIL_ADDRESS>(16)

where, one should notice that for rapid unloading (very fast scenario) the
pull-off stress will reach the cohesive strength by the following value of the
substrate thickness:

${b_{0\infty}}=\frac{4{a_{0}}}{k\pi}\;.$ (17)

However, for a thick layer, the halfplane limit will be always obtained at
$b_{1}=4a/\pi,$ irrespective of the unloading rate. Figure 3 schematically
displays the elastic limits at low and high retraction speed that constitute
the bounds for the possible viscoelastic solutions.

Form Eq.s (12) and (17) it follows that if $a/a_{0}<E_{\infty}/E_{0}$ then
$b_{0\infty}>b_{1}$. In other words, if $a<a_{0}/k$ at a high enough
retraction velocity, it is possible to reach the adhesive strength of the
interface. It is easy to find elastomers with $E_{\infty}/E_{0}\simeq
10^{3}\div 10^{4}$ [43]. This implies that punches with semiwidth $a$ much
larger than $a_{0}$ can still reach the theoretical interfacial strength if
unloading is performed fast enough. Also, this holds for layer thickness. If
we assume that the experiments with PDMS in Peng et al. [13] were conducted
relatively slowly, when $b_{0}$ was of the order of $0.1$ mm, then clearly
${b_{0\infty}}$ could reach very large values, provided the punch is also
large enough. This peculiarity may be exploited in the future as a technique
to improve/enhance interfacial adhesion in micro-structured interfaces by
optimizing not only the micro-pillar geometry, but also the unloading
protocols.

Figure 3: Schematic representation of the three possible detachment regimes.

##  4 Numerical implementation of the adhesive viscoelastic contact problem

In this section, the numerical scheme used to solve the adhesive viscoelastic
contact problem is introduced. We utilized the boundary element method, which
necessitates the discretization solely of the interface. A similar code has
been used in previous works for solving both adhesive elastic [44] and
viscoelastic [29] axisymmetric contact problems, therefore, in this section,
we will focus on the essential adaptations required to tailor the
implementation for a plane viscoelastic strip.

The flat punch interacts with the viscoelastic layer according to the Lennard-
Jones 3-9 force-separation law defined as

$\sigma\left(h\right)=\frac{8\Delta\gamma}{3h_{0}}\left[\left(\frac{h_{0}}{h}\right)^{3}-\left(\frac{h_{0}}{h}\right)^{9}\right]\;,$
(18)

where $\sigma$ is the traction ($\sigma>0,$ when it is tensile), $h$ is the
interfacial gap and $h_{0}$ is the equilibrium distance. The theoretical
strength of the interface (maximum tensile stress) is then equal to
$\sigma_{0}=\Delta\gamma/(\alpha_{LJ}h_{0})$ and takes place at a separation
equal to $h=3^{1/6}h_{0}$. The gap is a function of the in-plane coordinate
$x$ as

$h(x)=-\delta+h_{0}+v\left(x\right)\;,$ (19)

where, $v\left(x\right)$ is the deflection of the viscoelastic layer with
respect to the origin $\left(0,0\right)$ ($v\left(x\right)$ is positive as
shown in Fig. 2). Here, Eq. (19) is solved numerically in a discrete manner at
the $N=M+1$ nodes, being $M$ the number of equally spaced elements with the
length of $c=2a/M$. Following Bentall and Johnson [45] we implemented the
method of overlapping triangles, i.e. for the $n$-th node the pressure is 0 at
node $x_{n-1}$, rises linearly to $p_{n}$ at node $x_{n}$ and then falls
linearly to 0 at node $x_{n+1}$, which gives overall a linear variation of the
contact pressure $p\left(x\right)$ over the considered domain. With respect to
the case of constant pressure elements, a picewise-linear distribution of
normal tractions produces a displacement field which is everywhere smooth and
continuous. Hence, according to Bentall and Johnson [45] the vertical
deflection at node $m$ of an elastic layer relatively to the origin
$\left(x,y\right)=\left(0,0\right)$ due to a triangular distribution of
pressure centered in $x_{n}$ is

$v_{m}=aB\frac{4}{\pi
E^{\ast}}\left[I_{A0}+I_{A}\left[m-n\right]+4zI_{AR}\left[m-n\right]\right]p_{n}\;,$
(20)

where $\left\\{m,n\right\\}$ are integers numbers, $p_{n}$ is the pressure
acting on the $n$-th node determined using Eq. 18 ( $p_{n}>0$ when it is
tensile), $B=b/a$, $z=c/4b=1/2BM$ and $I_{A0}$, $I_{A}$, $I_{AR}$ are the
following integral functions222Care should be taken when integrating $I_{A0}$
which converges slowly. The Appendix 3 of Bentall and Johnson [45] suggests a
convenient integration strategy we have also adopted. Notice that Bentall and
Johnson [45] contains a misprint as the second part in which $I_{A0}$ is split
up should be integrated over the interval $\left[\delta,+\infty\right]$.

$\displaystyle I_{A0}$
$\displaystyle=\frac{2}{z}\int_{0}^{\infty}\left(\frac{1-\cosh\beta}{\beta+\sinh\beta}\right)\frac{\sin^{2}\left(\beta
z\right)}{\beta^{3}}d\beta\;,$ (21) $\displaystyle I_{A}\left[m-n\right]$
$\displaystyle=-\frac{4}{z}\int_{0}^{\infty}\left(1+\frac{1-\cosh\beta}{\beta+\sinh\beta}\right)\frac{\sin^{2}\left(\beta
z\right)}{\beta^{3}}\sin^{2}\left(\beta z\left(m-n\right)\right)d\beta\;,$
(22) $\displaystyle I_{AR}\left[m-n\right]$
$\displaystyle=\int_{0}^{\infty}\frac{\sin^{2}\left(\eta\right)\sin^{2}\left(\eta\left(m-n\right)\right)}{\eta^{3}}d\eta,\qquad\eta=\beta
z\;.$ (23)

By applying the superposition principle, the normal deflection $v_{m}$ at node
$m$ due to a piecewise linear distribution of pressure can be written as

$v_{m}=\frac{1}{E^{\ast}}{\displaystyle\sum\limits_{n=1}^{N}}G_{mn}p_{n},$
(24)

where each column of the compliance matrix
$\\{(1/E^{*})\overline{\overline{G}}\\}_{NxN}$ corresponds to the displacement
field due to a unity triangular pressure centered at node $n$ being all the
other nodes unloaded. Therefore, the displacement field and, correspondingly,
the compliance matrix can be readily computed using Eq. (20). Once the elastic
solution is obtained, the displacement field of the viscoelastic layer
${v}\left(x,t\right)$, can be determined by the elastic-viscoelastic
correspondence principle in the form of Boltzmann integrals [46] as

${v}\left(x,t\right){=\frac{1}{E_{0}^{\ast}}\int
G}{\left(x,s\right)\int_{-\infty}^{t}c(t-\tau)\frac{dp(s,\tau)}{d\tau}d\tau
ds}\;,$ (25)

where $c(t)$ is the dimensionless creep compliance function, the strain
variation after an application of a constant unit stress, and, in our discrete
formulation, the Green function ${G}{\left(x,s\right)}$ is replaced by the
appropriate tensor $\\{\overline{\overline{G}}\\}_{NxN}$, so that the
viscoelastic nodal displacements $\\{{v}(t)\\}_{Nx1}$ at time $t$ are

$\\{{v}(t)\\}_{Nx1}\boldsymbol{=}\\{\overline{\overline{G}}\\}_{NxN}\ast\biggl{\\{}\frac{1}{E_{0}^{\ast}}{\int_{-\infty}^{t}c(t-\tau)\frac{dp}{d\tau}d\tau}\biggl{\\}}_{Nx1},$
(26)

where the symbol “$\ast$” stands for the row by column product. For the linear
viscoelastic material, the standard model is assumed with a single relaxation
time $\tau$, composed by a spring placed in series with an element constituted
by a dashpot and a spring in parallel (see Fig. 2), for which the
dimensionless creep compliance function is as follows

$c(t)=\left[1+\left(k-1\right)\exp\left(-\frac{t}{\tau}\right)\right]\;.$ (27)

being $\tau$ the relaxation time of the material. Hence, by using a sequential
time-marching continuation, we solved Eq.s (18,19,26), where at each time step
an iterative scheme is used to determine the equilibrium solution.

##  5 Results

Here the results of the numerical investigations are shown by using the
following dimensionless parameters

$\widehat{a}=\frac{a}{h_{0}};\;\;\;\;\widehat{\sigma}(x)=\frac{\sigma(x)}{\sigma_{0}};\;\;\;\;\widehat{\overline{\sigma}}=\frac{P}{2aL\sigma_{0}};\;\;\;\;\widehat{\delta}=\frac{\delta}{h_{0}};\;\;\;\;\widehat{t}=\frac{t}{\tau},$
(28)

and $\widehat{\overline{\sigma}}_{po}$ is the (dimensionless) average stress
at pull-off and is defined as
$\widehat{\overline{\sigma}}_{po}=\textrm{max}(\widehat{\overline{\sigma}})$.
If not stated differently, in our simulations we considered $M=200$,
$\Sigma_{0}=0.05$, and $k=0.1$.

###  5.1 History dependence

Viscoelastic materials typically exhibit a “history-dependent” response and
this tremendously affects the detachment force in Hertzian indenters [30, 31].
Hence, first we aimed to explore how different loading scenarios affect the
detachment characteristics of the flat indenter we considered, while keeping
the unloading rate constant. The simulations were carried out under
displacement control using a trapezoidal function (see inset in Fig. 4.a)
which main parameters are shown in Table 1. We define the dwell time as
$\widehat{t}_{dwell}=\widehat{t}_{2}-\widehat{t}_{1}$, the unloading rate
$\widehat{r}=(\widehat{\delta}_{load}-\widehat{\delta}_{unload})/(\widehat{t}_{3}-\widehat{t}_{2})$,
and the loading rate as
$\widehat{r}_{load}=(\widehat{\delta}_{load}-\widehat{\delta}_{0})/(\widehat{t}_{1})$
with reference to the the inset of Figure 4.a. We endeavored to investigate
the impact of different loading protocols meticulously although the unloading
curves presented in Figure 4 are restricted to the six different protocols
described in Table 1.

Table 1: Description of parameters governing the loading protocols for the curves presented in Figs. 4 and 5 . Loading protocol | $\widehat{\delta}_{0}$ | $\widehat{\delta}_{load}$ | $\widehat{r}_{load}$ | $\widehat{r}$ | $\widehat{t}_{dwell}$
---|---|---|---|---|---
# 1 | 1 | 1 | very fast | 10 | 0
# 2 | 1 | 1 | very slow | 10 | 0
# 3 | 0 | 1 | 5 | 10 | 0
# 4 | 2 | 2 | very fast | 10 | 0
# 5 | 2 | 2 | very slow | 10 | 0
# 6 | 0 | 2 | 5 | 10 | 0

Figure 4: (a) Unloading curves for $\Sigma_{0}=0.05,k=0.1$ and punch of radius
$\widehat{a}/\widehat{a}_{0}=64.85$ from a fully relaxed viscoelastic surface
with different loading protocols. (b) The identical curves displayed in (a)
are reiterated here subsequent to a horizontal axis shift equal to
$\widehat{\delta}_{P=0}$.

These loading scenarios included: (i) Unloading from a fully relaxed substrate
after a slow loading process, indicated by the black curves (2,5). (ii)
Unloading following rapid loading, causing the substrate to exhibit an elastic
response with $E(t=0)=E_{\infty}$, as denoted by the blue curves (1,4). (iii)
Unloading after indenting the substrate at a constant loading rate
$\widehat{r}_{load}=5$, represented by the red curves (3,6). It’s important to
note that while the loading phase is not shown for curves (1, 2, 4, 5), we
accounted for the pre-loading effect in our simulations. Furthermore, the
maximum indentation depth $\widehat{\delta}_{load}$ was kept at
$\widehat{\delta}_{load}=1$ for curves (1, 2, 3) (solid lines) and set to
$\widehat{\delta}_{load}=2$ for curves [4, 5, 6] (dashed lines). The punch has
$a/a_{0}=64.85$, and the viscoelastic layer’s (dimensionless) thickness is
$B=b/a=0.005$ which is equivalent to $\widehat{b}/\widehat{a_{0}}=0.3242$.

Figure 5: Normalized pull-off stress as a function versus the normalized
layer’s thickness for different loading protocols.

The primary result from Fig. 4 is that, regardless of the significant
variations in loading conditions, the magnitude of the pull-off stress remains
consistent across various loading histories. Indeed, one can conclude that the
pull-off stress remains nearly unaffected by the loading history. In Figure
4.b, we present the same curves as in Fig. 4.a, but with a shift in the
horizontal axis by $\delta_{P}$. This shift corresponds to the indentation
depth at which the normal load vanishes during unloading. It helps to better
observe the slight changes in the unloading trajectories. Furthermore a
comprehensive investigation on the effect of the different loading scenario
was conducted for various layer thicknesses. In Figure 5, we illustrate
$\widehat{\overline{\sigma}}_{po}$ versus the variation in thickness
$b/a_{0}$. The legend in Figure 5 clarifies that the plot presents results
derived from all the loading conditions in Table 1, all with the same
unloading rate. Remarkably, these plots closely overlap, indicating very
similar values across the various simulations for all the thicknesses tested.
To further quantify the distinctions between these loading cases, we examined
the relative changes in pull-off stress with respect to one specific case that
serves as the foundation for our subsequent investigations. The results are
plotted in the inset of Figure 5 as a relative change for various thickness
values. It’s evident that within a certain accuracy we can state that the
detachment force of a flat indenter from a viscoelastic adhesive strip is
negligibly influenced by the loading history of the contact. For the remainder
of the paper we will consider unloading the viscoelastic strip from a fully
relaxed condition, unless explicitly stated otherwise, we assume
$\widehat{\delta}_{0}=\widehat{\delta}^{{load}}=1$, ${\widehat{t}_{dwell}}=0$.

###  5.2 Dependence on the unloading rate

After establishing that the loading history does not influence the pull-off
stress $\widehat{\overline{\sigma}}_{po}$, we examine how
$\widehat{\overline{\sigma}}_{po}$ varies with respect to the layer thickness
for four different unloading rates: $\widehat{r}=[0.1,1,10,100]$ represented
in Fig. 6 by black diamonds, green circles, red squares, and pink triangles,
respectively. Figure 6 shows a comprehensive analyses for the punch semi-width
$\widehat{a}/\widehat{a}_{0}=64.85$. The results are obtained starting from a
fully relaxed substrate. For $b<b_{0}\simeq 1.27a_{0}$, we reach the cohesive
limit, where the pull-off stress remains independent on both the unloading
rate and the layer thickness, approaching the theoretical value
$\widehat{\overline{\sigma}}_{po}=1$. Most importantly, for $b_{0}<b<b_{1}$,
the curves align well with the LEFM (Linear Elastic Fracture Mechanics)
solution we have derived in Section 3 showing a scaling of $\propto b^{-1/2}$.
Here, the pull-off stress increases with the unloading rate, and the pull-off
data consistently stay well by the “slow” and “fast” limits we derived,
represented by the blue dashed and solid black lines, respectively. For a
thickness larger than $b_{1}$ the curves align with the half-plane solution
and, for a given unloading rate, the pull-off stress remains constant
independently on the layer thickness.

Figure 6: Normalized pull-off stress as a function of the normalized layer’s
thickness for different unloading rates of a punch with
$\widehat{a}/\widehat{a}_{0}=64.85$.

Notice that in the theoretical elastic solution the detachment happens with no
propagation ($a_{c}=a$, being $a_{c}$ the semi-width of the crack ligament).
Clearly, this condition is never achieved in a more refined cohesive zone
model, as it is the one we have implemented numerically. This accounts for the
small deviations we found in the limiting case of very fast and very slow
unloading between numerical and theoretical results. Nevertheless, to
ascertain the correctness of the numerical viscoelastic results, the plots
also include the curves obtained unloading an elastic strip with modulus
$E_{0}$ (filled blue stars) and $E_{\infty}$ (empty blue stars). One easily
recognizes that the viscoelastic solutions are perfectly bounded between the
two limiting elastic cases.

To support our conclusion we focus on the mechanism of crack propagation and
stress distribution at the interface from the unloading onset up to pull-off.
Figures 7 show the stress distribution for three specific cases out of the 120
cases shown in Fig. 6. All the cases are for unloading rate $\widehat{r}=100$.
The punch radius in Figure 7 is $\widehat{a}/\widehat{a}_{0}=64.85$. The
corresponding points for these three cases are highlighted with red squares in
Fig. 6. According to Fig. 7, during unloading, the crack propagates at the
interface hence the semi-width of the crack ligament $a_{c}$ is smaller than
the punch semi-width $a$ when pull-off happens. This explains the difference
between the expected pull-off stress from the analytical limits and the actual
pull-off stress.

Figure 7 displays three distinct cases with different values of $B$ (0.001,
0.25, and 10), denoted as Figures 7a, 7b, and 7c, respectively. Figure 7a
pertains to the cohesive zone, where the detachment occurs at $a_{c}\simeq a$
and with a uniform distribution of tensile tractions at the interface. For a
more comprehensive understanding of crack propagation and the detachment
process, we have included gap plots on the right side of Figure 7. These plots
represent the gap $H(x)={h(x)}/{h_{0}}-1$ between the rigid punch and the
viscoelastic layer as a function of the in plane coordinate $x$. Figure 7a
illustrates that in the cohesive region the gap is uniformly distributed at
the interface, with no crack formation. This is, in fact, the reason why we
can achieve $\widehat{\overline{\sigma}}_{po}$ numerically with the same
results as the expected analytical results (see Figure 6 for $B=0.001$). In
contrast, for the other two cases, as depicted in Figures 7b and 7c, we
observe crack propagation, with detachment occurring at $a_{c}<a$. Notably,
for larger values of $B$, a reduction in the ratio $a_{c}/a$ is evident,
resulting in the small deviation observed in Fig. 6 between the numerical and
the analytical results.

Figure 7: Gap function (the right plots) and the stress distribution (the left
plots) on the surface of a layer with parameters $\Sigma_{0}=0.05$, $k=0.1$,
and $\widehat{a}/\widehat{a}_{0}=64.85$ are presented for various geometries:
(a) $B=0.001$, (b) $B=0.25$, and (c) $B=10$, all under an unloading rate of
$\widehat{r}=100$. Each plot displays results for three different moments,
with pull-off data highlighted in red.

A more detailed view on the dependence of the pull-off stress on rate effects
is shown in Fig. 8 that shows $\widehat{\overline{\sigma}}_{po}$ as a function
of the crack speed at pull-off, defined as

$V_{c}=-da_{c}/dt\;,$ (29)

where $a_{c}$ represents the crack ligament semi-width, which decreases as the
crack propagates at the interface. The results are closely related to the
interaction between adhesion and viscoelastic dissipation in the strip (see
also [33]), indeed it represents one of the major objective of viscoelastic
crack propagation theories [27, 25, 26, 28, 22, 23, 24]. Consequently, we
conducted additional analyses to examine this effect. We selected four
different cases with a punch radius of $\widehat{a}/\widehat{a}_{0}=64.85$ and
the corresponding thickness ratios of
$\widehat{b}/\widehat{a}_{0}=[3.24,6.48,12.97,32.42]$, corresponding to the
blue, orange, yellow, and purple curves, respectively. We conducted numerical
experiments with 20 different unloading rates ranging from $\widehat{r}=0.1$
to $\widehat{r}=100$ to obtain curves representing a wide range of the
dimensionless crack velocity $\widehat{V}_{c}=V_{c}\tau/h_{0}$ at pull-off.

Figure 8: Normalized pull-off stress as a function of the normalized crack
velocity at pull-off for four different values of $b/a_{0}$ for a punch with
$\widehat{a}/\widehat{a}_{0}=64.85$ unloaded from a fully relaxed viscoelastic
surface at different unloading rates.

The analysis of Fig. 8 illustrates clearly the trend: thin layers and high
retraction velocity favour high pull-off stress. Nevertheless, this effect is
mitigated when the $b\approx a_{0}\approx b_{0}$ as, in the cohesive region,
the detachment tends to happen at a uniform stress.

Figure 9: Normalized effective surface energy as a function of the crack
velocity at pull-off for four different values of $b/a_{0}$ for a punch with
$\widehat{a}/\widehat{a}_{0}=64.85$ unloaded from a fully relaxed viscoelastic
surface at different unloading rates. Dashed lies are obtained with Eq. (36).
Maximum enhancement for the dimensionless effective surface energy as a
function of $b/a_{0}$. Inset: maximum enhancement of the effective surface
energy obtained numerically at very high unloading rates for
$a/a_{0}=\left[32.42,64.84,129.69\right]$, respectively red circles, green
stars, blue triangles. The solid black line shows the prediction of Eq. (32).

Although, the pull-off force reduces increasing the layer thickness the
enhancement in terms of effective surface energy remains the same when moving
from low to high unloading rates, provided that $b>b_{0}$. Based on the thin
layer elastic solution Eq. (10) we define the effective surface energy as

${\Delta\gamma}_{eff}=\frac{{\overline{\sigma}}_{po}^{2}b}{2E^{*}_{0}}\left(\frac{a}{a_{c}}\right)^{2}\;,$
(30)

where we considered that in general the detachment happens at $a_{c}<a$.
Hence, in dimensionless form,

$\widehat{\Delta\gamma}_{eff}=\frac{\Delta\gamma_{eff}}{\Delta\gamma_{0}}=\frac{\widehat{\overline{\sigma}}_{po}^{2}\Sigma_{0}\widehat{b}}{2\alpha_{LJ}}\left(\frac{\widehat{a}}{\widehat{a}_{c}}\right)^{2}\;.$
(31)

Fig. 9 shows that the normalized effective surface energy
$\widehat{\Delta\gamma}_{eff}$ increases monotonically with respect to the
crack velocity $\widehat{V}_{c}$ at pull-off up to a certain plateau value. In
the case of $b/a_{0}=32.42$, the normalized effective surface energy reaches
its theoretical upper limit $\widehat{\Delta\gamma}_{eff}=1/k$ (all our
simulations are for k=0.1). Notice that for any thickness of the layer larger
than $b_{0\infty}$ one would get the maximum possible enhancement $1/k$. Care
should be taken when interpreting the data using Eq. 31 as, if the latter is
used for $b>b_{1}$ this may lead to unrealistic enhancements
$\widehat{\Delta\gamma}_{eff}$, which is due to the fact that for $b>b_{1}$
the halfplane solution should be considered.

It’s important to note that, due to the finite size effect, for cases with
$b<b_{0\infty}$, we observe the maximum enhancement of the normalized
effective surface energy to be lower than $1/k$. This happens because, for the
very thin layer, the cohesive region is approached, namely the DMT-type
failure rather than the JKR-type, in the Peng et al. [13] terminology. In the
latter case, if we assume $\widehat{\overline{\sigma}}_{po}=1$, and
$\widehat{a}=\widehat{a}_{c}$, and with the acquisition of Eq. (6), one can
obtain the following relation for the maximum enhancement that can be reached
at high retraction rates

$\widehat{\Delta\gamma}_{eff}|_{max}=\left(\frac{\pi}{4}\right)\left(\frac{b}{a_{0}}\right),$
(32)

which turns out to be solely dependent on the ratio $b/a_{0}$. In order to
validate our upper bound enhancement factor (Eq. (32)), we considered three
distinct values of punch semi-width $a/a_{0}$, on fully relaxed viscoelastic
substrate with varying $b/a_{0}$ ratios, unloaded at a high unloading rate
$\widehat{r}=100$. The inset of Fig. 9 shows that the maximum enhancement
obtained numerically compared very well with Eq. (32).

We incorporated this correction in Greenwood (2004) theory [21] for crack
propagation in viscoelastic semi-infinite media, which, in its original form
gives

$\widehat{\Delta\gamma}_{eff}=\left[k+\left(1-k\right)\frac{\alpha}{2}\int_{0}^{1}H\left(\xi\right)\exp\left(-\alpha\left(1-\xi\right)\right)d\xi\right]^{-1}\;,$
(33)

where

$\displaystyle H\left(\xi\right)$
$\displaystyle=2\xi^{1/2}-\left(1-\xi\right)\ln\left(\frac{1+\xi^{1/2}}{1-\xi^{1/2}}\right)\;,$
(34) $\displaystyle\alpha$
$\displaystyle=\frac{\pi}{4\Sigma_{0}}\frac{\widehat{\Delta\gamma}_{eff}\alpha_{LJ}}{\widehat{V}_{c}}\;.$
(35)

Equation (33) for very slow propagation gives
$\widehat{\Delta\gamma}_{eff}=1$, while at high speed provides the maximum
enhancement $\left.\widehat{\Delta\gamma}_{eff}\right|_{\max}=1/k$. This
picture, on which all present theories agree, is valid for semi-infinite
solids, nevertheless, in agreement with recent results [30, 31, 29], we have
found that due to finite size effects the maximum enhancement may be
consistently reduced. For the present problem, if $b<b_{0\infty}$, the maximum
enhancement will be given by
$\left.\widehat{\Delta\gamma}_{eff}\right|_{\max}=\left(\pi/4\right)\left(b/a_{0}\right)$,
so we propose here a generalization of Eq. (33) for $b_{0}<b<b_{0\infty}$

$\widehat{\Delta\gamma}_{eff}\left(\widehat{V}_{c},\frac{b}{a_{0}}\right)=\left[\frac{4}{\pi\left(b/a_{0}\right)}+\left(1-\frac{4}{\pi\left(b/a_{0}\right)}\right)\frac{\alpha}{2}\int_{0}^{1}H\left(\xi\right)\exp\left(-\alpha\left(1-\xi\right)\right)d\xi\right]^{-1}\;,$
(36)

where we have explicitly indicated that now the velocity dependent effective
surface energy depends not only on the crack speed, but also on the ratio
between the layer thickness and the fracture length $a_{0}$. Notice that for
$b>b_{0\infty}$ Eq. (33) remains valid, while for $b<b_{0}$ the effective
energy is velocity independent and equal to
$\left.\widehat{\Delta\gamma}_{eff}\right|_{\max}=\left(\pi/4\right)\left(b/a_{0}\right)$.
Figure 9 compares the predictions obtained with the finite size Greenwood
model (Eq. (36)) against the numerical results, which are found in fairly good
agreement. This result (Eq. (32)), obtained with the more general cohesive-
based theory, could also be used to correct the theory of Persson-Brener
(2005) [22] and Persson (2017) [23], as we shall do in Appendix-II.

One important parameter to examine is the work of separation, which is defined
as

$\widehat{w}_{sep}=\frac{{w}_{sep}}{2aL\Delta\gamma_{0}}=\int_{\widehat{\delta}_{P}}^{\infty}\widehat{\overline{\sigma}}\,d\widehat{\delta}\;.$
(37)

which indicates the energy spent during the unloading phase to separate the
contact. We calculated this parameter for four different layer thicknesses,
specifically $b/a_{0}=[3.24,6.48,12.97,32.42]$, which corresponds to the blue,
orange, yellow, and purple curves in Fig. 10. Similarly to previous research
works [30, 31, 29] the $\widehat{w}_{sep}$ has a typical bell shape; at low
and high velocity there is little energy expenditure to separate the contact
as the material behaves essentially as elastic, but for intermediate regimes
$\widehat{w}_{sep}$ presents a maximum related to the dissipative phenomena
happening in the viscoelastic layer.

Figure 10: Normalized work of separation as a function of the crack velocity
at pull-off for four different values of $b/a_{0}$ for a punch with
$\widehat{a}/\widehat{a}_{0}=64.85$ unloaded from a fully relaxed viscoelastic
surface at different unloading rates.

##  6 Conclusions

The plane problem of the detachment of a large flat punch from an adhesive
viscoelastic layer of finite thickness $b$ has been studied. First, we have
derived an elastic model based on the “thin strip” assumption by Johnson [38].
It was found that the pull-off stress decays as $\propto 1/\sqrt{b}$.
Nevertheless, this functional dependence is bounded (i) for very thin layer by
the cohesive limit where the pull-off stress equals the theoretical stress of
the material, (ii) for very thick layer by the halfplane limiting solution.
The elastic model provided the bounds for the viscoelastic analysis. We found
that if the layer is thin, particularly at high enough retraction velocity,
the theoretical limit of the material could be reached. This turns
particularly interesting as for soft polymers $E_{\infty}/E_{0}$ may easily be
of the order of $10^{3}\div 10^{4}$ and this amplifies the layer thickness for
which the theoretical strength can be observed. Clearly, this behavior will be
hindered by the fact that during unloading the crack starts to propagate at
the interface hence, at pull-off, the actual crack ligament $a_{c}$ is smaller
than the punch semi-width.

Theoretical predictions have been compared with boundary element numerical
simulations for a standard linear viscoelastic material and using a Lennard-
Jones force-interaction law. We have shown that the loading conditions have a
negligible effect on the pull-off force, in contrast with what was shown for a
Hertzian geometry. Instead, the pull-off force consistently increases with the
unloading rate up to a certain plateau given by the cohesive strength of the
interface.

Finally, we have shown that when the data are represented in terms of
effective surface energy, at high velocity the theoretical enhancement given
by $E_{\infty}/E_{0}$ is reached only when the layer thickness is larger than
a characteristic lengthscale $b_{0\infty}$. For $b_{0}<b<b_{0\infty}$ the
maximum adhesion enhancement is limited by finite size effect and in
particular we found
$\widehat{\Delta\gamma}_{eff}|_{\max}=\left(\pi/4\right)\left(b/a_{0}\right)$.
Hence, we have proposed an extension of Greenwood and Persson crack
propagation theories accounting for finite size effects which we found in good
agreement with numerical results.

##  7 Funding

A.P. and A.M. were supported by the European Union (ERC-2021-STG, “Towards
Future Interfaces With Tuneable Adhesion By Dynamic Excitation” - SURFACE,
Project ID: 101039198, Grant No. CUP: D95F22000430006). Views and opinions
expressed are however those of the authors only and do not necessarily reflect
those of the European Union or the European Research Council. Neither the
European Union nor the granting authority can be held responsible for them.
A.P. was partially supported by Regione Puglia (Italy), project ENOVIM (Grant
No. CUP: D95F21000910002) granted within the call “Progetti di ricerca
scientifica innovativi di elevato standard internazionale” (art. 22 della
legge regionale 30 novembre 2019, n. 52 approvata con A.D. n. 89 of
10-02-2021, BURP n. 25 del 18-02-2021). All the authors acknowledge support by
the Italian Ministry of University and Research (MUR) under the programme
“Department of Excellence” Legge 232/2016 (Grant No. CUP—D93C23000100001).

## Appendices

## Appendix-I. Perfectly bonded layer

For the case of a layer perfectly bonded to the rigid substrate, following
Johnson [38], we only need to correct the previous results for the Poisson
effect with $\zeta=\frac{\left(1-\nu\right)^{2}}{1-2\nu}$ (this requires
$\nu\lesssim 0.45$, [38]), so that we have

$\delta_{po}=\sqrt{\frac{2b\Delta\gamma}{\zeta
E^{\ast}}},\qquad\overline{\sigma}_{po}=\sqrt{\frac{2\zeta
E^{\ast}\Delta\gamma}{b}}\;,$ (38)

where $\overline{\sigma}_{po}$ equals the cohesive strength of the material
for the layer thickness

$b_{0}=\frac{2\zeta E^{\ast}\Delta\gamma}{\sigma_{0}^{2}}=\frac{4}{\pi}\zeta
a_{0}\;.$ (39)

## Appendix-II. An extension of Persson and Brener viscoelastic crack
propagation theory accounting for finite size effects

Using Persson and Brener theory [22] for a single relaxation time material
gives an implicit equation for the effective adhesion energy [32]

$\widehat{\Delta\gamma}_{eff}\left(\widehat{V}_{c}\right)=\left[1-\left(1-k\right)\frac{\widehat{\Delta\gamma}_{eff}}{\beta\widehat{V}_{c}}\left(\sqrt{1+\left(\frac{\beta\widehat{V}_{c}}{\widehat{\Delta\gamma}_{eff}}\right)^{2}}-1\right)\right]^{-1}\;,$
(40)

being $\beta=64\Sigma_{0}/\left(9\sqrt{3}\right)$. Eq. (40) can be extended to
finite size systems by accounting that for a very thin layer the maximum
enhancement will be reduced to
$\left.\widehat{\Delta\gamma}_{eff}\right|_{\max}=\left(\pi/4\right)\left(b/a_{0}\right)$,
so we propose here a generalization of Eq. (40) in order to take into account
finite size systems, i.e. for $b_{0}<b<b_{0\infty}$

$\widehat{\Delta\gamma}_{eff}\left(\widehat{V}_{c},\frac{b}{a_{0}}\right)=\left[1-\left(1-\frac{4}{\pi\left(b/a_{0}\right)}\right)\frac{\widehat{\Delta\gamma}_{eff}}{\beta\widehat{V}_{c}}\left(\sqrt{1+\left(\frac{\beta\widehat{V}_{c}}{\widehat{\Delta\gamma}_{eff}}\right)^{2}}-1\right)\right]^{-1}\;,$
(41)

where we have explicitly indicated that the normalized effective surface
energy depends not only on the crack speed, but also on the ratio between the
layer thickness and the fracture length $a_{0}$. A comparison between the
numerical results and Eq. (41) is shown in Fig. 11.

Figure 11: Normalized effective surface energy as a function of the crack
velocity at pull-off for four different values of $b/a_{0}$ for a punch with
$\widehat{a}/\widehat{a}_{0}=64.85$ unloaded from a fully relaxed viscoelastic
layer ($k=0.1$). Dashed lines are obtained with Eq. (41).

## References

  * [1] B. Lorenz, Y. R. Oh, S. K. Nam, S. H. Jeon, B. N. J. Persson, Rubber friction on road surfaces: Experiment and theory for low sliding speeds, The Journal of Chemical Physics 142 (19) (2015).
  * [2] B. Mazzolai, A. Mondini, F. Tramacere, G. Riccomi, A. Sadeghi, G. Giordano, S. Carminati, Octopus-inspired soft arm with suction cups for enhanced grasping tasks in confined environments, Advanced Intelligent Systems 1 (6) (2019) 1900041.
  * [3] G. Giordano, M. Carlotti, B. Mazzolai, A perspective on cephalopods mimicry and bioinspired technologies toward proprioceptive autonomous soft robots, Advanced Materials Technologies 6 (12) (2021) 2100437.
  * [4] G. Giordano, M. Gagliardi, Y. Huan, M. Carlotti, A. Mariani, A. Menciassi, B. Mazzolai, Toward mechanochromic soft material-based visual feedback for electronics-free surgical effectors, Advanced Science 8 (15) (2021) 2100418.
  * [5] J. Shintake, V. Cacucciolo, D. Floreano, H. Shea, Soft robotic grippers, Advanced Materials 30 (29) (2018) 1707035.
  * [6] C. Huon, A. Tiwari, C. Rotella, P. Mangiagalli, B. N. J. Persson, Air, helium and water leakage in rubber o-ring seals with application to syringes, Tribology Letters 70 (2) (2022) 35.
  * [7] F. Forsbach, M. Heß, A. Papangelo, A two-scale fem-bam approach for fingerpad friction under electroadhesion, Frontiers in Mechanical Engineering 8 (2023) 1074393.
  * [8] L. Felicetti, C. Sutter, E. Chatelet, A. Latour, L. Mouchnino, F. Massi, Tactile discrimination of real and simulated isotropic textures by friction-induced vibrations, Tribology International 184 (2023) 108443.
  * [9] S. A. Burke, M. Ritter-Jones, B. P. Lee, P. B. Messersmith, Thermal gelation and tissue adhesion of biomimetic hydrogels, Biomedical materials 2 (4) (2007) 203.
  * [10] B. M. Gumbiner, Cell adhesion: the molecular basis of tissue architecture and morphogenesis, Cell 84 (3) (1996) 345–357.
  * [11] A. Asbeck, S. Dastoor, A. Parness, L. Fullerton, N. Esparza, D. Soto, M. Cutkosky, Climbing rough vertical surfaces with hierarchical directional adhesion (2009) 2675–2680.
  * [12] C. A. Dahlquist, Pressure-sensitive adhesives (1969) 219–260.
  * [13] B. Peng, X.-Q. Feng, Q. Li, Decohesion of a rigid flat punch from an elastic layer of finite thickness, Journal of the Mechanics and Physics of Solids 139 (2020) 103937.
  * [14] R. Christensen, Theory of viscoelasticity: an introduction, Elsevier, 2012.
  * [15] C. Creton, M. Ciccotti, Fracture and adhesion of soft materials: a review, Reports on Progress in Physics 79 (4) (2016) 046601.
  * [16] M. Barquins, D. Maugis, Tackiness of elastomers, The Journal of Adhesion 13 (1) (1981) 53–65.
  * [17] A. N. Gent, Adhesion and strength of viscoelastic solids. is there a relationship between adhesion and bulk properties?, Langmuir 12 (19) (1996) 4492–4496.
  * [18] A. N. Gent, R. P. Petrich, Adhesion of viscoelastic materials to rigid substrates, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 310 (1502) (1969) 433–448.
  * [19] A. N. Gent, J. Schultz, Effect of wetting liquids on the strength of adhesion of viscoelastic material, Journal of Adhesion 3 (4) (1972) 281–294.
  * [20] M. L. Williams, R. F. Landel, J. D. Ferry, The temperature dependence of relaxation mechanisms in amorphous polymers and other glass-forming liquids, Journal of the American Chemical Society 77 (1955) 3701–3707.
  * [21] J. Greenwood, The theory of viscoelastic crack propagation and healing, Journal of Physics D: Applied Physics 37 (18) (2004) 2557.
  * [22] B. N. J. Persson, E. A. Brener, Crack propagation in viscoelastic solids, Physical Review E 71 (3) (2005) 036123.
  * [23] B. N. J. Persson, Crack propagation in finite-sized viscoelastic solids with application to adhesion, EPL (Europhysics Letters) 119 (1) (2017) 18002.
  * [24] B. N. J. Persson, A simple model for viscoelastic crack propagation, The European Physical Journal E 44 (2021) 1–10.
  * [25] R. A. Schapery, A theory of crack initiation and growth in viscoelastic media, Int. J. Fracture 11 (1975) 141–159.
  * [26] R. A. Schapery, A theory of crack initiation and growth in viscoelastic media ii. approximate methods of analysis, Int. J. Fracture 11 (3) (1975) 369–388.
  * [27] G. J. A, J. K. L, The mechanics of adhesion of viscoelastic solids, Phil. Mag. A 43 (1981) 697–711.
  * [28] R. Schapery, A theory of viscoelastic crack growth: revisited, International Journal of Fracture 233 (2022) 1–16.  doi:10.1007/s10704-021-00605-z.
  * [29] A. Papangelo, M. Ciavarella, Detachment of a rigid flat punch from a viscoelastic material, Tribology Letters 71 (2) (2023) 48.
  * [30] L. Afferrante, G. Violano, On the effective surface energy in viscoelastic hertzian contacts, Journal of the Mechanics and Physics of Solids 158 (2022) 104669.
  * [31] G. Violano, L. Afferrante, On the long and short-range adhesive interactions in viscoelastic contacts, Tribology Letters 70 (3) (2022) 1–5.
  * [32] M. Ciavarella, G. Cricrì, R. McMeeking, A comparison of crack propagation theories in viscoelastic materials, Theoretical and Applied Fracture Mechanics 116 (2021) 103113.
  * [33] M. Ciavarella, A. Papangelo, R. McMeeking, Crack propagation at the interface between viscoelastic and elastic materials, Engineering Fracture Mechanics 257 (2021) 108009.
  * [34] M. Ciavarella, A. Papangelo, R. M. McMeeking, Transient and steady state viscoelastic crack propagation in a double cantilever beam specimen, International journal of mechanical sciences 229 (2022).
  * [35] C. Müller, M. Samri, R. Hensel, E. Arzt, M. H. Müser, Revealing the coaction of viscous and multistability hysteresis in an adhesive, nominally flat punch: A combined numerical and experimental study, Journal of the Mechanics and Physics of Solids 174 (2023) 105260.
  * [36] L. Afferrante, G. Violano, G. Carbone, Exploring the dynamics of viscoelastic adhesion in rough line contacts, Scientific Reports 13 (1) (2023) 15060.
  * [37] F. Forsbach, A simple semi-analytical method for solving axisymmetric contact problems involving bonded and unbonded layers of arbitrary thickness, Machines 11 (4) (2023) 474.
  * [38] K. L. Johnson, Contact Mechanics, Cambridge University Press, New York, 1985.
  * [39] J. R. Barber, Contact Mechanics, Springer International Publishing, Berlin, 2018.
  * [40] D. Maugis, Contact, adhesion and rupture of elastic solids, Vol. 130, Springer Science & Business Media, Heidelberg, 2000.
  * [41] A. Jagota, S. J. Bennison, C. A. Smith, Int. j. fracture 104, 105–130 (2000), International Journal of Fracture 104 (2000) 105–130.
  * [42] T. Tang, A. Jagota, M. K. Chaudhury, C. Y. Hui, Thermal fluctuations limit the adhesive strength of compliant solids, The Journal of Adhesion 82 (7) (2006) 671–696.
  * [43] A. Bonfanti, J. L. Kaplan, G. Charras, A. Kabla, Fractional viscoelastic models for power-law materials, Soft Matter 16 (26) (2020) 6002–6020.
  * [44] A. Papangelo, M. Ciavarella, A numerical study on roughness-induced adhesion enhancement in a sphere with an axisymmetric sinusoidal waviness using lennard-jones interaction law, Lubricants 8 (9) (2020) 90.
  * [45] R. H. Bentall, K. L. Johnson, An elastic strip in plane rolling contact, International Journal of Mechanical Sciences 10 (8) (1968) 637–663.
  * [46] R. Christensen, Theory of viscoelasticity: an introduction, Elsevier, New York, 2012.