**1. Introduction**

During the last decades, fractional differential equations and their potential applications have gained a lot of importance, mainly because fractional calculus has become a powerful tool with more accurate and successful results when modeling several complex phenomena in numerous seemingly diverse and widespread fields of science and engineering [1]. It was found that various, especially interdisciplinary applications, can be elegantly modeled with the help of fractional derivatives which provide an excellent instrument for the description of memory and hereditary properties of various materials and processes [2,3]. Advanced analysis and numerical simulations of several fractional-order systems have been shown to be very interesting, producing more useful results in applied sciences [4,5].

Delay differential equations are a type of equations in which the derivative of the unknown function at a certain time is given in terms of the values of the function at previous times. It arises in many biological and physical applications, and it often forces us to consider variable or state-dependent delays [6–8]. Integer or fractional-order degenerate differential equations, i.e., evolution equations not solved with respect to the highest order derivative, are often used to describe various processes in science and engineering: in [9,10] certain classes of the time-fractional order partial differential equations with polynomials differential with respect to the spatial variables elliptic self-adjoint operator, which contain some equations from hydrodynamics and the filtration theory, are studied. In [11] approximate controllability issues for such models are investigated; the unique solvability of similar equations with distributed order time derivatives are researched in [12].

In applications, fractional-order degenerate evolution equations with a delay are often successful. Such kinds of equations with a degenerate operator at the highest-order fractional derivative describe the dynamics of some fractional models of viscoelastic fluids (see the application in the last section of

this work). There are very few papers dealing with essentially degenerate fractional-order equations with delay. Motivated by this fact, the purpose of this work is a step towards eliminating this gap.

We are concerned with the following fractional differential equations with delay

$$LD\_t^\kappa \mathbf{x}(t) = M\mathbf{x}(t) + \int\_{-r}^0 \mathcal{K}(\mathbf{s})\mathbf{x}(t+s)ds + \mathbf{g}(t), \quad t \in [0, T], \tag{1}$$

where X , Y are Banach spaces, *L*, *M* : X→Y are linear operators, *L* is continuous, ker *L* = {0} (for this reason such equations are called Sobolev type equations [13,14], or degenerate [15]), operator *M* is closed and densely defined in <sup>X</sup> , *<sup>D</sup><sup>α</sup> <sup>t</sup>* is the Gerasimov-Caputo derivative of the order *α* ∈ (*m* − 1, *m*], *<sup>m</sup>* <sup>∈</sup> <sup>N</sup>. Equation (1) is endowed by a given background

$$P\mathbf{x}(t) = h(t), \quad t \in [-r, 0], \tag{2}$$

and by the generalized Showalter–Sidorov conditions

$$(Px)^{(k)}(0) = \mathbf{x}\_k, \; k = 0, 1, \ldots, m - 1,\tag{3}$$

which are natural for degenerate evolution equations. Here, *P* is a projector along the degeneration space of the homogeneous equation *LD<sup>α</sup> <sup>t</sup> x*(*t*) = *Mx*(*t*), it will be defined below. By the contraction mappings method, the local unique solvability of problems (1)–(3) is established.

Degenerate first-order evolution equations in Banach spaces were studied in [16,17] under various conditions on the operators *L*, *M* and on the delay term. The unique solvability results for problems (1) and (2) with a strongly (*L*, *p*)-radial operator *M*, *g* ≡ 0 at *α* = 1 were obtained in [18]. Here we use a similar approach, which is adapted to the case of a fractional derivative. The second section contains the preliminary results which are needed for supporting our results, in particular, the theorem on unique solvability of the Cauchy problem to the inhomogeneous linear Equation (1) with K ≡ 0. In the third section, we obtain the proof of the main result by means of the Banach fixed point theorem. The fourth and fifth sections demonstrate the applications of the obtained abstract results to the study of the unique solvability of initial-boundary value problems for time-fractional systems of partial differential equations with delay.

#### **2. Solvability of Degenerate Inhomogeneous Equation**

$$\text{Let for } \delta > 0 \text{, } t > 0 \text{ g}\_{\delta}(t) := \Gamma(\delta)^{-1} t^{\delta - 1}, \text{ } l\_t^{\delta} h(t) := \underset{0}{\int}{\text{g}\_{\delta}} (t - s) h(s) ds, \text{ } m - 1 < a \le m \in \mathbb{N}\_{\delta}$$

*D<sup>m</sup> <sup>t</sup>* is the usual derivative of the order *<sup>m</sup>* <sup>∈</sup> <sup>N</sup>, *<sup>J</sup>*<sup>0</sup> *<sup>t</sup>* be the identical operator. The Gerasimov-Caputo derivative of a function *h* is defined as

$$D\_t^\alpha h(t) = D\_t^m J\_t^{m-\alpha} \left( h(t) - \sum\_{k=0}^{m-1} h^{(k)}(0) g\_{k+1}(t) \right).$$

**Lemma 1.** *Ref. [19]. Let* <sup>Z</sup> *be a Banach space, l* <sup>−</sup> <sup>1</sup> <sup>&</sup>lt; *<sup>β</sup>* <sup>≤</sup> *<sup>l</sup>* <sup>∈</sup> <sup>N</sup>*, t* <sup>&</sup>gt; <sup>0</sup>*. Then*

$$\exists \mathcal{C}\_{\mathcal{S}} > 0 \quad \forall h \in \mathcal{C}^{l}([0, t]; \mathcal{Z}) \quad \|D\_{t}^{\mathcal{S}}h\|\_{\mathcal{C}([0, t]; \mathcal{Z})} \leq \mathcal{C}\_{\mathcal{S}} \|h\|\_{\mathcal{C}^{l}([0, t]; \mathcal{Z})}.$$

For Banach spaces, X and Y denote as L(X ; Y) the Banach space of all linear continuous operators, acting from X to Y. Let C*l*(X ; Y) be the set of all linear closed operators, densely defined in X , with the image in Y.

In further consideration, we will assume that *L* ∈ L(X ; Y), ker *L* = {0}, *M* ∈ C*l*(X ; Y), *DM* is the domain of *<sup>M</sup>* with the graph norm ·*DM* :<sup>=</sup> ·X <sup>+</sup> *<sup>M</sup>* · Y. Denote *<sup>ρ</sup>L*(*M*) :<sup>=</sup> {*<sup>μ</sup>* <sup>∈</sup> <sup>C</sup> : (*μ<sup>L</sup>* <sup>−</sup> *<sup>M</sup>*)−<sup>1</sup> ∈ L(Y; <sup>X</sup> )}. An operator *<sup>M</sup>* is called (*L*, *<sup>σ</sup>*)-bounded, if a ball *Ba*(0) :<sup>=</sup> {*<sup>μ</sup>* <sup>∈</sup> <sup>C</sup> : <sup>|</sup>*μ*<sup>|</sup> <sup>&</sup>lt; *<sup>a</sup>*} with some *a* > 0 contains the set *ρL*(*M*). If *M* is (*L*, *σ*)-bounded, we have the projections

$$P := \frac{1}{2\pi i} \int\_{|\mu|=a} \int \left(\mu L - M\right)^{-1} L \, d\mu \in \mathcal{L}(\mathcal{X}), \ Q := \frac{1}{2\pi i} \int\_{|\mu|=a} L \left(\mu L - M\right)^{-1} d\mu \in \mathcal{L}(\mathcal{Y}).$$

(see [14] (pp. 89–90)). Set <sup>X</sup> <sup>0</sup> :<sup>=</sup> ker *<sup>P</sup>*, <sup>X</sup> <sup>1</sup> :<sup>=</sup> im*P*, <sup>Y</sup><sup>0</sup> :<sup>=</sup> ker *<sup>Q</sup>*, <sup>Y</sup><sup>1</sup> :<sup>=</sup> im*Q*. Denote by *Lk* (or *Mk*) the restriction of the operator *<sup>L</sup>* (or *<sup>M</sup>*) on <sup>X</sup> *<sup>k</sup>* (or *DMk* :<sup>=</sup> *DM* ∩ X *<sup>k</sup>* respectively), *<sup>k</sup>* <sup>=</sup> 0, 1.

**Theorem 1.** *Ref. [14] (pp. 90–91). Let an operator M be* (*L*, *σ*)*-bounded. Then*

*(i) M*<sup>1</sup> ∈ L <sup>X</sup> 1; <sup>Y</sup><sup>1</sup> *, M*<sup>0</sup> ∈ C*l* <sup>X</sup> 0; <sup>Y</sup><sup>0</sup> *, Lk* ∈ L <sup>X</sup> *<sup>k</sup>*; <sup>Y</sup>*<sup>k</sup> , k* = 0, 1; *(ii) there exist operators M*−<sup>1</sup> <sup>0</sup> ∈ L <sup>Y</sup>0; <sup>X</sup> <sup>0</sup> *, L*−<sup>1</sup> <sup>1</sup> ∈ L <sup>Y</sup>1; <sup>X</sup> <sup>1</sup> *.*

Denote <sup>N</sup><sup>0</sup> :<sup>=</sup> {0} ∪ <sup>N</sup>, *<sup>G</sup>* :<sup>=</sup> *<sup>M</sup>*−<sup>1</sup> <sup>0</sup> *<sup>L</sup>*0. For *<sup>p</sup>* <sup>∈</sup> <sup>N</sup><sup>0</sup> operator *<sup>M</sup>* is called (*L*, *<sup>p</sup>*)-bounded, if it is (*L*, *<sup>σ</sup>*)-bounded, *<sup>G</sup><sup>p</sup>* <sup>=</sup> 0, *<sup>G</sup>p*+<sup>1</sup> <sup>=</sup> 0.

Consider the degenerate inhomogeneous equation

$$LD\_t^a \mathbf{x}(t) = M\mathbf{x}(t) + f(t), \quad t \in [0, T]. \tag{4}$$

A solution of this equation is a function *<sup>x</sup>* <sup>∈</sup> *<sup>C</sup>*([0, *<sup>T</sup>*]; *DM*), such that *<sup>D</sup><sup>α</sup> <sup>t</sup> x* ∈ *C*([0, *T*]; X ) and equality (4) holds. A solution of the generalized Showalter–Sidorov problem

$$(Px)^{(k)}(0) = \mathbf{x}\_k \quad k = 0, 1, \ldots, m - 1,\tag{5}$$

to Equation (4) is a solution of the equation, such that conditions (5) are true.

Denote by *<sup>E</sup>α*,*β*(*z*) = <sup>∞</sup> ∑ *n*=0 *zn* <sup>Γ</sup>(*αn*+*β*) the Mittag-Leffler function.

**Theorem 2.** *Refs. [20,21]. Let <sup>p</sup>* <sup>∈</sup> <sup>N</sup>0*, an operator <sup>M</sup> be* (*L*, *<sup>p</sup>*)*-bounded, Q f* <sup>∈</sup> *<sup>C</sup>*([0, *<sup>T</sup>*]; <sup>Y</sup>)*, for all l* = 0, 1, ... , *p there exist* (*GD<sup>α</sup> t* )*l M*−<sup>1</sup> <sup>0</sup> (*<sup>I</sup>* <sup>−</sup> *<sup>Q</sup>*)*<sup>f</sup>* , *<sup>D</sup><sup>α</sup> <sup>t</sup>* (*GD<sup>α</sup> t* )*l M*−<sup>1</sup> <sup>0</sup> (*I* − *Q*)*f* ∈ *C*([0, *T*]; X ), *<sup>x</sup>*0, *<sup>x</sup>*1,..., *xm*−<sup>1</sup> ∈ X <sup>1</sup>*. Then, problems* (4) *and* (5) *have a unique solution*

$$\mathbf{x}\_f(t) = \sum\_{k=0}^{m-1} t^k E\_{\mathbf{a},k+1}(L\_1^{-1} M\_1 \mathbf{t}^\mathbf{a}) \mathbf{P} \mathbf{x}\_k + \int\_0^t E\_{\mathbf{a},\mathbf{a}}(L\_1^{-1} M\_1 (t-s)^\mathbf{a}) L\_1^{-1} Q f(s) ds - \sum\_{l=0}^p (G D\_l^\mathbf{a})^l M\_0^{-1} (I - Q) f(t). \tag{6}$$

**Remark 1.** *Due to Lemma <sup>1</sup> a function f* <sup>∈</sup> *<sup>C</sup>m*(*p*+1)([0, *<sup>T</sup>*]; <sup>Y</sup>) *satisfies the conditions of Theorem 2.*
