Dmitry Ponomarev

A GENERALISED TIME-EVOLUTION MODEL FOR CONTACT PROBLEMS WITH WEAR AND ITS ANALYSIS

In this paper, we revisit some classical and recent works on modelling sliding contact with wear and propose their generalisation. Namely, we upgrade the relation between the pressure and the wear rate by incorporating some non-local time-dependence. To this effect, we use a combination of fractional calculus and relaxation effects. Moreover, we consider a possibility when the load is not constant in time. The proposed model is analysed and solved. The results are illustrated numerically and comparison with similar models is discussed.

Introduction

Wear of material is a complicated process which involves different phenomena such as abrasion, adhesion and crack formation. Wear processes have been studied for decades and numerous empirical models have been proposed in attempt to fit experimental data for particular settings. One classical model of wear is due to Archard [START_REF] Archard | Contact and Rubbing of Flat Surfaces[END_REF]. According to this model, the wear rate is proportional to the load with a power-law dependence. In many settings, this reduces to a linear relation between the wear rate and the pressure (see e.g. [START_REF] Popov | Contact Mechanics and Friction[END_REF]) leading to several effectively solvable tribological models for sliding of an indented punch (stamp), see e.g. [START_REF] Argatov | Effective wear coefficient and wearing-in period for a functionally graded wear-resisting punch[END_REF][START_REF] Argatov | Wear contact problem with friction: Steady-state regime and wearing-in period[END_REF][START_REF] Argatov | A Macro-Scale Approximation for the Running-In Period[END_REF][START_REF] Feppon | Efficient steady-state computation for wear of multimaterial composites[END_REF][START_REF] Kovalenko | On an efficient method of solving contact problems with linearly deformable base with a reinforcing coating (in Russian)[END_REF][START_REF] Zhu | Simulation of Sliding Wear in Mixed Lubrication[END_REF]. It seems natural to explore a larger class of such models based on a more general relation between the wear rate and the contact pressure. A qualitative upgrade of this relation could be achieved by accounting for non-local (temporal) dependencies. Investigation of how a new wear model affects the solution of a problem would give a hope to extend the applicability of original models to contexts with different wear mechanisms and a broader set of materials.

Here, we consider the following instance of the two-dimensional punch problem: a rigid wearable punch, subject to a given normal (vertical) load, slides on a thick elastic layer or a half-plane, with a prescribed speed. The sliding speed is taken to be constant, but the normal load may be time-dependent. The contact area is assumed to be fixed (which is normally expected if the load is sufficiently large, see e.g. [START_REF] Komogortsev | Contact between a moving stamp and an elastic half-plane when there is wear[END_REF]).

Since the layer is homogeneous and material wear occurs on the interface, the problem can be effectively described by one-dimensional integral equations. On this level, the presence of wear in the problem manifests itself as an additional term in the integral equation for the pressure. This term stems from the linear Archard's wear law and brings a temporal dependence to the problem.

In the present work, we upgrade the wear term so that it corresponds to a more general differential relation between the wear rate and the contact pressure. This more general relation is meant to have two features (and combinations thereof). First, it incorporates a wear-relaxation effect which is consistent with typical observations that wear is the most intense in the beginning of the sliding process. Second, the time-derivative in this relation may be changed to a fractional order (in a Riemann-Liouville sense) which adds another non-locality and is motivated by recent success of fractional calculus in mechanics (viscoelasticity) and other applied contexts [START_REF] Gorenflo | Mittag-Leffler Functions, Related Topics and Applications[END_REF]. Importance of non-local relations between wear and contact pressure is also stressed in recent work [START_REF] Argatov | Contact Geometry Adaptation in Fretting Wear: A Constructive Review[END_REF].

As we shall see, newly introduced parameters can affect qualitative behavior of the solution. Namely, depending on a choice of the parameter values, the solution may exhibit exponential or algebraic decay in time that would range between monotone and arbitrary oscillatory. Therefore, such a generalised model is highly desirable and is expected to be useful for fitting experimentally observable data. On the other hand, this model is almost explicitly solvable, meaning that the solution can be written up in a closed form in terms of some auxiliary functions which could be precomputed (numerically or asymptotically). In particular, it is rather straightforward to analyse long-time behavior of this model, namely, its convergence speed to the stationary pressure distribution. Moreover, computing such a pressure profile itself is a problem of an essential practical interest.

The outline of the paper is as follows. We formulate a new model in Section 2. In Section 3, we derive the proposed form of the wear relation and recast the model in form of a single integral equation. Then, in Section 4, Geometry of the problem. A balance of displacements under the punch is governed by the equation: ηp (x, t) + ´a -a K (x -ξ) p (ξ, t) dξ = δ (t) -∆ (x) -w (x, t), x ∈ (-a, a), t ≥ 0.

we deal with the solution of the general model and its analysis. In particular, in Subsection 4.2, under appropriate conditions, we deduce the existence of a unique solution and provide its explicit form in terms of spectral functions of a pertinent integral operator. This is detailed even more for a concrete choice of the kernel function in Subsection 4.3. Furthermore, in Section 5, we show that the solution of the proposed model has the anticipated behavior consistent with previous models. A collection of mathematical results needed for Sections 4-5 is outsourced to Appendix. Section 6 is dedicated to the numerical solution of the proposed model: its validation and exploration. Finally, we conclude with Section 7 where we summarise and discuss the results outlining potential further work.

The model

According to [START_REF] Aleksandrov | Mathematical methods in problems with wear[END_REF][START_REF] Argatov | Effective wear coefficient and wearing-in period for a functionally graded wear-resisting punch[END_REF][START_REF] Argatov | Wear contact problem with friction: Steady-state regime and wearing-in period[END_REF][START_REF] Argatov | A Macro-Scale Approximation for the Running-In Period[END_REF][START_REF] Galin | Contact problems of the theory of elasticity in the presence of wear[END_REF][START_REF] Komogortsev | Contact between a moving stamp and an elastic half-plane when there is wear[END_REF][START_REF] Vorovich | Non-classical mixed contact problems in elasticity theory[END_REF], a balance of displacements at the contact interface (see Figure 2.1) yields an integral equation which involves two unknowns: the contact pressure p (x, t) and the indentation δ (t). This equation can be formulated as follows

(2.1) η (p (x, t) -p (x, 0)) + ˆa -a K (x -ξ) (p (ξ, t) -p (ξ, 0)) dξ + w [p] (x, t) = δ (t) -δ (0) , x ∈ (-a, a) , t > 0.
Furthermore, (2.1) should be complemented by the equilibrium condition

(2.2) ˆa -a p (x, t) dx = P (t) , t ≥ 0,
and the integral equation for the initial pressure distribution p (x, 0)

(2.3) ηp (x, 0) + ˆa -a K (x -ξ) p (ξ, 0) dξ = δ (0) -∆ (x) , x ∈ (-a, a) .
Here, K (x) is a given function which links the contact pressure to the vertical displacement, η ≥ 0 is a given constant, (-a, a) is the contact interval, P (t) > 0 is the contact load, w [p] (x, t) is the wear term which is a given mapping of the contact pressure p (x, t) and the punch profile ∆ (x) is a known function.

The main feature of our model is a new form of the wear term w [p] (x, t) entering (2.1). This wear term is given by

(2.4) w [p] (x, t) = -νµ 1/α-1 ˆt 0 E α µ 1/α (t -τ ) p (x, τ ) dτ,
where E α is an auxiliary function defined by (A.7) and µ > 0, α ∈ (0, 2) are given constants.

Relation to other models. The considered model is an amalgam as it attempts to treat different settings at once as we shall discuss here. First of all, we perform most of the analysis for a rather general form of the kernel function K (x) entering (2.1) and (2.3) while only making some assumptions on it. This allows dealing with an elastic half-plane or an elastic layer subject to different boundary conditions at its foundation, see e.g. [START_REF] Alblas | Contact problems of a rectangular block on an elastic layer of finite thickness: the thin layer[END_REF][START_REF] Alblas | Contact problems of a rectangular block on an elastic layer of finite thickness: the thick layer[END_REF][START_REF] Vorovich | Non-classical mixed contact problems in elasticity theory[END_REF]Par. 11]. In Subsection 4.3, we also provide a more detailed study for a concrete kernel function pertinent to the half-plane geometry.

Works [START_REF] Argatov | A Macro-Scale Approximation for the Running-In Period[END_REF][START_REF] Argatov | Effective wear coefficient and wearing-in period for a functionally graded wear-resisting punch[END_REF][START_REF] Argatov | Wear contact problem with friction: Steady-state regime and wearing-in period[END_REF][START_REF] Komogortsev | Contact between a moving stamp and an elastic half-plane when there is wear[END_REF] involve counterparts of (2.1) and (2.3) with η = 0. The case η > 0 accounts for the surface roughness effect or coating [START_REF] Aleksandrov | On the theory of contact problems in the presence of nonlinear wear[END_REF][START_REF] Kovalenko | On an efficient method of solving contact problems with linearly deformable base with a reinforcing coating (in Russian)[END_REF] which is due to an additional deformation proportional to the local pressure [START_REF] Shtaerman | Contact problem of the theory of elasticity[END_REF]Ch. 2 Par. 8]. We note that some aspects of the solution construction and its analysis differ depending on whether η = 0 or η > 0, and this will be reflected in our formulations of theorems and propositions given in Sections 4-5.

In most settings, P (t) ≡ P 0 was taken in (2.2) for some constant P 0 > 0. A possibility of a non-constant load is mentioned in [START_REF] Argatov | Effective wear coefficient and wearing-in period for a functionally graded wear-resisting punch[END_REF] but not studied. We study the more general load since, on the one hand, it appears to be of a direct physical interest and, on the other hand, the presented solution of the model can be adapted for this case, even though some computational details become more cumbersome. Despite this generality, we perform a more concrete analysis (see Section 5) for two particular type of loads: the most important case of a constant load and the case of a transitional load (when the load may vary in time but is eventually constant).

We shall now briefly discuss some relevant works that would motivate our new model based on a peculiar choice of the wear term. As discussed in Section 1, a simple but powerful version of the wear term is derived from the linear Archard's law. According to this relation, the wear is linearly proportional to the total sliding distance and the contact pressure. Consequently, we have

(2.5) ∂ t w (x, t) = νp (x, t)
for some constant ν > 0 related to the material and the sliding speed. Hence, since at the initial moment t = 0 the worn material is absent, i.e. w (x, 0) ≡ 0 for x ∈ (-a, a), (2.5) is equivalent to

(2.6) w [p] (x, t) ≡ w (x, t) = ν ˆt 0 p (x, τ ) dτ.
Recently, in [START_REF] Argatov | Effective wear coefficient and wearing-in period for a functionally graded wear-resisting punch[END_REF], the authors considered a generalised version of (2.6), namely,

w [p] (x, t) = k 0 (x) ˆt 0 p (x, τ ) dτ
for some suitable function k 0 (x). Another generalisation was proposed in [START_REF] Argatov | Artificial neural network modeling of sliding wear[END_REF] which, in its simple form, corresponds to

w [p] (x, t) = ˆt 0 k 1 (τ ) p (x, τ ) dτ,
but can also include additional dependencies of k 1 . Numerous models with nonlinear wear have been considered, in this case we have (see e.g. [START_REF] Zhu | Simulation of Sliding Wear in Mixed Lubrication[END_REF] and [21, eq. (6.82)])

w [p] (x, t) = ν ˆt 0 p γ (x, τ ) dτ,
with some nonlinearity parameter γ > 0.

Finally, we mention that, in addition to the wear term, the friction effects were also considered in [START_REF] Argatov | Wear contact problem with friction: Steady-state regime and wearing-in period[END_REF], by incorporating in (2.1) another term proportional to ´x -a p (ξ, t) dξ.

Derivation of the model

While the generalisation of the previous setting by including an extra term proportional to the solution and allowing load to vary in time is natural, choosing a complicated looking wear term is a nontrivial feature of the present model that we shall now discuss.

3.1.

Towards the wear term (2.4). We consider a generalisation of (2.5)-(2.6) that can be obtained in two independent steps. First, we add to the right-hand side of (2.5) a term -µw (x, t) with some constant µ ≥ 0. This would result in (2.6) being replaced by

(3.1) w [p] (x, t) = ν ˆt 0 e -µ(t-τ ) p (x, τ ) dτ,
where we took into account the condition w (x, 0) = 0 representing the initial absence of the worn material. We note that models with such a hereditary wear term has already been considered, see e.g. [START_REF] Yevtushenko | The applicability of a hereditary model of wear with an exponential kernel in the onedimensional contact problem taking frictional heat generation into account[END_REF].

At the second step, we replace the differentiation in time in the left-hand side of (2.5) with a more general notion of the derivative. In particular, we use the fractional Riemann-Liouville derivative of order α ∈ (0, 1] (for the values α > 1 see remark at the end of the paragraph) that we denote D α t , namely,

(3.2) D α t w (x, t) := 1 Γ (1 -α) ∂ ∂t ˆt 0 w (x, τ ) (t -τ ) α dτ,
where Γ denotes Gamma function defined in Appendix. Note that, due to the condition w (x, 0) = 0, this derivative would also coincide with the fractional Caputo derivative (see e.g. [20, eq. (1.19)]).

In other words, we choose to generalise (2.5) to the following non-homogeneous fractional order ODE (3.3) D α t w (x, t) = -µw (x, t) + νp (x, t) . We are now going to obtain the integral form of (3.3) and hence the corresponding relation w [p] (x, t). To this end, we set t := µ 1/α t, w x, t := w x, t/µ 1/α , p x, t := p x, t/µ 1/α . Then, relation (3.3) rewrites as

(3.4) D α t w x, t = D α t w x, t -w (x, 0) = -w x, t + ν µ p x, t ,
where we used that w (x, 0) = w (x, 0) = 0. A fractional differential equation in the form (3.4) is solved in [20, eq. (3.13)] using Laplace transforms. Taking into account that w (x, 0) = 0 and α ∈ (0, 1], the obtained solution is

w x, t = - ν µ ˆt 0 p x, t -τ E α (τ ) dτ .
Therefore, 3.2. Consistency. It is noteworthy that when µ = 0 in (3.3), we take the limit µ → 0 in (2.4). In doing so, we use the asymptotics (A.14) to obtain

(3.5) w [p] (x, t) = w x, µ 1/α t = -νµ 1/α-1 ˆt 0 p (x, t -τ ) E α µ 1/α τ
(3.6) w [p] (x, t) = αν Γ (1 + α) ˆt 0 p (x, t -τ ) τ α-1 dτ = ν Γ (α) ˆt 0 p (x, τ ) (t -τ ) 1-α dτ.
The right-hand side of (3.6) contains the Riemann-Liouville integral of order α which is consistent with the inversion of the respective fractional derivative appearing in (3.3) (see [20, eqs. (1.29), (1.2)]).

On the other hand, by taking α = 1 in (2.4), we have E 1 (z) = -e -z , and hence we recover the non-fractional order model corresponding to (3.1).

Hence, our present model can be succinctly characterised as a model with the wear term interpolating between the exponential relaxation and a new fractional order model.

3.3.

The integral equation. Equations (2.1)-(2.3) with the wear term (2.4) can be effectively reduced to a single equation. To achieve that, we integrate (2.1) in the x variable over (-a, a) and use (2.2) and (2.4). This yields (3.7)

η (P (t) -P (0))+ ˆa -a K 1 (ξ) (p (ξ, t) -p (ξ, 0)) dξ-νµ 1/α-1 ˆt 0 E α µ 1/α (t -τ ) P (τ ) dτ = 2a (δ (t) -δ (0)) , t > 0,
where we introduced the function K 1 (ξ) := ´a -a K (ζ -ξ) dζ. Dividing (3.7) over 2a and subtracting the resulting equation from (2.1), we eliminate δ (t) -δ (0) and obtain, for x ∈ (-a, a), t > 0, η p (x, t) -p (x, 0) -

P (t) -P (0) 2a + ˆa -a K (x -ξ) - 1 2a K 1 (ξ) (p (ξ, t) -p (ξ, 0)) dξ = (3.8) = ν µ 1-1/α ˆt 0 E α µ 1/α (t -τ ) p (x, τ ) - P (τ ) 2a dτ.
This is an integral equation featuring only one unknown p (x, t). Note that p (x, 0) can be assumed to be known from solving (2.3) with the unknown δ (0) that is found a posteriori from condition (2.2) imposed at t = 0. As we shall further discuss in Subsection 4.2 (see Corollary 8), the assumption on the knowledge of p (x, 0) does not restrict generality as long as η > 0.

The obtained relation (3.8) is an integral equation of the mixed type in two variables: it contains both Fredholm and Volterra integral operators in the spatial and temporal parts, respectively.

Solution by separation of variables

Since temporal and spatial operators appearing in equation (3.8) are not intertwined, it is natural to attempt solving the problem with the method of separation of variables. Therefore, we first focus on the spatial part of the problem and study the appropriate functional setting for its solution.

4.1. Notation, assumptions and their implications. We start by fixing the notation. We shall use L p (-a, a) to denote the space of all real-valued functions on (-a, a) whose p-th power is Lebesgue integrable on the interval (-a, a). In particular, L 2 (-a, a) is a Hilbert space endowed with the inner product and norm defined as

(4.1) f, g := ˆa -a f (x) g (x) dx, f := f, f , f, g ∈ L 2 (-a, a) .
We will also need the space L 2 0 (-a, a) a subspace of L 2 (-a, a) that consists of all square-integrable functions with vanishing mean value on (-a, a), i.e. (4.2)

L 2 0 (-a, a) := f ∈ L 2 (-a, a) : ˆa -a f (x) dx = 0 .
Note that L 2 0 (-a, a) is a closed subspace of L 2 (-a, a), and hence is also a Hilbert space. Given T > 0, let us denote C b ([0, T ]), C b (R + ) the spaces of bounded continuous functions on the closed interval [0, T ] and the positive half-line R + := [0, ∞).

Assumption 1 (Parity and real-valuedness). Suppose that K (x) is a real-valued even function on (-a, a):

K (-x) = K (x) ∈ R, x ∈ (-a, a) .
Assumption 2 (Hilbert-Schmidt regularity). Suppose that K (x) is sufficiently regular, namely, we assume validity of the following integrability condition:

ˆa -a ˆa -a [K (x -ξ)] 2 dxdξ < ∞.
Assumption 3 (Positive semidefiniteness). We suppose that, for any function f ∈ L 2 (-a, a), the following condition holds true:

ˆa -a f (x) ˆa -a K (x -ξ) f (ξ) dξdx ≥ 0.
Because of Assumption 2, the condition of Lemma 17 is satisfied, and therefore the integral transformation given by

(4.3) f → K [f ] (x) := ˆa -a K (x -ξ) f (ξ) dξ defines a compact linear operator K : L 2 (-a, a) → L 2 (-a, a)
. Moreover, due to Assumption 1, this operator is self-adjoint. Consequently, by the spectral theory for compact self-adjoint operators (Lemma 18), it follows that, for any f ∈ L 2 (-a, a), we can write

(4.4) f (x) = ∞ n=1 f n ϕ n (x) + f ⊥ (x) ,
where f ⊥ (x) ∈ Ker K and ϕ n (x), n ≥ 1, are normalised eigenfunctions of the operator K with non-zero eigenvalues.

In other words, we have

(4.5) ˆa -a K (x, ξ) f ⊥ (ξ) dξ = 0, x ∈ (-a, a) , and
(4.6) ˆa -a K (x, ξ) ϕ n (ξ) dξ = λ n ϕ n (x) , x ∈ (-a, a) ,
with some λ n = 0, n ≥ 1, which are eigenvalues (in general, repeated sequentially according to their multiplicity). Setting (4.7)

K 2 (x, ξ) := K (x -ξ) - 1 2a K 1 (x) - 1 2a K 1 (ξ) , K 1 (x) := ˆa -a K (ζ -x) dζ, we observe that (4.8) ˆa -a K 2 (x, ξ) dx = - 1 2a ˆa -a K 1 (x) dx =: c 0 ,
The right-hand side of (4.8) is constant (independent of ξ), and hence, we obtain a relation

(4.9) ˆa -a ˆa -a K 2 (x, ξ) dξdx = - ˆa -a K 1 (x) dx.
Moreover, for any f ∈ L 2 0 (-a, a), we have

ˆa -a ˆa -a K 2 (x, ξ) f (ξ) dξdx = 0.
Therefore, using Lemma 17, we can define K 2 : L 2 0 (-a, a) → L 2 0 (-a, a), a compact linear operator that corresponds to the integral transformation with the kernel function K 2 (x, ξ):

(4.10) f → K 2 [f ] (x) := ˆa -a K 2 (x, ξ) f (ξ) dξ.
By the symmetry of K 2 (x, ξ), the operator K 2 is self-adjoint and hence, using Lemma 18, we deduce that, for any f ∈ L 2 0 (-a, a), we can write

(4.11) f (x) = ∞ n=1 f n φ n (x) + f ⊥ (x) ,
where

f n := f, φ n , f ⊥ := f - ∞ n=1 f n φ n ∈ Ker K 2 , and φ n , n ≥ 1, are normalised (i.e. φ n = 1, n ≥ 1) eigenfunctions of the operator K 2 with non-zero eigenvalues. That is, we have (4.12) ˆa -a K 2 (x, ξ) f ⊥ (ξ) dξ = 0, x ∈ (-a, a) , and
(4.13) ˆa -a K 2 (x, ξ) φ n (ξ) dξ = σ n φ n (x) , x ∈ (-a, a) ,
with some σ n = 0, n ≥ 1, which are eigenvalues, repeated consequently if multiple. Note that ϕ k ⊥ ϕ n (meaning that ϕ k , ϕ n = 0) for any k = n ≥ 1, and similarly, φ k ⊥ φ n , k = n ≥ 1. This orthogonality is automatic when eigenfunctions correspond to different eigenvalues. When belonging to the same eigensubspace, we assume that they have already been orthogonalised (e.g. by the Gram-Schmidt procedure). Moreover, since each φ n ∈ L 2 0 (-a, a), we have

(4.14) 1, φ n = ˆa -a φ n (x) dx = 0, n ≥ 1.
We shall now show deduce more information about eigenvalues of the operators K and K 2 .

Proposition 4. The operator K 2 : L 2 0 (-a, a) → L 2 0 (-a, a) defined by (4.10) is positive semidefinite, and thus σ n ≥ 0.

Proof. First of all, note that according to (4.7), for any f ∈ L 2 0 (-a, a), we have (4.15)

K [f ] -K 2 [f ] , f = 0.
Let us assume that there is at least one eigenvalue σ 1 < 0 (if there are few, we assume that σ 1 is the largest in absolute value) and we shall derive a contradiction.

Observe that, by the variational characterisation of the smallest negative eigenvalue (Rayleigh principle) given by Lemma 19, we have (4. [START_REF] Gakhov | Boundary Value Problems[END_REF])

σ 1 = min f ∈L 2 0 (-a,a) K 2 [f ], f f 2 = min f ∈L 2 0 (-a,a) K[f ], f f 2 ,
where the second equality is due to (4.15).

By positive semidefiniteness of K, we have K[f ], f ≥ 0, and hence equation (4.16) shows that σ 1 ≥ 0 which contradicts the assumption σ 1 < 0. Therefore, all eigenvalues σ n , n ≥ 1, are non-negative, and hence K 2 is a positive semidefinite operator, that is, we have

K 2 [f ], f ≥ 0 for any f ∈ L 2 0 (-a, a).
Proposition 5. Eigenvalues λ n and σ n , n ≥ 1, defined by (4.6) and (4.13), respectively, satisfy the following inequalities:

(4.17)

λ n+1 ≤ σ n ≤ λ n , n ≥ 1.
Proof. By Proposition 4, we have σ n ≥ 0, n ≥ 1. Furthermore, the Weyl-Courant-Fischer min-max principle for characterisation of positive eigenvalues (Lemma 20) gives, for n ≥ 2,

σ n = min q1,...,qn-1∈L 2 (-a,a) max f ∈L 2 0 (-a,a) f ⊥(q1,...,qn-1) K 2 [f ], f f 2 ≤ max f ∈L 2 0 (-a,a) f ⊥(ϕ1,...,ϕn-1) K 2 [f ], f f 2 (4.18) = max f ∈L 2 0 (-a,a) f ⊥(ϕ1,...,ϕn-1) K[f ], f f 2 ≤ max f ∈L 2 (-a,a) f ⊥(ϕ1,...,ϕn-1) K[f ], f f 2 = λ n ,
where we used (4.15) and, in the last equality on the second line, we employed the Rayleigh's variational characterisation for positive eigenvalues (Lemma 19). For n = 1, the situation is even simpler:

σ 1 = max f ∈L 2 0 (-a,a) K 2 [f ], f f 2 = max f ∈L 2 0 (-a,a) K[f ], f f 2 ≤ max f ∈L 2 (-a,a) K[f ], f f 2 = λ 1 . (4.19)
We can also obtain the lower bound estimate for eigenvalues σ n . To this effect, we apply the Weyl-Courant-Fischer min-max principle twice: for λ n and σ n-1 . Namely, we have, for n ≥ 3,

λ n = min q1,...,qn-1∈L 2 (-a,a) max f ∈L 2 (-a,a) f ⊥(q1,...,qn-1) K[f ], f f 2 ≤ min q1,...,qn-2∈L 2 (-a,a) max f ∈L 2 (-a,a) f ⊥(q1,...,qn-2,1) K[f ], f f 2 (4.20) = min q1,...,qn-2∈L 2 (-a,a) max f ∈L 2 0 (-a,a) f ⊥(q1,...,qn-2) K[f ], f f 2 = min q1,...,qn-2∈L 2 (-a,a) max f ∈L 2 0 (-a,a) f ⊥(q1,...,qn-2) K 2 [f ], f f 2 = σ n-1 .
Similarly, we can also obtain:

λ 2 = min q1∈L 2 (-a,a) max f ∈L 2 (-a,a) f ⊥q1 K[f ], f f 2 ≤ max f ∈L 2 0 (-a,a) K[f ], f f 2 = max f ∈L 2 0 (-a,a) K 2 [f ], f f 2 = σ 1 . (4.21)
Combining the estimates obtained in (4.18)-(4.21), we arrive at (4.17). 4.2. Solution existence, uniqueness and construction. Let us set (4.22) q (x, t) := p (x, t) -

P (t) 2a , x ∈ (-a, a) , t ≥ 0,
and observe that

(4.23) ˆa -a q (x, t) dx = 0, t ≥ 0.
Then, equation (3.8) can be equivalently rewritten as

(4.24) ηq (x, t) + K 2 [q] (x, t) - ν µ 1-1/α ˆt 0 E α µ 1/α (t -τ ) q (x, τ ) dτ = F (x, t) , x ∈ (-a, a) , t > 0,
where

K 2 [q] (x, t) := ˆa -a K (x -ξ) - 1 2a K 1 (ξ) - 1 2a K 1 (x) q (ξ, t) dξ, F (x, t) :=ηq (x, 0) + K 2 [q] (x, 0) - P (t) -P (0) 2a (K 1 (x) + c 0 ) , (4.25)
and c 0 is a constant defined in (4.8).

With this preparation, our main results concerning the solution of the model can be formulated as two theorems below. In their statements and proofs, we will, according to the aforementioned convention, employ the simplified notation •, • , • to denote the inner product and the norm in L 2 (-a, a), respectively. Theorem 6. Assume that µ ≥ 0, η, ν > 0, α ∈ (0, 2), P ∈ C b (R + ) and p (•, 0) ∈ L 2 (-a, a). Suppose that K satisfying Assumptions 1-3 is such that the equation

(4.26) K 2 [φ] (x) := ˆa -a K (x -ξ) - 1 2a ˆa -a (K (ζ -x) + K (ζ -ξ)) dζ φ (ξ) dξ = 0, x ∈ (-a, a) ,
has at most one non-zero solution φ ⊥ ∈ L 2 0 (-a, a) with φ ⊥ = 1, i.e. dim (Ker K 2 ) ≤ 1. Then, equations (3.8) and (4.24) have unique solutions in C b R + ; L 2 (-a, a) and C b R + ; L 2 0 (-a, a) , respectively. These solutions are given by

(4.27) p (x, t) = P (t) 2a + ∞ k=1 d k (t) φ k (x) + d ⊥ (t) φ ⊥ (x) , x ∈ (-a, a) , t ≥ 0, (4.28) q (x, t) = ∞ k=1 d k (t) φ k (x) + d ⊥ (t) φ ⊥ (x) , x ∈ (-a, a) , t ≥ 0,
where φ k , σ k , k ≥ 1, are eigenfunctions and eigenvalues of the operator K 2 defined as in (4.13), and

d k (t) :=d 0 k 1 + ν µ (η + σ k ) + ν E α -µ + ν η + σ k t α -1 - l k 2a (η + σ k ) (P (t) -P (0)) (4.29) - νl k 2a (η + σ k ) 2 µ + ν η + σ k 1/α-1 ˆt 0 E α µ + ν η + σ k 1/α (t -τ ) [P (τ ) -P (0)] dτ, k ≥ 1, d ⊥ (t) :=d 0 ⊥ 1 + ν µη + ν E α -µ + ν η t α -1 - l ⊥ 2aη (P (t) -P (0)) (4.30) - νl ⊥ 2aη 2 µ + ν η 1/α-1 ˆt 0 E α µ + ν η 1/α (t -τ ) [P (τ ) -P (0)] dτ, with (4.31) d 0 ⊥ := p (•, 0) , φ ⊥ , l ⊥ := K 1 , φ ⊥ , d 0 k := p (•, 0) , φ k , l k := K 1 , φ k , k ≥ 1.
Proof. Since equations (3.8) and (4.24) are equivalent and simply related, we only consider the latter.

The assumption p (•, 0) ∈ L 2 (-a, a) and (4.22)-(4.23) imply that q (•, 0) ∈ L 2 0 (-a, a). Therefore, following (4.11) and using that dim (Ker K 2 ) ≤ 1, we can write (4.32) q (x, 0) :

= ∞ k=1 d 0 k φ k (x) + d 0 ⊥ φ ⊥ (x) , x ∈ (-a, a) , with d 0 k := q (•, 0) , φ k , k ≥ 1, and d 0 ⊥ := q (•, 0) , φ ⊥ .
Note that these definitions coincide with those in (4.31) due to the mean-zero property (4.14).

For an arbitrary integer m ≥ 1, let us define

(4.33) q m (x, t) := m k=1 d k (t) φ k (x) + d ⊥ (t) φ ⊥ (x) ,
where the functions d k , d ⊥ solve the following integral equations

(4.34) (η + σ k ) d k (t) - ν µ 1-1/α ˆt 0 E α µ 1/α (t -τ ) d k (τ ) dτ = (η + σ k ) d 0 k - l k 2a (P (t) -P (0)) , t > 0, k ≥ 1, (4.35) ηd ⊥ (t) - ν µ 1-1/α ˆt 0 E α µ 1/α (t -τ ) d ⊥ (τ ) dτ = ηd 0 ⊥ - l ⊥ 2a (P (t) -P (0)) , t > 0, with d k (0) = 1 η + σ k ηq (•, 0) + K 2 [q (•, 0)] , φ k = q (•, 0) , φ k = d 0 k , l k := K 1 + c 0 , φ k = K 1 , φ k , k ≥ 1, d ⊥ (0) = d 0 ⊥ , l ⊥ := K 1 + c 0 , φ ⊥ = K 1 , φ ⊥ .
Here, we performed some simplifications of the above expressions due to decomposition (4.32) and used (4.12)-( 4

d k (t) =d 0 k - l k 2a (η + σ k ) (P (t) -P (0)) (4.36) + ν η + σ k µ + ν η + σ k 1/α-1 d 0 k + P (0) 2a l k η + σ k ˆt 0 E α µ + ν η + σ k 1/α (t -τ ) dτ - l k 2a (η + σ k ) ˆt 0 E α µ + ν η + σ k 1/α (t -τ ) P (τ ) dτ , t ≥ 0, k ≥ 1, d ⊥ (t) =d 0 ⊥ - l ⊥ 2aη (P (t) -P (0)) (4.37) + ν η µ + ν η 1/α-1 d 0 ⊥ + P (0) 2a l ⊥ η ˆt 0 E α µ + ν η 1/α (t -τ ) dτ - l ⊥ 2aη ˆt 0 E α µ + ν η 1/α (t -τ ) P (τ )
sup t>0 q (•, t) -q m (•, t) -→ m→∞ 0.
To prove the convergence (4.38), we are going to show that (q m ) ∞ m=1 , with q m given by (4.33), is a Cauchy sequence in C b R + ; L 2 0 (-a, a) . To this effect, let us set (4.39) S mn (x, t) := q n (x, t) -q m (x, t)

= n k=m+1 d k (t) φ k (x) , x ∈ (-a, a) , t ∈ R + , n > m,
and use the orthonormality of φ k to write

S mn (•, t) = n k=m+1 |d k (t)| 2 1/2 = (d k (t)) n k=m+1 l 2 ,
where in the last expression we used • l 2 for designating the Euclidean vector norm. Taking into account that σ k > 0 (recall Proposition 4), we can estimate

ν µ (η + σ k ) + ν E α -µ + ν η + σ k t α -1 ≤ ν µη + ν (C 1 + 1) ,
where

C 1 := sup τ >0 |E α (-τ )| is finite since the function E α (t)
is continuous and decaying for large negative values of the argument t (see (A.5), (A.8), (A.13) and (A.15)). Also, due to the continuity and boundedness of P (t) on R + , it is immediate to see that

1 2a (η + σ k ) |P (t) -P (0)| ≤ 1 aη P L ∞ (R+) , t ∈ R + .
Let us deal with the term on the second line of (4.29). To this end, we have, for any t 0 > 0,

µ + ν η + σ k 1/α-1 ˆt 0 E α µ + ν η + σ k 1/α (t -τ ) [P (τ ) -P (0)] dτ = µ + ν η + σ k 1/α-1 ˆt 0 E α µ + ν η + σ k 1/α τ [P (t -τ ) -P (0)] dτ ≤ ˆt0 0 . . . + ˆt t0 . . . with µ + ν η + σ k 1/α-1 ˆt0 0 E α µ + ν η + σ k 1/α τ [P (t -τ ) -P (0)] dτ ≤ 2C 1,t0 P L ∞ (R+) ˆt0 0 dτ τ 1-α = 2C 1,t0 P L ∞ (R+) t α 0 α , and
µ + ν η + σ k 1/α-1 ˆt0 0 E α µ + ν η + σ k 1/α τ [P (t -τ ) -P (0)] dτ ≤ 2C 2,t0 µ + ν η + σ k -1/α-1 P L ∞ (R+) ˆ∞ t0 dτ τ 1+α = 2C 2,t0 P L ∞ (R+) (µ + ν/ (η + σ 1 )) 1/α+1 αt α 0 .
Here, the constants C 1,t0 := sup τ ∈(0,t0) τ 1-α E α (τ ) and C 2,t0 := sup τ ∈(t0,∞) τ 1+α E α (τ ) are finite due to (A.7), (A.14) and (A.17) since α ∈ (0, 2). Consequently, applying the triangle inequality to (4.36)-(4.37), we estimate

S mn (•, t) ≤ 1 + ν µη + ν (C 1 + 1) d 0 k n k=m+1 l 2
(4.40)

+ 1 + νC 1,t0 t α 0 ηα + νC 2,t0 η (µ + ν/ (η + σ 1 )) 1/α+1 αt α 0 1 aη P L ∞ (R+) (l k ) n k=m+1 l 2 , t ∈ R + .
Now, recalling (4.31), we have, by the Bessel's inequality, , (l k ) n k=m+1 l 2 can be made arbitrary small for large m. Therefore, we deduce from (4.40) that S mn (•, t) is guaranteed to be arbitrary small, uniformly for all t ∈ R + , once sufficiently large value of m is chosen. In other words, recalling (4.39), we have obtained that (q m ) ∞ m=1 is a Cauchy sequence in C b R + ; L 2 0 (-a, a) . Since C b R + ; L 2 0 (-a, a) is a closed subspace of a Banach space (see e.g. [START_REF] Hunter | Applied Analysis[END_REF]Thm 6.28] for the standard fact that L ∞ R + ; L 2 (-a, a) is a Banach space), it is also a Banach space, and hence complete. This implies the desired convergence (4.38).

∞ k=1 d 0 k 2 ≤ q (•, 0) 2 < ∞, ∞ k=1 |l k | 2 ≤ K 1 2 < ∞, i.e. both series ∞ k=1 d 0
Remark 7. The condition dim (Ker K 2 ) ≤ 1 imposed in the formulation of Theorem 6 is essential for the uniqueness of the solution. Otherwise, one could, wihout changing the validity of (3.8), (4.24), add to the solution (4.27)-(4.28) the term d ⊥ (t) φ ⊥ (x) with any φ ⊥ ∈ Ker K 2 , φ ⊥ ⊥ q 0 ⊥ , and d ⊥ solving the integral equation

η d ⊥ (t) - ν µ 1-1/α ˆt 0 E α µ 1/α (t -τ ) d ⊥ (τ ) dτ = - l ⊥ 2a (P (t) -P (0)) , t > 0.
If, on the other hand, K is such that (4.26) admits only the zero solution, i.e. Ker K 2 = 0, then φ ⊥ ≡ 0. Proof. By the standard Fredholm theory applied to the positive compact self-adjoint operator K defined in (4.3), the assumptions ∆ ∈ L 2 (-a, a) and η > 0 entail the existence of a unique p (•, 0) ∈ L 2 (-a, a) satisfying (2.3). This makes Theorem 6 applicable yielding the claim of this Corollary.

Theorem 9. Assume that µ ≥ 0, α ∈ (0, 2), η = 0, ν > 0, P (t) ≡ P (0) =: P 0 for t ∈ [0, T ] with some T > 0, and p (•, 0) ∈ L 2 (-a, a). Suppose that K satisfies Assumptions 1-3 and q (•, 0) := p (•, 0) -P0 2a is orthogonal to any ψ ∈ Ker K 2 . Then, equations (3.8) and (4.24) have unique solutions in C [0, T ] ; L 2 (-a, a) and C [0, T ] ; L 2 0 (-a, a) , respectively. These solutions are given by

(4.41) p (x, t) = P 0 2a + ∞ k=1 d k (t) φ k (x) , q (x, t) = ∞ k=1 d k (t) φ k (x) , x ∈ (-a, a) , t ∈ (0, T ) ,
where φ k , σ k , k ≥ 1, are eigenfunctions and eigenvalues of the operator K 2 defined as in (4.13), and

d k (t) :=d 0 k 1 + ν µσ k + ν E α -µ + ν σ k t α -1 , k ≥ 1, (4.42) with d 0 k , k ≥ 1,
as in (4.31). Proof. The proof generally goes along the same lines as that for Theorem 6. Similarly to (4.33), we introduce (4.43)

q m (x, t) := m k=1 d k (t) φ k (x) + d ⊥ (t) φ ⊥ (x) , m ≥ 1,
where φ ⊥ ∈ L 2 0 (-a, a) denotes any solution of the equation (4.26). By substitution of (4.43) in (4.24) and using orthogonality of φ k , k ≥ 1, and φ ⊥ , we deduce that d k , d ⊥ are the functions solving the following integral equations (4.44)

σ k d k (t) - ν µ 1-1/α ˆt 0 E α µ 1/α (t -τ ) d k (τ ) dτ = σ k d 0 k , t ∈ (0, T ) , k ≥ 1, (4.45) - ν µ 1-1/α ˆt 0 E α µ 1/α (t -τ ) d ⊥ (τ ) dτ = 0, t ∈ (0, T ) .
We note immediately that application of Corollary 27 or Lemma 28 to (4.45), depending on whether α ∈ (0, 1) or α ∈ [1, 2), yields d ⊥ (t) ≡ 0, t ∈ (0, T ). For the solution of (4.44), we invoke Corollary 25 which gives

d k (t) =d 0 k 1 + ν/σ k (µ + ν/σ k ) 1-1/α ˆt 0 E α µ + ν σ k 1/α τ dτ , t ∈ (0, T ) . (4.46)
Making use of (A.11), equation (4.46) simplifies into (4.42). The proof of the convergence of (4.43) as m → ∞ is identical to that for Theorem 6 with an additional simplification that there is now no necessity to deal with the l k terms.

Finally, we note that the orthogonality condition q (•, 0) ⊥ ψ for any ψ ∈ Ker K 2 imposed in the formulation of the Theorem is due to the requirement of the continuity of the solution (4.28). Indeed, if this condition was violated, we would have

q (x, 0) = ∞ k=1 d k (t) φ k (x) + ψ 0 (x) , x ∈ (-a, a) ,
for some ψ 0 ∈ Ker K 2 . This would not be consistent with the limit of (4.28) as t → 0 + (which itself is a consequence of our conclusion that d ⊥ (t) ≡ 0, t ∈ (0, T ), in (4.43)).

Remark 10. It is noteworthy that, in case η = 0, the regularity requirement p (•, 0) ∈ L 2 (-a, a) may be a rather strong one, when viewed in the context of the entire problem (2.1)-(2.3). As we shall see on an example of a particular kernel function K in Proposition 12, this may restrict ∆ in (2.3) and P 0 . On the other hand, the same example shows that the orthogonality condition in the statement of Theorem 9 may be trivially satisfied.

4.3.

A concrete form of the kernel function. We now make results of the previous subsection more precise by focussing on a concrete kernel function that is often used for contact mechanical problems under consideration, namely, problems with a half-space geometry or for an elastic layer of a large thickness (see e.g. [START_REF] Aleksandrov | Mathematical methods in problems with wear[END_REF][START_REF] Argatov | Effective wear coefficient and wearing-in period for a functionally graded wear-resisting punch[END_REF][START_REF] Shtaerman | Contact problem of the theory of elasticity[END_REF]). Namely, we consider the kernel function given by (4.47)

K 0 (x) := -log |x| + C K , x ∈ (-2a, 2a) ,
where C K > log a is a constant. It is straightforward to see that K = K 0 satisfies Assumptions 1-2. Verification of Assumption 3 requires a change of variable x = ax, the additive property of logarithms and the positive-definiteness result of [START_REF] Reade | Asymptotic behavior of eigen-values of certain integral equations[END_REF] valid for the integral operator with purely logarithmic kernel (i.e. C K = 0) on the interval (-1, 1). Here, we have also made use of the assumption C K > log a.

We now claim that the eigenvalues of the corresponding K 2 operator can be characterised as follows.

Proposition 11. Let K = K 0 with K 0 defined in (4.47). Then, the eigenvalues of the operator K 2 defined in (4.10) satisfy the following estimates

(4.48) 0 < σ 1 ≤ aπ log 2 + 2a (C K + |log a|) , 0 < σ n = O 1 n , n ≫ 1.
Proof. First, we claim that λ n , the eigenvalues of the operator K, defined in (4.3), decrease to zero as O (1/n) for large n. This follows from the corresponding result of [START_REF] Reade | Asymptotic behavior of eigen-values of certain integral equations[END_REF] obtained for the case of purely logarithmic kernel (i.e. (4.47) with C K = 0). Indeed, since the asymptotic decrease of the eigenvalues of a positive integral operator is related to the regularity of the kernel function (see, in addition to [START_REF] Reade | Asymptotic behavior of eigen-values of certain integral equations[END_REF], also [START_REF] Reade | Eigen-values of Lipschitz kernels[END_REF]), the presence of an extra constant C K > 0 does not affect this asymptotic behavior. The final asymptotic result given in (4.48) now follows from λ n = O (1/n), n ≫ 1, by employing Proposition 5.

To deduce the upper bound for σ 1 , we shall first get the one for λ 1 . To this effect, we transform (4.6) to an equivalent problem on the interval (-1, 1). Namely, setting ψ (x) := φ (ax), we have

(4.49) ˆ1 -1 [-log |x -ξ| + log a + C K ] ψ (ξ) dξ = λ a ψ (x) , x ∈ (-1, 1) .
Rayleigh's variational characterisation for positive eigenvalues (Lemma 19) now yields

λ 1 a ≤ max f ∈L 2 (-1,1) ´1 -1 -log |x -ξ| f (ξ) f (x) dξdx f L 2 (-1,1) + (C K + |log a|) max f ∈L 2 (-1,1) ´1 -1 |f (ξ)| |f (x)| dξdx f L 2 (-1,1) ≤ π log 2 + 2 (C K + |log a|) ,
where we used π log 2 as an upper bound for the first eigenvalue of the logarithmic kernel due to [START_REF] Reade | Asymptotic behavior of eigen-values of certain integral equations[END_REF], and we used the Cauchy-Schwarz inequality to trivially estimate the second quotient. Finally, using Proposition 5, we obtain the bound for σ 1 in (4.48).

When η = 0, the particular form of the kernel function K 0 makes it possible to work with an explicit form of the initial data p (x, 0). Indeed, we have the following constructive result. Proposition 12. Let η = 0, a = 2 and K = K 0 with K 0 defined in (4.47). Suppose that ∆ ∈ C 2 ([-a, a]). Then, a solution of (2.3) subject to (2.2) with t = 0 is given by (4.50)

p (x, 0) = 1 (a 2 -x 2 ) 1/2 a -a a 2 -ξ 2 1/2 ∆ ′ (ξ) ξ -x dξ + 1 log (a/2) ˆa -a ∆ (ξ) (a 2 -ξ 2 ) 1/2 dξ + π (C K P (0) -δ (0)) ,
and (4.51)

δ (0) = 1 π ˆa -a ∆ (ξ) (a 2 -ξ 2 ) 1/2 dξ-P (0) 1 π 2 log a 2 -C K + 1 π 2 log a 2 ˆa -a 1 (a 2 -x 2 ) 1/2 a -a a 2 -ξ 2 1/2 ∆ ′ (ξ) ξ -x dξdx.
Proof. Inserting (4.47) and taking into account (2.2) with t = 0, (2.3) rewrites as

- ˆa -a log |x -ξ| p (ξ, 0) dξ = -C K P (0) + δ (0) -∆ (x) , x ∈ (-a, a) .
Application of Lemma 29 using the elementary identity

(4.52) ˆa -a dx (a 2 -x 2 ) 1/2 = ˆ1 -1 dx (1 -x 2 ) 1/2 = π, now yields p (x, 0) = 1 (a 2 -x 2 ) 1/2 a -a a 2 -ξ 2 1/2 ∆ ′ (ξ) ξ -x dξ (4.53) + 1 (a 2 -x 2 ) 1/2 1 log (a/2) ˆa -a ∆ (ξ) (a 2 -ξ 2 ) 1/2 dξ + π (C K P (0) -δ (0)) ,
where ffl a -a stands for the Cauchy principal value integral. Integrating (4.53) on (-a, a), we reuse, in the left-hand side, (2.2) with t = 0. Rearranging the terms, we obtain (4.51).

Remark 13. We see, from a form of the solution (4.50), that, in general, p (•, 0) / ∈ L 2 (-a, a). However, the squareintegrability condition can be achieved for some particular profiles ∆ and values P (0), namely, those that annihilate the square bracket in (4.50) at x = ±a, see e.g. [17, pp. 47-48] for examples of finite pressure distributions for quadratic and quartic symmetric shapes of ∆.

Finally, by focussing on the concrete kernel function K 0 , we will show that the auxiliary condition Ker K 2 ≤ 1 in Theorem 6 and the orthogonality condition in Theorem 9 are not difficult to verify. Proposition 14. Assume that a = 2. Let K = K 0 with K 0 defined in (4.47). Then, the only solution of (4.26) in L 2 0 (-a, a) is φ ≡ 0. Proof. Since we look for the solution in L 2 0 (-a, a), we use zero-mean condition (4.23) to rewrite (4.26) as

ˆa -a K 0 (x -ξ) φ (ξ) dξ = 1 2a ˆa -a K 0 (ζ -ξ) φ (ξ) dζdξ, x ∈ (-a, a) ,
and, furthermore,

(4.54) - ˆa -a log |x -ξ| φ (ξ) dξ = - 1 2a ˆa -a ˆa -a log |ζ -ξ| φ (ξ) dζdξ, x ∈ (-a, a) .
Since the right-hand side of (4.54) is just a constant, application of Lemma 29 yields a particularly simple result

φ (x) = 1 (a 2 -x 2 ) 1/2 π 2a log (a/2) ˆa -a ˆa -a log |ζ -ξ| φ (ξ) dξdζ, x ∈ (-a, a) ,
where we used (4.52). Upon further integration over (-a, a) and use of (4.23), we conclude that we must have ˆa -a ˆa -a log |ζ -ξ| φ (ξ) dζdξ = 0.

Getting back to (4.54), we see that this condition entails that

ˆa -a log |x -ξ| φ (ξ) dξ = 0, x ∈ (-a, a) .
Hence, by applying Lemma 29 again, we conclude that φ ≡ 0.

Analysis of the solution

We are going to show that the solution to the proposed model has features reflecting the expected physical behavior. In particular, consider two settings: a constant load and a transitional load (i.e. the one which stabilises to a constant value after a finite time). We show that, in both cases, for large times, the solution stabilises to a stationary pressure distribution that can be found explicitly. Moreover, in some cases, this pressure distribution is simply constant (uniform), an aspect which is consistent with previous models [START_REF] Argatov | Wear contact problem with friction: Steady-state regime and wearing-in period[END_REF][START_REF] Argatov | A Macro-Scale Approximation for the Running-In Period[END_REF] but is certainly not a general feature (see e.g. [START_REF] Goryacheva | Contact Mechanics in Tribology[END_REF]Ch. 6]).

5.1. Constant load. First, let us consider the most commonly investigated case of a constant load, i.e. where P (t) ≡ P (0) =: P 0 for t ≥ 0. Proposition 15. Assume that µ, η ≥ 0, α ∈ (0, 2), ν > 0, P (t) ≡ P (0) =: P 0 , t ≥ 0, and p (•, 0) ∈ L 2 (-a, a). Suppose that K satisfying Assumptions 1-3 is such that equation (4.26) has at most one solution φ ⊥ ∈ L 2 0 (-a, a) with φ ⊥ = 1. Moreover, if η = 0, we additionally assume that p (•, 0) -P0 2a , φ ⊥ = 0. Then, for the solution of (3.8), we have the following asymptotic results (4.31). Proof. First, let us consider η > 0. Application of Theorem 6, yields

(5.1) p (•, t) -p (1) ∞ = O exp -µ + ν η + σ 1 t , t ≫ 1, α = 1, (5.2) p (•, t) -p (1) ∞ = O 1 t α , t ≫ 1, α ∈ (0, 1) ∪ (1, 2) , with (5.3) p (1) ∞ (x) := P 0 2a + ∞ k=1 µ (η + σ k ) µ (η + σ k ) + ν d 0 k φ k (x) + µη µη + ν d 0 ⊥ φ ⊥ (x) , x ∈ (-a, a) ,
and d 0 k , k ≥ 1, d 0 ⊥ , as in
p (x, t) - P 0 2a = ∞ k=1 d k (t) φ k (x) + d ⊥ (t) φ ⊥ (x) , x ∈ (-a, a) , t ≥ 0,
or equivalently, rearranging the terms so that the right-hand side contains only those proportional to E α ,

p (x, t) -p (1) ∞ (x) = ∞ k=1 d k (t) φ k (x) + d ⊥ (t) φ ⊥ (x) , x ∈ (-a, a) , t ≥ 0, with p (1)
∞ defined as in (5.3), and

(5.4)

d k (t) := d k (t) - µ (η + σ k ) µ (η + σ k ) + ν d 0 k = ν µ (η + σ k ) + ν d 0 k E α -µ + ν η + σ k t α , t ≥ 0, k ≥ 1, (5.5) d ⊥ (t) := d ⊥ (t) - µη µη + ν d 0 ⊥ = ν µη + ν d 0 ⊥ E α -µ + ν η t α , t ≥ 0.
Note that the series in (5.3) converges in L 2 (-a, a) due to the Parseval's identity, since

∞ k=1 µ (η + σ k ) µ (η + σ k ) + ν d 0 k φ k 2 = ∞ k=1 η + σ k µ (η + σ k ) + ν d 0 k 2 ≤ sup n≥1 µ (η + σ n ) µ (η + σ n ) + ν 2 ∞ k=1 d 0 k 2 < ∞ k=1 d 0 k 2 ≤ q (•, 0) 2 < ∞,
and hence we have p

(1) ∞ ∈ L 2 (-a, a).
By the orthonormality of φ k , k ≥ 1, and φ ⊥ , we have

(5.6) p (•, t) -p (1) ∞ = ∞ k=1 d k (t) 2 + d ⊥ (t) 2 1/2 .
We can estimate

d k (t) ≤ d 0 k sup n≥1 ν µ (η + σ n ) + ν sup n≥1 E α -µ + ν η + σ n t α ≤ d 0 k E α -µ + ν η + σ 1 t α , t ≫ 1, k ≥ 1,
where, in the second line, we used 0 < σ n ≤ σ 1 , n ≥ 1, together with the fact that |E α (-τ )| is monotonically decreasing for sufficiently large τ (as evident from asymptotic expansion (A.15)). Consequently, (

) p (•, t) -p (1) ∞ 2 ≤ ∞ k=1 d 0 k 2 E α -µ + ν η + σ 1 t α 2 + d 0 ⊥ 2 E α -µ + ν η t α 2 , t ≫ 1. 5.6) implies (5.7
When α ∈ (0, 1)∪(1, 2), the use of (A.15) in (5.7) immediately gives (5.2). When α = 1, we have E 1 (-τ ) = exp (-t) (see (A.8)), and we observe that the first term in the square bracket of (5.7) is dominant as it decays slower for t ≫ 1. This yields (5.1). Now, when η = 0, we use Theorem 9 which gives

p (x, t) -p (1) ∞ (x) = ∞ k=1 d k (t) φ k (x) , x ∈ (-a, a) , t ≥ 0, with p (1)
∞ and d k defined as before, according to (5.3) and (5.4), respectively. This leads to an analog of (5.7), namely,

p (•, t) -p (1) ∞ 2 ≤ ν 2 (µη + ν) 2 ∞ k=1 d 0 k 2 E α -µ + ν η + σ 1 t α 2 , t ≫ 1,
and hence estimate (5.1) or (5.2) follows depending on the value of α.

5.2. Transitional load. We now consider the second scenario, when the load is of transitional type, i.e. P is taken to be a continuous function on R + with P (0) =: P 0 and such that P (t) ≡ P 1 for t ≥ t 0 with some t 0 > 0.

Proposition 16. Assume that µ ≥ 0, α ∈ (0, 2), η, ν > 0, P ∈ C b (R + ) with P (0) =: P 0 and such that P (t) ≡ P 1 , t ≥ t 0 , for some t 0 > 0, and p (•, 0) ∈ L 2 (-a, a). Suppose that K satisfying Assumptions 1-3 is such that equation (4.26) has at most one solution φ ⊥ ∈ L 2 0 (-a, a) with φ ⊥ = 1. Then, for the solution of (3.8), we have the following asymptotic results

(5.8) p (•, t) -p (2) ∞ = O exp -µ + ν η + σ 1 t , t ≫ 1, α = 1, (5.9
) p (•, t) -p (2) ∞ = O 1 t α , t ≫ 1, α ∈ (0, 1) ∪ (1, 2) , with (5.10) p (2) ∞ (x) := P 1 2a + ∞ k=1 c (2)
k φ k (x) + c (2)
⊥ φ ⊥ (x) , x ∈ (-a, a) , where (5.11) c (2)
k := µ (η + σ k ) µ (η + σ k ) + ν d 0 k + ν (η + σ k ) (µ (η + σ k ) + ν) - 1 η + σ k (P 1 -P 0 ) l k 2a , k ≥ 1, (5.12) c (2)
⊥ := µη µη + ν d 0 ⊥ + ν µη + ν -1 (P 1 -P 0 ) l ⊥ 2aη ,
and d 0 k , k ≥ 1, d 0 ⊥ are as in (4.31). Proof.
The proof is very similar to the one of Proposition 15 with the main difference that the separation of terms in (4.29)-(4.30) into the constant and the time-decaying components is now slightly more complicated. Namely, for the present choice of the load function P , from Theorem 6, we have

d k (t) =d 0 k 1 + ν µ (η + σ k ) + ν E α -µ + ν η + σ k t α -1 - l k 2a (η + σ k ) (P (t) -P 0 ) (5.13) - νl k 2a (η + σ k ) 2 µ + ν η + σ k 1/α-1 ˆt0 0 E α µ + ν η + σ k 1/α (t -τ ) P (τ ) dτ + P 0 2a νl k (η + σ k ) (µ (η + σ k ) + ν) E α -µ + ν η + σ k t α -1 - P 1 2a νl k (η + σ k ) (µ (η + σ k ) + ν) E α -µ + ν η + σ k (t -t 0 ) α -1 , k ≥ 1, d ⊥ (t) =d 0 ⊥ 1 + ν µη + ν E α -µ + ν η t α -1 - l ⊥ 2aη (P (t) -P 0 ) (5.14) - νl ⊥ 2aη 2 µ + ν η 1/α-1 ˆt0 0 E α µ + ν η 1/α (t -τ ) P (τ ) dτ + P 0 2a νl ⊥ η (µη + ν) E α -µ + ν η t α -1 - P 1 2a νl ⊥ η (µη + ν) E α -µ + ν η (t -t 0 ) α -1 .
Adding and subtracting l k 2a(η+σ k ) P 1 and l ⊥ 2aη P 1 in (5.13) and (5.14), respectively, we rearrange the terms to arrive at p (x, t) -p (2) ∞ (x) =

P (t) -P 1 2a + ∞ k=1 dk (t) φ k (x) + d⊥ (t) φ ⊥ (x) , x ∈ (-a, a) , t ≥ 0, with p (2)
∞ defined as in (5.10), and

dk (t) := d 0 k + P 0 2a l k (η + σ k ) ν µ (η + σ k ) + ν E α -µ + ν η + σ k t α (5.15) - ν (η + σ k ) (µ (η + σ k ) + ν) E α -µ + ν η + σ k (t -t 0 ) α P 1 l k 2a - (P (t) -P 1 ) l k 2a (η + σ k ) - νl k 2a (η + σ k ) 2 µ + ν η + σ k 1/α-1 ˆt0 0 E α µ + ν η + σ k 1/α (t -τ ) P (τ ) dτ, t ≥ 0, k ≥ 1, d⊥ (t) := d 0 ⊥ + P 0 2aη l ⊥ ν µη + ν E α -µ + ν η t α (5.16) - ν µη + ν E α -µ + ν η (t -t 0 ) α P 1 l ⊥ 2aη - (P (t) -P 1 ) l ⊥ 2aη - νl ⊥ 2aη 2 µ + ν η 1/α-1 ˆt0 0 E α µ + ν η 1/α (t -τ ) P (τ ) dτ, t ≥ 0.
We note that here again the series in (5.10) converges in L 2 (-a, a) because of

∞ k=1 c (2) k 2 < ∞ which, in turn, follows from ∞ k=1 d 0 k 2 < ∞, ∞ k=1 |l k | 2 < ∞, since µ(η+σ k )
µ(η+σ k )+ν < 1 and the square-bracketed term in (5.11) is uniformly bounded for all k ≥ 1 and t ≥ 0.

To deduce estimates (5.8)-(5.9), we consider

p (•, t) -p (2) ∞ - P (t) -P 1 2a 2 = ∞ k=1 dk (t) 2 + d⊥ (t) 2 (5.17)
≤ 2 ∞ k=1 |D k (t)| 2 d 0 k 2 + |L k (t)| 2 |l k | 2 + 2 |D 0 (t)| 2 d 0 ⊥ 2 + 2 |L 0 (t)| 2 |l ⊥ | 2 ,
where we used the elementary inequality (a + b) 2 ≤ 2 a 2 + b 2 , a, b ∈ R, and we introduced

(5.18) D 0 (t) := ν µη + ν E α -µ + ν η t α , D k (t) := ν µ (η + σ k ) + ν E α -µ + ν η + σ k t α , k ≥ 1, L 0 (t) := P 0 2aη ν µη + ν E α -µ + ν η t α - P 1 2aη ν µη + ν E α -µ + ν η (t -t 0 ) α (5.19) - P (t) -P 1 2aη - ν 2aη 2 µ + ν η 1/α-1 ˆt0 0 E α µ + ν η 1/α (t -τ ) P (τ ) dτ, L k (t) := P 0 2a ν (η + σ k ) (µ (η + σ k ) + ν) E α -µ + ν η + σ k t α
(5.20)

- P 1 2a ν (η + σ k ) (µ (η + σ k ) + ν) E α -µ + ν η + σ k (t -t 0 ) α - P (t) -P 1 2a (η + σ k ) - ν 2a (η + σ k ) 2 µ + ν η + σ k 1/α-1 ˆt0 0 E α µ + ν η + σ k 1/α (t -τ ) P (τ ) dτ, k ≥ 1.
From asymptotics (A.15), (A.17), it follows that |E α (-τ )| and |E α (τ )| are monotonically decreasing functions for sufficiently large τ . Since 0 < σ k ≤ σ 1 , k ≥ 1, and P (t) ≡ P 1 for t ≥ t 0 , we can estimate from (5.18)-(5.20)

|D 0 (t)| ≤ E α -µ + ν η t α , |D k (t)| ≤ E α -µ + ν η + σ 1 t α , t ≫ 1, k ≥ 1,
|L 0 (t)| ≤ P 0 2aη E α -µ + ν η t α + P 1 2aη E α -µ + ν η (t -t 0 ) α + ν 2aη 2 µ + ν η 1/α-1 E α µ + ν η 1/α (t -t 0 ) P L 1 (0,t0) , t ≫ 1, |L k (t)| ≤ P 0 2a (η + σ 1 ) E α -µ + ν η + σ 1 t α + P 1 2a (η + σ 1 ) E α -µ + ν η + σ 1 (t -t 0 ) α + ν 2aη (η + σ 1 ) sup n≥1 µ + ν η + σ n 1/α-1 E α µ + ν η + σ 1 1/α (t -t 0 ) P L 1 (0,t0) , t ≫ 1, k ≥ 1,
For α ∈ (0, 1) ∪ (1, 2), we reuse (A.15), (A.17) to get (5.9). For α = 1, employing (A.6)-(A.8) which entails that E 1 (-τ ) = |E 1 (τ )| = exp (-τ ), τ > 0, we deduce (5.8). In both cases, we again used that P (t) ≡ P 1 for t ≥ t 0 which entailed identical vanishing of the third term in the left-hand side of (5.17).

Numerical illustrations

To verify and illustrate some of the obtained results numerically, we consider the particular kernel function K 0 as given by (4.47), with C K = 1.6. We use a specially written MATLAB code which employs an external function [START_REF] Podlubny | Mittag-Leffler function[END_REF] for computing E α . 6.1. Verification of the solution for a variable load. Let us fix the following set of parameters a = 1, ν = 2, η = 1, µ = 1.2, α = 0.6. Consider the load profile given by (6.1)

P (t) = P 1 χ (t0,∞) (t) + P 0 + P 1 2 - P 1 -P 0 2 cos πt t 0 χ (0,t0) (t)
, t ≥ 0, P 0 = 6, P 1 = 10, t 0 = 0.5 consistent with Subsection 5.2. This describes a smooth load switch from P 0 to P 1 occuring over time t 0 and remaining constant afterwards. For the sake of simplicity, let us suppose that the function ∆ in (2.3) is chosen such that (6.2)

p 0 (x) = 2P 0 a 2 π a 2 -x 2 , x ∈ (-a, a) ,
and note that, for such a choice, (2.2) is automatically satisfied. We verify the approach presented in Section 4 by comparing the solution formula given in Theorem 6, with d ⊥ (t) ≡ 0, d 0 ⊥ = l ⊥ = 0 (according to Proposition 14), to a numerical finite-difference method for solving integral equation (4.24) in time. In particular, we truncate the series in (4.28) at 60 terms (however, much less would already be sufficient). For the numerical approach, we combine Nyström collocation method with a finite-difference scheme in time using singularity subtraction. In Figure 6.1, we see that both solutions almost coincide for all shown instances of time. The solutions have an oscillatory mismatch, especially in small regions close to the endpoints x = ±1, which is a typical phenomenon for a spectral approach. Let us consider the same setup as in Subsection 6.1, but instead of a fixed value of µ, we explore the range of this model parameter starting from the degenerate case µ = 0 (purely fractional order model the wear term (2.4) given by (3.6)) up to µ = 6. This is done in order to investigate the effect of such a parameter variation on an essential output of the model: the stationary pressure distribution p

∞ (x) given by (5.10). As Figure 6.2 shows, the increase of µ amounts to steepening of the curve p (2) ∞ (x) making a deviation from the uniform pressure distribution more pronounced. Note that, according to (5.10), the stationary pressure distribution is independent of the model parameter α. 6.3. Illustration of the convergence under a constant load. We now consider a simple setting which is more classical for analysis: the constant load case, i.e. we replace (6.1) with P (t) ≡ P 0 , t ≥ 0. We take (6.2) and numerical values of P 0 and other parameters as in Subsection 6.1, except for the model "order" parameter α which we now vary. The goal here is to demonstrate qualitatively different behaviour of the model depending on this parameter. In particular, we investigate three different values of α, by looking at the pressure evolution at two particular points. For α = 1, the convergence to the stationary value occurs fast, and hence Figure 6.3 corroborates exponential bound (5.1). For all other values of α inside the interval (0, 2), the stabilisation occurs only at an algebraic rate. In particular, this is illustrated in Figures 6.4 and 6.5 where values α = 0.6 and α = 1.2 are taken, respectively. We see that in the first case convergence happens monotonically whereas in the second one it is accompanied by oscillations. Figure 6.6 shows, on extended time interval, that for a larger value (α = 1.8), oscillations increase even more. Several formulations of the punch-sliding problem have been considered at once. These included rather general relation between contact pressure and displacement generalising the half-space geometry, a potential presence of coating (or additional surface microstructure), time-varying load and a linear material wear. In particular, the wear was modelled through a novel non-local relation between the contact pressure and the wear generalising the classical one. One advantage of transient models with such a wear relation is that they admit a closed-form solution as long as the spatial part of the formulation can be easily resolved and, most importantly, analysed. This possibility should not be underestimated in contexts of inverse design of mechanical materials.

Another relevant advantage is that the solution may exhibit predictably different transient behavior and the stationary pressure distribution depending on newly introduced parameters entering the formulation. As it was shown numerically, these parameters can be chosen to ensure a slow algebraic stabilisation towards the stationary distribution whether its monotone or oscillatory. This complements the fast exponential stabilisation (which happens in the particular case of the model parameter α = 1), a typically occuring phenomenon for such problems, and paves a way for description of new materials such as polymers within the same framework.

The non-uniform stationary pressure distributions are also possible and correspond to non-zero values of the model parameter µ. At the next stage, a practical validation of the obtained results and a fit of the model to experimental data is highly desirable. In doing this, the parameter µ should be chosen from an observation of the stationary pressure distribution whereas the parameter α should be determined from temporal observations (the convergence speed), according to what was described above.

We thus conclude that the proposed model has a potential to describe materials whose temporal convergence to a stationary state is slow as well as those for which the stationary pressure distribution is not uniform.

As a continuation of this work, analysis of the long-time behavior under a periodic load P has already been considered in [START_REF] Ponomarev | Generalised model of wear in contact problems: the case of oscillatory load[END_REF]. More challenging would be to deal rigorously with more general kernel functions K such as those that do not satisfy the parity assumption. However, it seems that the positivity assumption is not essential, and, in view of that, some technical results in Appendix are already provided for a slightly more general setting than the one being considered here.

Lemma 19. [35, Sect. 95 Thm p. 237] Let X be a Hilbert space with the inner product •, • X , and A : X → X is a compact self-adjoint operator. Then, λ + n , the n-th largest positive eigenvalue of A can be characterised as

λ + 1 = max f ∈X A [f ] , f X f, f X , λ + n = max f ∈X f ⊥(ϕ + 1 ,...ϕ + n-1 ) A [f ] , f X f, f X , n ≥ 2,
where ϕ + 1 , . . . ϕ + n-1 are the eigenfunctions corresponding to the n-1 largest positive eigenvalues, and the orthogonality condition f ⊥ ϕ + 1 , . . . ϕ + n-1 means that f, ϕ + j X = 0, j = 1, . . ., n -1. Similarly, λ - n , the n-th smallest negative eigenvalue of A can be characterised as

λ - 1 = min f ∈X A [f ] , f X f, f X , λ - n = min f ∈X f ⊥(ϕ - 1 ,...ϕ - n-1 ) A [f ] , f X f, f X , n ≥ 2,
where ϕ - 1 , .

A [f ] , f X f, f X , n ≥ 2,
where the orthogonality condition f ⊥ (f 1 , . . . f n-1 ) means that f, f j X = 0, j = 1, . . ., n -1.

A.2. Some special functions and their properties.

Gamma function:

(A.2) Γ (z) := ˆ∞ 0 x z-1 e -x dx, z > 0.

This function satisfies the following fundamental relation Note that when λ > 0, we can write (A.9) e α (z; λ) = λ 1/α E α λ 1/α z , z > 0, α > 0, λ > 0.

We have a useful relation (A.10) E α,α (z) = αE ′ α,1 (z) = αE ′ α (z) , z ∈ C, α > 0, which is straightforward to obtain using definitions (A.5) and (A.3). Also, from (A.9) and (A.6), it follows (A.11) ˆz0 0 E α λ 1/α z dz = λ -1/α ˆz0 0 e α (z; λ) dz = λ -1/α [E α (-λz α 0 ) -1] , z 0 > 0, α > 0, λ > 0.

If f ′ is absolutely continuous on [0, T ], we can apply the same procedure to the integral on the right-hand side before another differentiation and then we arrive at

d 2 dt 2 ˆt 0 f (τ ) (t -τ ) α-1 dτ = - α -1 t α f (0) + f ′ (0) t α-1 + ˆt 0 f ′′ (τ ) (t -τ ) α-1 dτ,
which, upon insertion into (A.26), gives (A.27). Now, we consider the situation when α = 1. Equation (A.23) degenerates and its unique solution u (t) = f ′ (t) follows immediately by differentiation of both sides of the equation. Taking into account that sin[π(α-1)] π(α-1) → 1 as α → 1, upon an elementary simplification, we see that this solution coincides precisely with (A.27) for α = 1. (µ + λ) 1-1/α ˆt 0 E α (µ + λ) 1/α (t -τ ) f (τ ) dτ, t > 0.

Proof. For µ > 0, the result follows from Lemma 24 applied to the integral equation for u t/µ 1/α with λ/µ and f t/µ 1/α in place of λ and f (t), respectively. To treat the case µ = 0, we observe, using (A.14) and (A.3), that in the limit µ → 0, equation (A.34) becomes (A.21) whose solution, given by Lemma 21, coincides exactly with (A.35) for µ = 0. Lemma 26. Let α ∈ (0, 1), and assume that f is an absolutely continuous function on every subinterval of R + . Then, the integral equation (A. [START_REF] Samko | Fractional Integrals and Derivatives -Theory and Applications[END_REF] ˆt 0 E α (t -τ ) u (τ ) dτ = f (t) , t > 0, has a unique solution given by

u (t) = - f (0) Γ (1 -α) t α - ˆt 0 1 + 1 Γ (1 -α) 1 (t -τ ) α f ′ (τ ) dτ (A.37) = -f (t) - 1 Γ (1 -α) t α -1 f (0) - 1 Γ (1 -α) ˆt 0 f ′ (τ ) (t -τ ) α dτ, t > 0.
Proof. The proof is ideologically similar to that of Lemma 24. Let us denote U (s) := L [u] (s), F (s) := L [f ] (s) the Laplace transforms of u (t) and f (t). Upon Laplace transformation of (A.36), using the convolution theorem and (A.18), we have

-U (s) 1 s α + 1 = F (s) =⇒ U (s) = -F (s) (1 + s α ) = -[sF (s) -f (0)] 1 s + 1 s 1-α -f (0) 1 s + 1 s 1-α .
Recalling (A.14) and (A.3), we observe that in the limit µ → 0, equation (A.38) becomes (A.23). Solution of the latter equation, for α ∈ [1, 2), is given by Corollary 23. Using the identities Γ (α -1) Γ (2 -α) = π/ sin [π (α -1)] and Γ (α) = (α -1) Γ (α -1), this solution can be seen to coincide with (A.40) for µ = 0. Proof. See e.g. [12, p. 428] for the case a = 1 and [16, p. 591] for arbitrary a > 0, a = 2.

Figure 2 .

 2 Figure 2.1. Geometry of the problem. A balance of displacements under the punch is governed by the equation: ηp(x, t) + ´a -a K (x -ξ) p (ξ, t) dξ = δ (t) -∆ (x) -w (x, t), x ∈ (-a, a), t ≥ 0.

Corollary 8 .

 8 Assume that µ ≥ 0, α, η, ν > 0, P ∈ C b (R + ) and ∆ ∈ L 2 (-a, a). Then, the problem given by (2.1)-(2.3), with K satisfying Assumptions 1-3 and such that dim (Ker K 2 ) ≤ 1, and w [p] defined as in (2.4), is uniquely solvable and the solution p ∈ C b R + , L 2 (-a, a) is furnished by (4.27).

Figure 6 . 1 .

 61 Figure 6.1. Time evolution of the pressure profile computed from equation (4.27) (left) and by a numerical finite-difference method (right)

Figure 6 . 2 .Figure 6 . 3 . α = 1 : 6 . 2 .

 6263162 Figure 6.2. Effect of the model parameter µ on the stationary state p(2) ∞ (x)

Figure 6 . 4 . α = 0. 6 :Figure 6 . 5 .Figure 6 . 6 . 7 .

 64665667 Figure 6.4. α = 0.6: Time evolution of the pressure at x = 0 (left) and x = 0.5 (right)

(A. 3 )(- 1 )

 31 Γ (z + 1) = zΓ (z) , Γ (n + 1) = n!, z > 0, n ∈ N 0 , as well as Euler's reflection formula(A.4) Γ (1 -z) Γ (z) = π sin (πz) k kλ k z αk Γ (αk + 1) , z > 0, α > 0, λ ∈ C, (A.7) E α (z) := e α (z; 1) = α z ∞ k=1 k kz αk Γ (αk + 1) , z > 0, α > 0.The following two identities are direct consequences of definitions (A.5) and (A.3)(A.8) E α,1 (z) = E α (z) , E 1 (z) = exp (z) , z ∈ C.

Lemma 24 .Corollary 25 .

 2425 Let α > 0, λ ∈ R. Then, the integral equation(A.31) u (t) -λ ˆt 0 E α (t -τ ) u (τ ) dτ = f (t) , t > 0,has a unique solution given by(A.32) u (t) = f (t) + λ 1 + λ ˆt 0 e α (t -τ ; 1 + λ) f (τ ) dτ, t > 0. Proof. Let us denote U (s) := L [u] (s), F (s) := L [f ] (s)the Laplace transforms of u (t) and f (t), respectively. Application of the Laplace transformation to the both sides of (A.31) yields, upon using the convolution theorem and (A.18), (A.33)U (s) + λ s α + 1 U (s) = F (s) =⇒ U (s) = s α + 1 s α + 1 + λ F (s) = 1 -λ s α + 1 + λ F (s) .Using (A.[START_REF] Gorenflo | Mittag-Leffler Functions, Related Topics and Applications[END_REF]), inversion of Laplace transformation of the rightmost equation above leads to (A.32). Let α > 0, µ ≥ 0, λ ∈ R. Then, the integral equation(A.34) u (t) -λ µ 1-1/α ˆt 0 E α µ 1/α (t -τ ) u (τ ) dτ = f (t) , t > 0,has a unique solution given by u (t) = f (t) + λ µ + λ ˆt 0 e α (t -τ ; µ + λ) f (τ ) dτ (A.35) = f (t) + λ

Lemma 29 . 1 (a 2 -x 2 )a 2 -

 291222 Suppose that f ∈ C 2 ([-a, a]), then, if a = 2, the integral equation -ˆa -a log |x -ξ| u (ξ) dξ = f (x) , x ∈ (-a, a) ,admits a solution given byu (x) = -ξ 2 1/2 f ′ (ξ) ξ -x dξ + 1 log (a/2) ˆa -a f (ξ) (a 2 -ξ 2 ) 1/2 dξ , x ∈ (-a, a) .Moreover, this solution is unique in the class of Hölder continuous functions with possible integrable singularities at the endpoints of the interval [-a, a].

  .13) as well as (4.14). We observe that equations (4.34)-(4.35) can be solved in a closed form by application of Corollary 25. Using (A.11), this yields

  . . ϕ - n-1 are the eigenfunctions corresponding to the n -1 smallest negative eigenvalues.Lemma 20. [35, Sect. 95 Thm p. 237] Let X be a Hilbert space with the inner product •, • X , and A : X → X is a compact self-adjoint operator. Then, λ n , the n-th largest eigenvalue of A can be characterised as

	λ n =	min f1,...,fn-1∈X	max f ∈X
			f ⊥(f1,...fn-1)

FACTAS team, Centre Inria d'Université Côte d'Azur, France

St. Petersburg Department of Steklov Mathematical Institute of Russian Academy of Sciences, Russia

Contact: dmitry.ponomarev@inria.fr

Acknowledgements

The author is grateful to the support of the bi-national FWF-project I3538-N32 of the Austrian Science Fund (FWF) used for his employment at Vienna University of Technology. The paper has benefited from fruitful discussions with Ivan Argatov on mechanical aspects of the matter.

Appendix A.

We collect here some auxiliary results that are needed in the paper.

A.1. Basic theoretical facts about compact linear integral operators. Lemma 17. Let K ∈ L 2 ((-a, a) × (-a, a)), i.e. such that

Then, the integral transformation f → K [f ] (x) := ´a -a K (x, ξ) f (ξ) dξ is a compact linear operator from L 2 (-a, a) to L 2 (-a, a).

Proof. Under assumption of the validity of (A.1), it follows that K [f ] ∈ L 2 (-a, a) by an application of the Cauchy-Schwarz inequality (see also [START_REF] Hackbusch | Integral Equations: Theory and Numerical Treatment[END_REF]Lem 3.2.3]). The fact that the operator K is compact can be shown [START_REF] Hackbusch | Integral Equations: Theory and Numerical Treatment[END_REF]Thm 3.2.7] using Weierstrass approximation theorem and the characterisation of a compact operator as the limit of finite rank operators.

Lemma 18. Let X be a Hilbert space and A : X → X is a compact self-adjoint operator. Then, the eigenvalues of A form a real non-increasing in absolute value sequence, and every eigenvalue different from zero has finite multiplicity. Moreover, the eigenfunctions of A form an orthonormal basis in Ran A, the closure of the range of A.

Proof. The statement of the lemma is a collection of standard results of the spectral theory of compact self-adjoint operators. A close result in the form of a single theorem is given in [START_REF] Atkinson | Theoretical Numerical Analysis -A Functional Analysis Framework[END_REF]Thm 2.8.15]. See also [START_REF] Riesz | Functional Analysis[END_REF]Sect. 97 Small-argument asymptotics of the above functions follow directly from definitions (A.5)-(A.7)

Large-argument asymptotics of E α (z; λ) can be derived from that of E α (z) given, for example, in [19, eqs. (3.4.14)-(3.4. [START_REF] Feppon | Efficient steady-state computation for wear of multimaterial composites[END_REF]] (see also [29, eq. (1.

2)] for the same results for α ∈ (0, 1))

We also need some Laplace transforms of the above functions. We denote the Laplace transform of f as L [f ] (s) := ´∞ 0 f (t) e -st dt. In particular, we have (see [20, eq. (3.14)] or [27, eq. ( 4)])

which implies that, for λ > 0,

To obtain similar result for λ < 0, we use the formula [36, eq. (1.93)]

which implies that

Here, we used the identity (A.10) and the definition of the function E a . Note that (A.20) shows that (A. [START_REF] Gorenflo | Mittag-Leffler Functions, Related Topics and Applications[END_REF]) is actually valid for all λ ∈ R.

A.3. Some singular integral equations and their solutions. We consider here Abel's integral equations of the first and the second kinds and two other related equations.

Lemma 21. [19, Thm 4.2] Let α > 0, λ ∈ C and f ∈ L 1 (0, T ) for any T > 0. Then, the integral equation

has a unique solution u ∈ L 1 (0, T ) given by

Note that the integral form of the solution on the second line of (A. [START_REF] Hackbusch | Integral Equations: Theory and Numerical Treatment[END_REF]) can be seen due to the identity (A.10), and it is also given elsewhere (see [19, eq. (8.1.20)], [20, eq. (2.12)]).

Lemma 22. [36, Thm 2.1 & Lem. 2.1] (see also [START_REF] Gorenflo | Mittag-Leffler Functions, Related Topics and Applications[END_REF]Sect. 8.1.1]) Let α ∈ (0, 1), and assume that f is such that ´t 0 f (τ ) (t-τ ) α dτ is an absolutely continuous function on [0, T ] for some T > 0 and, moreover, lim t→0 + ´t 0 f (τ ) (t-τ ) α dτ = 0. Then, the integral equation

has a unique solution u ∈ L 1 (0, T ) given by

In particular, to satisfy the aforementioned conditions on f , it is sufficient that f is an absolutely continuous function on [0, T ]. In this case, the solution (A.24) can be written in the alternative form

, and assume that f is such that ´t 0 f (τ ) (t-τ ) α-1 dτ is an absolutely continuous function on [0, T ] for some T > 0 and, moreover, lim t→0 + ´t 0 f (τ ) (t-τ ) α-1 dτ = 0. Then, integral equation (A.23) has a unique solution u such that ´• 0 u (τ ) dτ ∈ L 1 (0, T ) and it is given by

Moreover, if f ′ is absolutely continuous on [0, T ], solution formula (A.26) can be rewritten as

Proof. Let us first assume that α ∈ (1, 2). Denoting u (t) := ´t 0 u (τ ) dτ , integration by parts yields

where the boundary terms vanish due to α > 1 and u (0) = 0. Lemma 22 now applies to the equation

and gives the existence of a unique solution ũ ∈ L 1 (0, T ):

Hence, differentiating and using the identity

)), we obtain (A.26). Suppose now that f is absolutely continuous on [0, T ], we can then integrate by parts, taking into account that α < 2,

and thus

Hence, employing L t β-1 (s) = Γ (β) /s β , β > 0 (see e.g. [13, pp. 318-319]), and, in particular, L [START_REF] Alblas | Contact problems of a rectangular block on an elastic layer of finite thickness: the thin layer[END_REF] (s) = 1/s, we invert the transform to obtain (A.37).

Corollary 27. Let α ∈ (0, 1), µ ≥ 0 and assume that f is an absolutely continuous function on every bounded subinterval of R + . Then, the integral equation

has a unique solution given by

Proof. For µ > 0, the result follows from Lemma 26 applied to the integral equation for u t/µ 1/α with -µf t/µ 1/α in place of f (t).

To treat the case µ = 0, we observe, using (A.14) and (A.3), that in the limit µ → 0, equation (A.38) becomes (A.23) whose solution, given by Lemma 22, coincides exactly with (A.39) for µ = 0. Lemma 28. Let α ∈ [1, 2), µ ≥ 0, and assume that f ′ is an absolutely continuous function on every bounded subinterval of R + . Then, integral equation (A.38) has a unique solution given by

First, let us consider α ∈ (1, 2). Let u (t) := ´t 0 u (τ ) dτ . Since E α (0) = 0 for α > 1 and u (0) = 0, integration by parts of (A.38) leads to

the Laplace transforms of u (t) and f (t), respectively, we apply Laplace transformation to the above equation and use the convolution theorem and L E ′ α µ 1/α • (s) = -µ 1-2/α s s α +µ (which easily follows from (A.18)) to arrive at

Therefore, upon inversion of Laplace transformation, using that L t 1-α (s) = Γ (2 -α) /s 2-α , α < 2, (and, in particular, L [t] (s) = 1/s 2 ), we obtain

Note that, similarly to (A.29)-(A.30), we have, for α ∈ (1, 2),

and hence, by differentiation of (A.41), we deduce (A.40). It remains to consider the situation when α = 1. To this effect, we recall (A.9) and (A.6) which imply that E 1 (µt) = -e -µt , and hence equation (A.38) can be recast as ˆt 0 e µτ u (τ ) dτ = e µt f (t) , t > 0, which is then immediately solved by the differentiation to give u (t) = µf (t) + f ′ (t). This solution coincides with that given by (A.40), upon substitution α = 1 and relevant simplifications (integration and cancellations). Case µ = 0: