Solution Operators for Fractional Order Coupled Systems Governed by Distributed Order Equations

1. Introduction

Fractional-order systems have attracted considerable attention due to their ability to incorporate memory and hereditary effects that are often observed in real-world phenomena. Unlike classical integer-order models, fractional differential equations account for the influence of the entire past history of a process, providing a more realistic description of dynamics in complex media. Such memory-coupled models arise naturally in anomalous diffusion [16,27,32], viscoelasticity [19,20], heat conduction with memory [11,26], biological and ecological systems [8,12,18], finance [5,25], and control theory [22,23]. The nonlocal nature of fractional operators allows one to capture long-range temporal correlations and multiscale relaxation mechanisms, making fractional-order systems a powerful framework for modeling processes whose evolution depends not only on their current state but also on their accumulated history.

In the study of linear systems of evolution equations, solution operators play a central role by mapping initial data and forcing terms to the corresponding solution at later times. For a classical linear evolution equation, such as $u_t=Au$ in a Banach space, the solution can be expressed formally as $u\left(t\right)=S\left(t\right)u_0$, where $S\left(t\right)$ is a strongly continuous semigroup generated by the operator A [6,14]. These operators encapsulate the temporal dynamics of the system, allowing one to analyze stability, long-time behavior, and regularity properties of solutions through the properties of the semigroup itself. In particular, semigroup theory provides powerful tools for the characterization of solutions in both abstract and concrete settings, ranging from the heat equation to wave propagation problems.

In the fractional setting, however, a direct extension of the classical semigroup framework is generally not available due to the nonlocal nature of fractional derivatives. Consequently, the traditional treatment of fractional evolution equations relies on two distinct solution operators: one governing the homogeneous problem and another describing the contribution of the inhomogeneous forcing term; e.g. [1,2,3,13,15,17] and references therein. The fractional Duhamel principle overcomes this difficulty by restoring a unified representation of solutions, allowing both the initial data and the forcing term to be incorporated through a single solution operator; see [29,30,31].

More precisely, the fractional Duhamel principle represents solutions of inhomogeneous problems in terms of homogeneous solution operators and fractional impulse terms. Within this framework, the influence of past states is transmitted through convolution of the forcing term with suitable fractional kernels, thereby capturing the memory effects inherent in the system. As a result, the associated solution operators can be viewed as natural extensions of classical semigroups, with the standard composition law replaced by convolution structures reflecting the accumulated influence of the system history. The construction of such operators is therefore fundamental to the representation and analysis of both linear and nonlinear fractional evolution systems. Furthermore, this framework reveals intricate couplings among different fractional derivatives in multi-term and coupled systems, providing a natural mechanism for describing memory effects, interactions, and the temporal evolution of fractional dynamics.

In this paper, we construct and analyze solution operators for a broad class of fractional-order systems. The considered models include uniformly and heterogeneously distributed-order differential equation (DODE) systems, memory-decoupled systems, as well as fully coupled systems. Such models arise naturally in the study of diffusion processes exhibiting multiple dynamical regimes, as well as in systems whose components interact through both differential operators and memory effects. Section 2 presents the preliminary material used throughout the paper. In Section 3, Section 4 and Section 5, we derive representations of solution operators for uniformly and heterogeneously distributed-order systems, memory-decoupled systems, and fully coupled systems, and investigate their properties. The illustrative examples are also provided.

2. Preliminaries

In this section we recall briefly some basic notions and results from fractional calculus that will be used throughout the paper. The details can be found e.g. in [24,28]. Let $f:[0,\infty )\to \mathbb{R}$ be a locally integrable function and let $\alpha >0$. The Riemann-Liouville fractional integral of order $\alpha$ is defined by

$$I^{\alpha }f\left(t\right)={\displaystyle \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}}{\int }_0^t{(t-\tau )}^{\alpha -1}f\left(\tau \right)d\tau ,t>0,$$

where $\mathsf{\Gamma }(·)$ denotes the Euler gamma function. For convenience, we set $I^{0}f=f$. The family ${\left\{I^{\alpha }\right\}}_{\alpha \ge 0}$ satisfies the semigroup property

$$I^{\alpha }{I}^{\beta }f=I^{\alpha +\beta }f,\alpha ,\beta \ge 0.$$

Let $\alpha \in (0,1)$. The Riemann–Liouville fractional derivative of order $\alpha$ is defined by

$$D_{+}^{\alpha }f\left(t\right)={\displaystyle \frac{d}{dt}}\left(I^{1-\alpha }f\right)\left(t\right)={\displaystyle \frac{1}{\mathsf{\Gamma }(1-\alpha )}}{\displaystyle \frac{d}{dt}}{\int }_0^t{(t-\tau )}^{-\alpha }f\left(\tau \right)d\tau .$$

The Caputo–Djrbashian fractional derivative of order $\alpha$ is given by

$$D_{*}^{\alpha }f\left(t\right)=I^{1-\alpha }\left({\displaystyle \frac{df}{dt}}\right)\left(t\right)={\displaystyle \frac{1}{\mathsf{\Gamma }(1-\alpha )}}{\int }_0^t{(t-\tau )}^{-\alpha }{f}^{\prime }\left(\tau \right)d\tau .$$

More generally, if $n-1<\alpha <n$ with $n\in \mathbb{N}$, then

$$D_{+}^{\alpha }f\left(t\right)={\displaystyle \frac{d^n}{d{t}^n}}{I}^{n-\alpha }f\left(t\right),$$

and

$$D_{*}^{\alpha }f\left(t\right)=I^{n-\alpha }{\displaystyle \frac{d^{n}f}{d{t}^n}}\left(t\right).$$

The Caputo–Djrbashian derivative is particularly convenient in applications because it allows the use of classical initial conditions involving integer-order derivatives.

Let $\mu$ be a finite positive Borel measure supported on $[0,1]$. The corresponding distributed-order Caputo differential operator is defined by

$${\mathbb{D}}_{\mu }f\left(t\right)={\int }_0^1{D}_{*}^{\alpha }f\left(t\right)\mu \left(d\alpha \right),$$

provided the integral exists. If the measure $\mu$ is absolutely continuous with respect to the Lebesgue measure, namely $\mu \left(d\alpha \right)=w\left(\alpha \right)d\alpha ,$ where $w\in L^1\left([0,1]\right)$ and $w\left(\alpha \right)\ge 0$, then

$${\mathbb{D}}_{\mu }f\left(t\right)={\int }_0^{1}w\left(\alpha \right)D_{*}^{\alpha }f\left(t\right)d\alpha .$$

Distributed-order operators naturally arise in models exhibiting a continuum of memory scales and are widely used in anomalous diffusion and relaxation phenomena. Associated with $\mu$ we introduce the function

$$\psi \left(s\right)={\int }_0^1{s}^{\alpha }\mu \left(d\alpha \right),s>0,$$

(1)

which plays the role of the symbol of the distributed-order operator.

For a suitable function f, its Laplace transform is defined by

$$\widehat{f}\left(s\right)=\mathcal{L}\left\{f\right\}\left(s\right)={\int }_0^{\infty }{e}^{-st}f\left(t\right)dt,s>0.$$

The Laplace transform of the Riemann-Liouville fractional integral is

$$\mathcal{L}\left\{I^{\alpha }f\right\}\left(s\right)=s^{-\alpha }\widehat{f}\left(s\right),\alpha >0.$$

For $\alpha \in (0,1)$, the Laplace transform of the Caputo-Djrbashian derivative is

$$\mathcal{L}\left\{D_{*}^{\alpha }f\right\}\left(s\right)=s^{\alpha }\widehat{f}\left(s\right)-s^{\alpha -1}f\left(0\right).$$

(2)

Similarly, the Riemann-Liouville derivative satisfies

$$\mathcal{L}\left\{D_{+}^{\alpha }f\right\}\left(s\right)=s^{\alpha }\widehat{f}\left(s\right)-\left(I^{1-\alpha }f\right)\left(0^{+}\right).$$

For the distributed-order operator, Fubini’s theorem yields

$$\mathcal{L}\left\{{\mathbb{D}}_{\mu }f\right\}\left(s\right)=\psi \left(s\right)\widehat{f}\left(s\right)-{\displaystyle \frac{\psi \left(s\right)}{s}}f\left(0\right),$$

where $\psi \left(s\right)$ defined in (1).

The Mittag-Leffler function plays a fundamental role in the theory of fractional differential equations. For $\alpha >0$, it is defined by

$$E_{\alpha }\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \frac{z^k}{\mathsf{\Gamma }(\alpha k+1)}},z\in \mathbb{C}.$$

More generally, for $\alpha >0$ and $\beta \in \mathbb{R}$, the two-parameter Mittag–Leffler function is given by

$$E_{\alpha ,\beta }\left(z\right)={\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \frac{z^k}{\mathsf{\Gamma }(\alpha k+\beta )}},z\in \mathbb{C}.$$

Both series define entire functions on the complex plane. In the special case $\alpha =1$, one has $E_1\left(z\right)=e^z.$ The Mittag-Leffler functions satisfy the Laplace transform formulas

$$\mathcal{L}\left\{t^{\beta -1}{E}_{\alpha ,\beta }\left(\lambda t^{\alpha }\right)\right\}\left(s\right)={\displaystyle \frac{s^{\alpha -\beta }}{s^{\alpha }-\lambda }},\Re \left(s\right)>{\left|\lambda \right|}^{1/\alpha },$$

and, in particular,

$$\mathcal{L}\left\{E_{\alpha }\left(\lambda t^{\alpha }\right)\right\}\left(s\right)={\displaystyle \frac{s^{\alpha -1}}{s^{\alpha }-\lambda }}.$$

Furthermore, for $\lambda >0$ and $0<\alpha <1$, the function $t\mapsto E_{\alpha }(-\lambda t^{\alpha })$ is completely monotone on $(0,\infty )$ and satisfies the asymptotic relation

$$E_{\alpha }(-\lambda t^{\alpha })\sim {\displaystyle \frac{1}{\lambda \mathsf{\Gamma }(1-\alpha )}}{t}^{-\alpha },t\to \infty .$$

We will use the following lemmas in the proofs of main results.

Lemma 1

([9,21,28]). Let $0<\alpha <1$ and $\lambda <0$. Then

$$E_{\alpha }\left(\lambda t^{\alpha }\right)={\int }_0^{\infty }{e}^{-rt}{K}_{\alpha }(r;\lambda )dr,t>0,$$

where

$$K_{\alpha }(r;\lambda )=-{\displaystyle \frac{\lambda sin\left(\pi \alpha \right)}{\pi }}{\displaystyle \frac{r^{\alpha -1}}{r^{2\alpha }-2\lambda r^{\alpha }cos\left(\pi \alpha \right)+{\lambda }^2}},r>0.$$

Moreover,

$$K_{\alpha }(r;\lambda )\ge 0,{\int }_0^{\infty }{K}_{\alpha }(r;\lambda )dr=1.$$

Consequently, the function $t\mapsto E_{\alpha }\left(\lambda t^{\alpha }\right)$ is completely monotone on $(0,\infty )$.

Lemma 2

([10,14]). Let $T>0$ and let K be the Volterra operator defined by

$$\left(Ku\right)\left(t\right)={\int }_0^{t}K(t,s)u\left(s\right)ds,0\le t\le T,$$

acting on $C\left(\right[0,T\left]\right)$. Assume that $\left|K\right(t,s\left)\right|\le M$, $0\le s\le t\le T$, for some constant $M>0$. Let $f\in C\left(\right[0,T\left]\right)$. Then the integral equation $u=f+Ku$ has a unique solution $u\in C\left(\right[0,T\left]\right)$, given by the Neumann series

$$u={\displaystyle \sum _{n=0}^{\infty }}{K}^{n}f,$$

where the series converges uniformly on $[0,T]$. Moreover, for every $n\in \mathbb{N}$, we have

$$|\left(K^{n}f\right)\left(t\right){|\le \parallel f\parallel }_{\infty }{\left(Mt\right)}^n/n!,0\le t\le T,$$

and consequently ${\parallel u\parallel }_{\infty }\le e^{MT}{\parallel f\parallel }_{\infty }$.

Lemma 3

([4,7]). Let $f:(0,\infty )\to (0,\infty )$ be a locally integrable function such that $f\left(t\right)=t^{\rho -1}L\left(t\right),$ where $\rho \ge 0$ and $L\left(t\right)$ is a slowly varying function at infinity in the sense of Karamata. Then, as $s\to 0^{+}$, the Laplace transform of $f\left(t\right)$ satisfies

$$\widehat{f}\left(s\right)\sim \mathsf{\Gamma }\left(\rho \right)s^{-\rho }L\left({\displaystyle \frac{1}{s}}\right).$$

Conversely, if $\widehat{f}\left(s\right)$ has this asymptotic form at $s\to 0^{+}$, then $f\left(t\right)$ satisfies

$$f\left(t\right)\sim {\displaystyle \frac{t^{\rho -1}}{\mathsf{\Gamma }\left(\rho \right)}}L\left(t\right),t\to \infty .$$

3. Solution Operators Associated with DODE Systems

Fractional derivatives of multiple orders arise naturally in the modeling of various dynamical processes. In this section, we study solution operators for systems with distributed-order time-fractional derivatives in an abstract setting. In general, such systems take the form

$$D_{\mathit{\mu }}U\left(t\right):={\int }_0^1{D}_{*}^{\beta }U\left(t\right)d\mathit{\mu }\left(\beta \right)=\mathcal{F}\left(A\right)U\left(t\right)+H\left(t\right),t>0,$$

(3)

subject to the initial condition

$$U\left(0\right)=\Phi ,$$

(4)

where $\mathit{\mu }=({\mu }_1,\dots ,{\mu }_m)$ is a family of bounded Borel measures on $[0,1]$, $A:X\to X$ is a closed operator with domain $\mathfrak{D}\left(A\right)$ in a Banach space X, $\mathcal{F}\left(A\right):X^m\to X^m$ is an operator matrix, $\Phi \in X^m$, and $H:[0,\infty )\to X^m$. Throughout this section, we assume that the initial data $\Phi$ and the forcing term H possess sufficient regularity to ensure that all expressions appearing below are well defined. We associate the matrix-valued operator $\mathcal{F}\left(A\right)$ with the matrix-valued function $\mathcal{F}\left(z\right),z\in \mathbb{C},$ which serves as a symbol of the operator. In componentwise form, system (3) can be written as

$${\int }_0^1{D}_{*}^{\beta }{u}_i\left(t\right)d{\mu }_i\left(\beta \right)={\displaystyle \sum _{j=1}^m}{f}_{ij}\left(A\right)u_j\left(t\right)+h_i\left(t\right),t>0,i=1,\dots ,m.$$

In the particular case ${\mu }_j={\delta }_{{\beta }_j}$, where ${\delta }_a$ denotes the Dirac measure concentrated at a, the system reduces to the incommensurate fractional system with vector order $\mathcal{B}=({\beta }_1,\dots ,{\beta }_m)$.

Throughout this paper, our objective is to derive formal representations of the solution operators. We do not address the general questions of existence and uniqueness of solutions. Instead, we assume that the initial datum $\Phi$ and the forcing term H belong to a class of admissible data for which the Laplace transform method is applicable and all operator expressions arising below are well defined. In particular, whenever resolvent-type expressions of the form ${(\psi \left(s\right)I-\mathcal{F}\left(z\right))}^{-1}\Phi$ occur, it is assumed that $\Phi$ belongs to the corresponding admissible subspace so that these expressions exist and depend analytically on s in the domain under consideration. Under these assumptions, the formal manipulations involving Laplace transforms, contour deformations, and inverse transforms are justified, and the resulting formulas provide representations of the associated solution operator.

We first consider the uniformly distributed-order case, meaning that all the components of the measure $\mathit{\mu }$ are equal, that is ${\mu }_1=\dots ={\mu }_m=\mu .$ We note that for convenience $\psi \left(r{e}^{\pm i\pi }\right)$ we write as

$$\psi \left(r{e}^{\pm i\pi }\right)={\int }_0^1{r}^{\beta }{e}^{\pm i\pi \beta }d\mu \left(\beta \right)=A\left(r\right)\pm iB\left(r\right)$$

where

$$A\left(r\right)={\int }_0^1{r}^{\beta }cos\left(\pi \beta \right)d\mu \left(\beta \right)\mathrm{and}B\left(r\right)={\int }_0^1{r}^{\beta }sin\left(\pi \beta \right)d\mu \left(\beta \right).$$

Theorem 1.

Let $\mathit{\mu }=(\mu ,\dots ,\mu )$ with μ a bounded Borel measure on $[0,1].$ Then the solution operator $S(t,A):X^m\to X^m$ corresponding to problem (3)-(4) admits the series representation

$$S(t,A)=I+{\displaystyle \sum _{k=1}^{\infty }}{\left[\mathcal{F}\left(A\right)\right]}^k{\mathsf{\Psi }}_k\left(t\right),$$

(5)

where

$${\mathsf{\Psi }}_k\left(t\right)={\mathcal{L}}^{-1}\left\{{\displaystyle \frac{1}{s\psi {\left(s\right)}^k}}\right\}\left(t\right).$$

Moreover, $S(t,A)$ admits the integral representation

$$S(t,A)={\int }_0^{\infty }{e}^{-rt}\mathit{K}(r,\mathcal{F}\left(A\right))dr,$$

(6)

where

$$\mathit{K}(r,\mathcal{F}\left(A\right))={\displaystyle \frac{1}{\pi r}}Im\left[\psi \left(r{e}^{i\pi }\right){\left(\psi \left(r{e}^{i\pi }\right)I-\mathcal{F}\left(A\right)\right)}^{-1}\right].$$

(7)

Proof. 

To construct the solution operator, we consider the homogeneous problem $H\left(t\right)=0$. Applying the Laplace transform $\mathcal{L}\left\{U\right(t\left)\right\}=\mathit{U}\left(s\right)$, and formally replacing A by a parameter $z\in G$, and using the Laplace transform formula, we obtain

$${\int }_0^1\left(s^{\alpha }\mathit{U}\left(s\right)-s^{\alpha -1}\Phi \right)d\mu \left(\alpha \right)=\mathcal{F}\left(z\right)\mathit{U}\left(s\right),$$

where $\mathcal{F}\left(z\right)$ is the symbol of the matrix-valued operator $\mathcal{F}\left(A\right)$. Since $\mu$ is a bounded measure, we may interchange integration and the Laplace transform, yielding

$$\psi \left(s\right)\mathit{U}\left(s\right)-s^{-1}\psi \left(s\right)\Phi =\mathcal{F}\left(z\right)\mathit{U}\left(s\right),$$

where $\psi \left(s\right)$ is defined in (1). Rearranging gives

$$(\psi \left(s\right)I-\mathcal{F}\left(z\right))\mathit{U}\left(s\right)=s^{-1}\psi \left(s\right)\Phi ,$$

and hence

$$\mathit{U}\left(s\right)=\widehat{S}(s,z)\Phi ,\widehat{S}(s,z)={\displaystyle \frac{\psi \left(s\right)}{s}}{(\psi \left(s\right)I-\mathcal{F}\left(z\right))}^{-1}.$$

(8)

To derive the series representation, we write

$${(\psi \left(s\right)I-\mathcal{F}\left(z\right))}^{-1}={\displaystyle \frac{1}{\psi \left(s\right)}}{\left(I-{\displaystyle \frac{\mathcal{F}\left(z\right)}{\psi \left(s\right)}}\right)}^{-1}.$$

For s with sufficiently large real part $\Re \left(s\right)$, the Neumann series converges yielding

$${\left(I-{\displaystyle \frac{\mathcal{F}\left(z\right)}{\psi \left(s\right)}}\right)}^{-1}={\displaystyle \sum _{k=0}^{\infty }}{\left({\displaystyle \frac{\mathcal{F}\left(z\right)}{\psi \left(s\right)}}\right)}^k.$$

Substituting this into $\widehat{S}\left(s\right)$ in (8) gives

$$\widehat{S}(s,z)={\displaystyle \frac{1}{s}}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \frac{{\left[\mathcal{F}\left(z\right)\right]}^k}{\psi {\left(s\right)}^k}},\Re \left(s\right)>\sigma ,$$

where $\sigma >0$ is sufficiently large. Taking the inverse Laplace transform term-wise, we obtain

$$S(t,z)=I+{\displaystyle \sum _{k=1}^{\infty }}{\left[\mathcal{F}\left(z\right)\right]}^k{\mathsf{\Psi }}_k\left(t\right),t\ge 0,z\in \mathbb{C},$$

where

$${\mathsf{\Psi }}_k\left(t\right)={\mathcal{L}}^{-1}\left\{{\displaystyle \frac{1}{s\psi {\left(s\right)}^k}}\right\}\left(t\right).$$

Replacing the symbol $\mathcal{F}\left(z\right)$ by the operator $\mathcal{F}\left(A\right)$ yields representation (5).

To derive the integral representation of the solution operator, we compute the inverse Laplace transform of the symbol $\widehat{S}(s,z)$ defined in (8). Thus,

$$S(t,z)={\displaystyle \frac{1}{2\pi i}}{\int }_{\mathsf{\Gamma }}{e}^{st}{\displaystyle \frac{\psi \left(s\right)}{s}}{\left(\psi \left(s\right)I-\mathcal{F}\left(z\right)\right)}^{-1}ds,$$

where $\mathsf{\Gamma }=\{s\in \mathbb{C}:\Re (s)=\sigma \}$ is the Bromwich contour. Since $\psi \left(s\right)$ contains fractional powers, it is analytic in $\mathbb{C}\setminus (-\infty ,0]$ and has a branch point at $s=0$. For admissible data, we assume that the resolvent expression ${(\psi \left(s\right)I-\mathcal{F}\left(z\right))}^{-1}$ is analytic in $\mathbb{C}\setminus (-\infty ,0]$. Consequently, the integrand

$$e^{st}{\displaystyle \frac{\psi \left(s\right)}{s}}{\left(\psi \left(s\right)I-\mathcal{F}\left(z\right)\right)}^{-1}$$

is analytic in the slit plane $\mathbb{C}\setminus (-\infty ,0]$. Hence, by Cauchy’s theorem, the Bromwich contour may be deformed to the sectorial Hankel contour.

$${\mathsf{\Gamma }}_{R,\epsilon ,\theta }=\{r{e}^{-i\theta }:\epsilon \le r\le R\}\cup C_R\cup \{r{e}^{i\theta }:\epsilon \le r\le R\}\cup C_{\epsilon },$$

where

$$C_R=\{R{e}^{i\phi }:\theta \le \left|\phi \right|\le \pi \},C_{\epsilon }=\{\epsilon e^{i\phi }:-\pi \le \phi \le \pi \}.$$

Since $\Re \left(s\right)\le Rcos\theta <0$ on $C_R$, we have

$$|e^{st}|\le e^{-cRt},c=-cos\theta >0.$$

Using the resolvent estimate

$${\parallel (\psi \left(s\right)I-\mathcal{F}\left(z\right))}^{-1}\parallel \le {\displaystyle \frac{C}{\left|\psi \right(s\left)\right|}},$$

it follows that

$$∥e^{st}\psi \left(s\right){(\psi \left(s\right)I-\mathcal{F}\left(z\right))}^{-1}{\displaystyle \frac{1}{s}}∥\le {\displaystyle \frac{C{e}^{-cRt}}{R}},s\in C_R.$$

Hence

$$∥{\int }_{C_R}{e}^{st}\psi \left(s\right){(\psi \left(s\right)I-\mathcal{F}\left(z\right))}^{-1}{\displaystyle \frac{ds}{s}}∥\le \pi C{e}^{-cRt}\to 0,R\to \infty .$$

Therefore the contribution from the large arc vanishes. Letting $\epsilon \to 0$ and subsequently $\theta \to {\pi }^{-}$, Cauchy’s theorem yields that the contour integral reduces to the jump across the branch cut $(-\infty ,0]$. Parameterizing the contour by $s=r{e}^{\pm i\pi }$, $r>0$, we obtain

$$\begin{array}{cc}\hfill S(t,z)& ={\displaystyle \frac{1}{2\pi i}}{\int }_0^{\infty }{e}^{-rt}[\psi \left(r{e}^{i\pi }\right){\left(\psi \left(r{e}^{i\pi }\right)I-\mathcal{F}\left(z\right)\right)}^{-1}\hfill \\ \hfill & -\psi \left(r{e}^{-i\pi }\right){\left(\psi \left(r{e}^{-i\pi }\right)I-\mathcal{F}\left(z\right)\right)}^{-1}]{\displaystyle \frac{dr}{r}}.\hfill \end{array}$$

Since $\psi \left(r{e}^{-i\pi }\right)=\overline{\psi \left(r{e}^{i\pi }\right)}$, and the corresponding matrix-valued resolvents are conjugate, we may use the identity $Z-\overline{Z}=2iIm\left(Z\right)$ to obtain

$$S(t,z)={\int }_0^{\infty }{e}^{-rt}\left[{\displaystyle \frac{1}{\pi r}}Im\left(\psi \left(r{e}^{i\pi }\right){\left(\psi \left(r{e}^{i\pi }\right)I-\mathcal{F}\left(z\right)\right)}^{-1}\right)\right]dr.$$

Replacing the symbol $\mathcal{F}\left(z\right)$ by the operator $\mathcal{F}\left(A\right)$ yields the representation (6) with kernel

$$\mathit{K}(r,\mathcal{F}\left(A\right))={\displaystyle \frac{1}{\pi r}}Im\left(\psi \left(r{e}^{i\pi }\right){\left(\psi \left(r{e}^{i\pi }\right)I-\mathcal{F}\left(A\right)\right)}^{-1}\right).$$

This completes the proof. □

Corollary 1.

Let the measure be concentrated at a single order $\beta \in (0,1)$, i.e., $\mu ={\delta }_{\beta }$. Then the solution operator has the representation

$$S(t,A)=E_{\beta }\left(\mathcal{F}\left(A\right)t^{\beta }\right).$$

Proof. 

In this particular case the system reduces to the commensurate fractional system

$$D_{*}^{\beta }U\left(t\right)=\mathcal{F}\left(A\right)U\left(t\right).$$

It is easy to see that, in this case $\psi \left(s\right)=s^{\beta }$ and the solution operator admits the integral representation

$$S(t,A)={\int }_0^{\infty }{e}^{-rt}\mathit{K}(r,\mathcal{F}\left(A\right))dr.$$

Without loss of generality, we can sssume that the symbol $\mathcal{F}\left(z\right)$ is diagonalizable, i.e.,

$$\mathcal{F}\left(z\right)=P\mathsf{\Lambda }{P}^{-1},\mathsf{\Lambda }=\mathrm{diag}({\lambda }_1\left(z\right),\dots ,{\lambda }_m\left(z\right)),$$

where $P_j\left(z\right)$ are the spectral projectors of $\mathcal{F}\left(z\right)$. Then

$$S(t,z)=P\left(z\right)\mathrm{diag}\left({\int }_0^{\infty }{e}^{-rt}\mathcal{K}(r,{\lambda }_j\left(z\right))dr\right)P^{-1}\left(z\right),$$

where

$$\mathcal{K}(r,{\lambda }_j\left(z\right))={\displaystyle \frac{1}{\pi r}}Im\left({\displaystyle \frac{r^{\beta }{e}^{i\pi \beta }}{r^{\beta }{e}^{i\pi \beta }-{\lambda }_j}}\right)={\displaystyle \frac{1}{\pi }}{\displaystyle \frac{-{\lambda }_j\left(z\right)r^{\beta -1}sin\left(\pi \beta \right)}{r^{2\beta }-2{\lambda }_j\left(z\right)r^{\beta }cos\left(\pi \beta \right)+{\lambda }_j^2\left(z\right)}}.$$

In accordance with the integral representation of the Mittag-Leffler function (Lemma 1) the latter yields

$${\int }_0^{\infty }{e}^{-rt}\mathcal{K}(r,{\lambda }_j\left(z\right))dr=E_{\beta }\left({\lambda }_j\left(z\right)t^{\beta }\right),$$

and therefore

$$S(t,z)=P\left(z\right)\mathrm{diag}\left(E_{\beta }\left({\lambda }_j\left(z\right)t^{\beta }\right)\right)P^{-1}\left(z\right).$$

Since $E_{\beta }\left(0\right)=1$, it follows that $S(0,z)=I$. Moreover, by the spectral calculus,

$$S(t,z)={\displaystyle \sum _{j=1}^m}P\left(z\right)E_{\beta }\left({\lambda }_j\left(z\right)t^{\beta }\right)P_j\left(z\right)=E_{\beta }\left(\mathcal{F}\left(z\right)t^{\beta }\right),$$

implying that the solution operator is $S(t,A)=E_{\beta }\left(\mathcal{F}\left(A\right)t^{\beta }\right).$ This confirms that the solution operator reduces to the well-known representation via the classical Mittag-Leffler function valid for commensurate systems. □

Corollary 2.

Let the measure consist of two discrete components,

$$\mu ={\delta }_{\beta }+a{\delta }_{\gamma },0<\gamma <\beta <1,a>0.$$

Then the solution operator admits the integral representation

$$S(t,A)={\displaystyle \sum _{k=0}^{\infty }}\left[t^{k\beta }{E}_{\beta -\gamma ,k\beta +1}^k\left(-a{t}^{\beta -\gamma }\right)\right]{\left[\mathcal{F}\left(A\right)\right]}^k.$$

(9)

Proof. 

Notice that in this particular case $\psi \left(s\right)=s^{\beta }+a{s}^{\gamma }.$ Again for simplicity, assume that the symbol $\mathcal{F}\left(z\right)$ is diagonalizable,

$$\mathcal{F}\left(z\right)=P\left(z\right)\mathsf{\Lambda }\left(z\right)P^{-1}\left(z\right),\mathsf{\Lambda }\left(z\right)=\mathrm{diag}({\lambda }_1\left(z\right),\dots ,{\lambda }_m\left(z\right)).$$

Then

$$S(t,z)=P\left(z\right)\mathrm{diag}\left(s({\lambda }_j\left(z\right),t)\right)P^{-1}\left(z\right),$$

where each scalar component admits the series representation

$$s(\lambda \left(z\right),t)={\displaystyle \sum _{k=0}^{\infty }}{\left[\lambda \left(z\right)\right]}^k{t}^{k\beta }{E}_{\beta -\gamma ,k\beta +1}^k\left(-a{t}^{\beta -\gamma }\right).$$

(10)

Indeed, this representation follows from the Laplace transform formula

$$\widehat{S}(s,z)={\displaystyle \frac{1}{s}}{\left(I-\psi {\left(s\right)}^{-1}\mathcal{F}\left(z\right)\right)}^{-1}={\displaystyle \sum _{k=0}^{\infty }}{\left[\mathcal{F}\left(z\right)\right]}^k{\displaystyle \frac{1}{s\psi {\left(s\right)}^k}}.$$

For $\psi \left(s\right)=s^{\beta }+a{s}^{\gamma }$, we rewrite

$${\displaystyle \frac{1}{s\psi {\left(s\right)}^k}}=s^{-(k\beta +1)}{\left(1+a{s}^{-(\beta -\gamma )}\right)}^{-k}.$$

Applying the inverse Laplace transform formula

$${\mathcal{L}}^{-1}\left\{s^{-\rho }{(1-w{s}^{-\alpha })}^{-\delta }\right\}=t^{\rho -1}{E}_{\alpha ,\rho }^{\delta }\left(w{t}^{\alpha }\right),$$

we obtain

$${\mathsf{\Psi }}_k\left(t\right)=t^{k\beta }{E}_{\beta -\gamma ,k\beta +1}^k\left(-a{t}^{\beta -\gamma }\right),$$

which yields (10). Thus, the solution operator admits the representation (9). □

Remark 1.

It follows from Theorem 1 that the integral representation of the solution operator under the condition of Corollary 2 takes the form

$$S(t,A)={\int }_0^{\infty }{e}^{-rt}\mathit{K}(r,\mathcal{F}\left(A\right))dr,$$

where

$$\mathit{K}(r,\mathcal{F}\left(A\right))={\displaystyle \frac{1}{\pi r}}Im\left[\left(r^{\beta }{e}^{i\pi \beta }+a{r}^{\gamma }{e}^{i\pi \gamma }\right){\left(\psi \left(r{e}^{i\pi }\right)I-\mathcal{F}\left(A\right)\right)}^{-1}\right].$$

(11)

Example 1.

Consider an example of the uniformly distributed-order case $\mu \left(d\beta \right)=d\beta$ on $[0,1]$. The choice of a uniform measure corresponds to a system in which no single fractional order dominates the dynamics. Instead, all fractional exponents contribute equally, leading to a model with a continuous spectrum of memory effects. In contrast to single-order fractional models, this setting captures a broad range of relaxation mechanisms that arise in highly disordered media. The distributed-order operator in this case takes the form

$$D_{\mathit{\mu }}{U}_{(}t)={\int }_0^1{D}_{*}^{\beta }{U}_{(}t)d\beta .$$

(12)

It admits the convolution representation

$$D_{\mathit{\mu }}{U}_{(}t)={\int }_0^t\varphi (t-\tau )U_i^{\prime }\left(\tau \right)d\tau ,\varphi \left(t\right)={\int }_0^1{\displaystyle \frac{t^{-\beta }}{\mathsf{\Gamma }(1-\beta )}}d\beta .$$

(13)

The memory kernel satisfies the asymptotic behavior

$$\varphi \left(t\right)\sim {\displaystyle \frac{1}{t{(lnt)}^2}},t\to \infty .$$

(14)

The corresponding symbol function $\psi \left(s\right)$ is given by

$$\psi \left(s\right)={\int }_0^1{s}^{\beta }d\beta ={\displaystyle \frac{s-1}{lns}},\Re \left(s\right)>0,\psi \left(1\right)=1.$$

(15)

Then the associated auxiliary functions $A\left(r\right)$ and $B\left(r\right)$ are defined by

$$A\left(r\right)={\displaystyle \frac{-(r+1)lnr}{{(lnr)}^2+{\pi }^2}},B\left(r\right)={\displaystyle \frac{(r+1)\pi }{{(lnr)}^2+{\pi }^2}}.$$

(16)

Indeed, setting $s=r{e}^{\pm i\pi }$ and $ln\left(r{e}^{\pm i\pi }\right)=lnr\pm i\pi$, one obtains

$$\psi \left(r{e}^{\pm i\pi }\right)=-{\displaystyle \frac{(r+1)(lnr\mp i\pi )}{{(lnr)}^2+{\pi }^2}},$$

(17)

which yields (16). Applying Theorem 1, the solution operator admits the representation

$$S(t,A)=I+{\displaystyle \sum _{k=1}^{\infty }}{\left[\mathcal{F}\left(A\right)\right]}^k{\mathsf{\Psi }}_k\left(t\right),{\mathsf{\Psi }}_k\left(t\right)={\mathcal{L}}^{-1}\left\{{\displaystyle \frac{{(lns)}^k}{s{(s-1)}^k}}\right\}\left(t\right),$$

(18)

or equivalently

$$S(t,A)={\int }_0^{\infty }{e}^{-rt}\mathit{K}(r,\mathcal{F}\left(A\right))dr,$$

(19)

with

$$\mathit{K}(r,\mathcal{F}\left(A\right))={\displaystyle \frac{1}{\pi r}}Im\left(\psi \left(r{e}^{i\pi }\right){\left(\psi \left(r{e}^{i\pi }\right)I-\mathcal{F}\left(A\right)\right)}^{-1}\right).$$

(20)

The integral representation (19) shows that the long-time behavior of the solution operator is determined by the behavior of $\mathit{K}(r,\mathcal{F}(A\left)\right)$ as $r\to 0^{+}$. Indeed, it follows from (1) that

$$\psi \left(r{e}^{i\pi }\right)=-{\displaystyle \frac{lnr-i\pi }{{(lnr)}^2+{\pi }^2}}=O\left({\displaystyle \frac{1}{|lnr|}}\right),r\to 0+.$$

Since $\psi \left(r\right)\to 0$, and assuming $0\notin \sigma \left(\mathcal{F}\right(A\left)\right)$, in terms of symbols $\mathcal{F}\left(z\right),$ we have

$${(\psi \left(r\right)I-\mathcal{F}\left(z\right))}^{-1}\to {(-\mathcal{F}\left(z\right))}^{-1},z\in \sigma \left(\mathcal{F}\left(A\right)\right).$$

Therefore,

$$\mathit{K}(r,\mathcal{F}\left(z\right))={\displaystyle \frac{1}{\pi r}}Im\left(\psi \left(r{e}^{i\pi }\right){(\psi \left(r{e}^{i\pi }\right)I-\mathcal{F}\left(z\right))}^{-1}\right)\sim {\displaystyle \frac{1}{\pi r}}\psi \left(r{e}^{i\pi }\right){(-\mathcal{F}\left(z\right))}^{-1}.$$

Consequently, the kernel admits the representation

$$\mathit{K}(r,\mathcal{F}\left(z\right))={\displaystyle \frac{1}{r|lnr|}}R(r,z),$$

where $R(r,z)$ remains bounded as $r\to 0^{+}$ for each fixed z. As a consequence, by Karamata-type Tauberian arguments for Laplace transforms (Lemma 3) with slowly varying corrections, the solution operator satisfies the logarithmic decay asymptote

$$\parallel S(t,z)\parallel =O\left({\displaystyle \frac{1}{lnt}}\right),t\to \infty .$$

This establishes that the evolution generated by $S(t,A)$ exhibits ultra-slow relaxation, in the sense that its decay is slower than any algebraic rate $t^{-\beta }$, $\beta \in (0,1)$, and is governed entirely by logarithmic corrections arising from the distributed-order structure with measure $\mu \left(d\beta \right)=d\beta .$ In particular, the absence of a dominant fractional exponent leads to a cumulative memory effect across all orders in $[0,1]$, which manifests in a logarithmically moderated long-time decay characteristic of ultra-slow diffusion.

Corollary 2 illustrates a specific case of multi-term systems where the equations share a uniform memory structure (the same two fractional derivatives). In the following section, we extend this result to the general multi-term case, where each equation is governed by its own distinct multi-term fractional operator.

4. Solution Operators Associated with Heterogeneous Multi-term Systems

In the previous section, we examined the uniform distributed order case, where every equation in the system shared the identical memory measure $\mu$. This homogeneity allowed the memory operator to commute with the coupling matrix $\mathcal{F}$, resulting in a solution operator $S(t,A)$ that could be expressed as a simple matrix function of the form

$$I+\sum {\mathsf{\Psi }}_k\left(t\right){\mathcal{F}}^k\left(A\right).$$

In that setting, the spatial coupling and the temporal relaxation were perfectly decoupled, meaning the system’s memory was a global property shared by all components.

In contrast, the system considered in this section represents a heterogeneous multi-term case. Here, each equation i is governed by its own distinct set of fractional orders ${\beta }_i$ and coefficients $a_{ij}$, resulting in unique memory profiles $P_i\left(s\right)$ for each component. Because the left-hand side operators now vary from one equation to another, they no longer commute with the coupling operators $f_{ij}\left(A\right)$. This structural change implies that the evolution of any single component $u_i\left(t\right)$ is filtered through a unique “path” of differing memory structures as it interacts with other components, leading to a path-dependent representation of the solution operator.

The memory-decoupled heterogeneous multi-term fractional system has the form

$$D^{{\beta }_i}{u}_i\left(t\right)+\sum _{j\ne i}{a}_{ij}{D}^{{\beta }_{ij}}{u}_i\left(t\right)={\displaystyle \sum _{j=1}^m}{f}_{ij}\left(A\right)u_j\left(t\right),u_i\left(0\right)={\phi }_i\in X,$$

(21)

where $0<{\beta }_{ij}\le {\beta }_i\le 1,i,j=1,\dots ,m,$ and coefficients $a_{i,j},i\ne j,$ are arbitrary numbers.

Theorem 2.

The formal solution operator $S(t,A)$ for system (21) admits the series representation

$$S(t,A)=I+{\displaystyle \sum _{k=1}^{\infty }}{\mathcal{M}}_k(t,A),$$

(22)

where the entries of the k-th order interaction operator ${\mathcal{M}}_k(t,A)$ are

$${\left[{\mathcal{M}}_k(t,A)\right]}_{ij}={\displaystyle \sum _{i_1,\dots ,i_{k-1}=1}^m}\left(f_{i{i}_1}\left(A\right)\circ f_{i_1{i}_2}\left(A\right)\circ \dots \circ f_{i_{k-1}j}\left(A\right)\right){\psi }_{(i,i_1,\dots ,i_{k-1},j)}\left(t\right),$$

(23)

with the scalar weights ${\psi }_{(i,i_1,\dots ,i_{k-1},j)}\left(t\right)$ defined by

$${\psi }_{(i,i_1,\dots ,j)}\left(t\right)={\mathcal{L}}^{-1}\left\{{\displaystyle \frac{1}{s{P}_i\left(s\right)P_{i_1}\left(s\right)\dots P_j\left(s\right)}}\right\}\left(t\right),$$

(24)

and

$$P_n\left(s\right)=s^{{\beta }_n}+\sum _{\ell \ne n}{a}_{n\ell }{s}^{{\beta }_{n\ell }}.$$

The series (22) converges in the topology of X for each $t\ge 0$ and admissible datum $\Phi \in {\bigcap }_{k=1}^{\infty }\mathcal{D}\left({\mathcal{F}}^k\left(A\right)\right)$, provided that all operator compositions are well-defined.

Proof. 

Applying the Laplace transform to the i-th equation in (21), we obtain

$$\mathcal{L}\left\{D^{{\beta }_i}{u}_i\left(t\right)\right\}+\sum _{j\ne i}{a}_{ij}\mathcal{L}\left\{D^{{\beta }_{ij}}{u}_i\left(t\right)\right\}={\displaystyle \sum _{j=1}^m}{f}_{ij}\left(A\right){\widehat{U}}_j\left(s\right).$$

Using the formula (2) and defining

$$P_i\left(s\right)=s^{{\beta }_i}+\sum _{j\ne i}{a}_{ij}{s}^{{\beta }_{ij}},$$

we obtain

$$P_i\left(s\right){\widehat{U}}_i\left(s\right)-s^{-1}{P}_i\left(s\right){\phi }_i={\displaystyle \sum _{j=1}^m}{f}_{ij}\left(A\right){\widehat{U}}_j\left(s\right).$$

Further, let $\mathit{P}\left(s\right)=\mathrm{diag}(P_1\left(s\right),\dots ,P_m\left(s\right))$ and $\mathcal{F}\left(z\right)$ is the symbol of the operator $\mathcal{F}\left(A\right)$. Then

$$(\mathit{P}\left(s\right)I-\mathcal{F}\left(z\right))\widehat{U}\left(s\right)=s^{-1}\mathit{P}\left(s\right)\Phi .$$

Thus,

$$\widehat{U}\left(s\right)=s^{-1}I-\mathit{P}{\left(s\right)}^{-1}{\mathcal{F}\left(z\right))}^{-1}\Phi .$$

For s with $\Re \left(s\right)$ sufficiently large, $\mathit{P}{\left(s\right)}^{-1}$ is well-defined and ${\parallel \mathit{P}\left(s\right)}^{-1}\mathcal{F}\left(z\right)\parallel$ is sufficiently small, so the Neumann series converges

$$\widehat{S}(s,z)={\displaystyle \frac{1}{s}}{\displaystyle \sum _{k=0}^{\infty }}{\left(\mathit{P}{\left(s\right)}^{-1}\mathcal{F}\left(z\right)\right)}^k.$$

Since $\Phi \in {\bigcap }_{k=1}^{\infty }\mathcal{D}\left({\mathcal{F}}^k\right)$, all compositions are well-defined. The $(i,j)$-entry of the k-th term is

$${\left[{\widehat{S}}_k(s,z)\right]}_{ij}=\sum _{i_1,\dots ,i_{k-1}}{\displaystyle \frac{1}{s{P}_i\left(s\right)P_{i_1}\left(s\right)\dots P_j\left(s\right)}}(f_{i{i}_1}\left(z\right)\circ \dots \circ f_{i_{k-1}j}\left(z\right)).$$

(25)

The inverse Laplace transform defines the symbol of the solution operator $S(t,z)={\mathcal{L}}^{-1}\left\{\widehat{S}(s,z)\right\}.$ Termwise inversion is justified by uniform convergence for sufficiently large $\Re \left(s\right)$, yielding

$$S(t,z)=I+{\displaystyle \sum _{k=1}^{\infty }}{\mathcal{M}}_k(t,z).$$

Finally, for large $\left|s\right|$, one has $P_n\left(s\right)\sim s^{{\beta }_n}$, and therefore we obtain the following asymptotic relation

$${\displaystyle \frac{1}{s{P}_i\left(s\right)\dots P_j\left(s\right)}}\sim s^{-(1+\sum {\beta }_{i_n})}.$$

Hence,

$${\psi }_{(i,\dots ,j)}\left(t\right)\sim {\displaystyle \frac{t^{\sum {\beta }_{i_n}}}{\mathsf{\Gamma }(\sum {\beta }_{i_n}+1)}}.$$

This decay ensures convergence of the series (25) by the dominating convergent multivariate Mittag-Leffler series for each fixed t and admissible $\Phi .$ This complets the proof. □

Example 2.

To better feel the construction of the solution operator, consider a $2\times 2$ system with distinct memory polynomials $P_1\left(s\right)$ and $P_2\left(s\right)$ and coupling operators $f_{ij}\left(A\right)$. The solution operator in this case has the form

$$S(t,A)=I+{\displaystyle \sum _{k=1}^{\infty }}{\mathcal{M}}_k(t,A)$$

and exhibits a path-dependent weighting structure. To illustrate this, we examine the first two interaction terms. The first-order term is

$${\mathcal{M}}_1(t,A)=\left(\begin{array}{cc}{\psi }_{\left(1\right)}\left(t\right)& 0\\ 0& {\psi }_{\left(2\right)}\left(t\right)\end{array}\right)\mathcal{F}\left(A\right),{\psi }_{\left(i\right)}\left(t\right)={\mathcal{L}}^{-1}\left\{{\displaystyle \frac{1}{s{P}_i\left(s\right)}}\right\}.$$

(26)

The second-order term involves compositions of operators and mixed memory effects. For example, the $(1,2)$-entry is

$${\left[{\mathcal{M}}_2(t,A)\right]}_{12}=f_{11}\left(A\right)f_{12}\left(A\right){\mathcal{L}}^{-1}\left\{{\displaystyle \frac{1}{s{P}_1{\left(s\right)}^2}}\right\}+f_{12}\left(A\right)f_{22}\left(A\right){\mathcal{L}}^{-1}\left\{{\displaystyle \frac{1}{s{P}_1\left(s\right)P_2\left(s\right)}}\right\}.$$

(27)

The full solution operator $S(t,A)={\left(S_{ij}(t,A)\right)}_{i,j=1}^2$ can be expressed entrywise. The diagonal term $S_{11}\left(t\right)$, representing internal dynamics and feedback loops, has the expansion

$$S_{11}(t,A)=I+f_{11}\left(A\right){\psi }_{(1,1)}\left(t\right)+\left(f_{11}^2\left(A\right){\psi }_{(1,1,1)}\left(t\right)+f_{12}\left(A\right)f_{21}\left(A\right){\psi }_{(1,2,1)}\left(t\right)\right)+\dots ,$$

(28)

where

$${\psi }_{(1,2,1)}\left(t\right)={\mathcal{L}}^{-1}\left\{{\displaystyle \frac{1}{s{P}_1{\left(s\right)}^2{P}_2\left(s\right)}}\right\}.$$

Similarly, the off-diagonal term $S_{12}(t,A)$, describing the influence of ${\phi }_2$ on $u_1$, is

$$S_{12}(t,A)=f_{12}\left(A\right){\psi }_{(1,2)}\left(t\right)+\left(f_{11}\left(A\right)f_{12}\left(A\right){\psi }_{(1,1,2)}\left(t\right)+f_{12}\left(A\right)f_{22}\left(A\right){\psi }_{(1,2,2)}\left(t\right)\right)+\dots ,$$

(29)

with

$$\begin{array}{cc}\hfill {\psi }_{(1,2)}\left(t\right)& ={\mathcal{L}}^{-1}\left\{{\displaystyle \frac{1}{s{P}_1\left(s\right)P_2\left(s\right)}}\right\},{\psi }_{(1,1,2)}\left(t\right)={\mathcal{L}}^{-1}\left\{{\displaystyle \frac{1}{s{P}_1{\left(s\right)}^2{P}_2\left(s\right)}}\right\},\hfill \\ \hfill {\psi }_{(1,2,2)}\left(t\right)& ={\mathcal{L}}^{-1}\left\{{\displaystyle \frac{1}{s{P}_1\left(s\right)P_2{\left(s\right)}^2}}\right\}.\hfill \end{array}$$

In general, each entry admits a path expansion of the form

$$S_{ij}(t,A)={\delta }_{ij}I+{\displaystyle \sum _{k=1}^{\infty }}\sum _{i_1,\dots ,i_{k-1}\in \{1,2\}}(f_{i{i}_1}\left(A\right)\dots f_{i_{k-1}j}\left(A\right)){\psi }_{(i,i_1,\dots ,j)}\left(t\right),$$

(30)

where

$${\psi }_{(n_1,\dots ,n_k)}\left(t\right)={\mathcal{L}}^{-1}\left\{{\displaystyle \frac{1}{s{\prod }_{j=1}^k{P}_{n_j}\left(s\right)}}\right\},$$

(31)

with indices $n_j\in \{1,2\}$, so that each factor $P_{n_j}\left(s\right)$ is either $P_1\left(s\right)$ or $P_2\left(s\right)$.

This example shows that in heterogeneous systems, the temporal evolution is not governed by a single global memory kernel, but is determined by the sequence of memory scales encountered along each interaction path.

Example 3.

If the memory structures become uniform, so that $P_1\left(s\right)=P_2\left(s\right)=P\left(s\right),$ the heterogeneous solution operator $S\left(t\right)$ reduces to a homogeneous representation. In this case, the path-dependent weights no longer depend on the specific index sequence and simplify to

$${\psi }_{(i_1,i_2,\dots ,i_k)}\left(t\right)={\mathsf{\Psi }}_k\left(t\right)={\mathcal{L}}^{-1}\left\{{\displaystyle \frac{1}{s{\left[P\left(s\right)\right]}^k}}\right\},$$

(32)

which is a scalar function depending only on the path length k. Under this condition, the $(1,2)$-entry of the solution operator becomes

$$S_{12}(t,A)={\displaystyle \sum _{k=1}^{\infty }}\left(\sum _{i_1,\dots ,i_{k-1}\in \{1,2\}}{f}_{1i_1}\left(A\right)\circ f_{i_1{i}_2}\left(A\right)\circ \dots \circ f_{i_{k-1}2}\left(A\right)\right){\mathsf{\Psi }}_k\left(t\right).$$

(33)

The inner sum represents all length-k paths from index 1 to 2, and coincides with the $(1,2)$-entry of the k-th power of the operator matrix $\mathcal{F}\left(A\right)=\left[f_{ij}\left(A\right)\right]$. Hence,

$$S_{12}(t,A)={\displaystyle \sum _{k=1}^{\infty }}{\left[\mathcal{F}{\left(A\right)}^k\right]}_{12}{\mathsf{\Psi }}_k\left(t\right).$$

(34)

Similarly, the diagonal entries satisfy

$$S_{ii}(t,A)=I+{\displaystyle \sum _{k=1}^{\infty }}{\left[\mathcal{F}{\left(A\right)}^k\right]}_{ii}{\mathsf{\Psi }}_k\left(t\right),i=1,2.$$

(35)

Therefore, the full solution operator factorizes as

$$S(t,A)=I+{\displaystyle \sum _{k=1}^{\infty }}\mathcal{F}{\left(A\right)}^k{\mathsf{\Psi }}_k\left(t\right),$$

(36)

recovering the homogeneous matrix-valued expansion in which the temporal evolution is governed by a single family of scalar kernels ${\left\{{\mathsf{\Psi }}_k\left(t\right)\right\}}_{k\ge 1}$.

Systems involving multi-term fractional derivatives can be further generalized to the heterogeneous case in which each component is governed by a distributed-order operator. In this setting, the j-th component is associated with a finite Borel measure ${\mu }_j$ on $[0,1]$, and the system takes the form

$${\int }_0^1{D}_{*}^{\beta }{u}_j\left(t\right){\mu }_j\left(d\beta \right)={\displaystyle \sum _{k=1}^m}{f}_{jk}\left(A\right)u_k\left(t\right),j=1,\dots ,m,$$

(37)

with initial conditions $u_j\left(0\right)={\phi }_j\in X.$ This formulation provides a unified framework encompassing:

• single-order models (${\mu }_j={\delta }_{\beta }$),
• multi-term models (${\mu }_j={\sum }_k{a}_{jk}{\delta }_{{\beta }_k}$),
• continuous-order models ($d{\mu }_j={\omega }_j\left(\beta \right)d\beta$),
• and the homogeneous case (${\mu }_j=\mu$ for all j).

Applying the Laplace transform to the system and taking into account the initial data, we obtain

$${\psi }_j\left(s\right){\widehat{U}}_j\left(s\right)={\displaystyle \sum _{k=1}^m}{f}_{jk}\left(A\right){\widehat{U}}_k\left(s\right)+G_j\left(s\right),$$

where $G_j\left(s\right)$ collects the contributions of the initial conditions and

$${\psi }_j\left(s\right)={\int }_0^1{s}^{\beta }{\mu }_j\left(d\beta \right).$$

The obtained system can be written in matrix form as

$$Q(s,A)\widehat{U}\left(s\right)=G\left(s\right),Q(s,A)=\mathrm{diag}({\psi }_1\left(s\right),\dots ,{\psi }_m\left(s\right))-\mathcal{F}\left(A\right).$$

Whenever the resolvent exists and satisfies a suitable smallness condition, assuming that

$${\parallel \mathcal{F}\left(A\right)\psi \left(s\right)}^{-1}{\parallel <1,\psi \left(s\right)}^{-1}=\mathrm{diag}({\psi }_1\left(s\right),\dots ,{\psi }_m\left(s\right))$$

we obtain the Neumann series expansion

$$Q{(s,A)}^{-1}={(\psi \left(s\right)-\mathcal{F}\left(A\right))}^{-1}=\psi {\left(s\right)}^{-1}{\displaystyle \sum _{k=0}^{\infty }}{\left(\mathcal{F}\left(A\right)\psi {\left(s\right)}^{-1}\right)}^k.$$

The Laplace symbol of the solution operator is therefore given by

$$\widehat{S}(s,A)={\displaystyle \frac{1}{s}}Q{(s,A)}^{-1}={\displaystyle \frac{1}{s}}\psi {\left(s\right)}^{-1}{\displaystyle \sum _{k=0}^{\infty }}{\left(\mathcal{F}\left(A\right)\psi {\left(s\right)}^{-1}\right)}^k.$$

Now applying the inverse Laplace transform, we obtain

$$S(t,A)={\displaystyle \sum _{k=0}^{\infty }}{K}_k(t,A),$$

where the operators $K_k(t,A)$ are defined as the inverse Laplace transforms of

$${\widehat{K}}_k(s,A)={\displaystyle \frac{1}{s}}\psi {\left(s\right)}^{-1}{\left(\mathcal{F}\left(A\right)\psi {\left(s\right)}^{-1}\right)}^k.$$

Each operator $K_k(t,A)$ admits a spectral representation of the form

$$K_k(t,A)={\int }_0^{\infty }{e}^{-rt}{\mathsf{\Psi }}_k(r,A)dr,$$

(38)

where ${\mathsf{\Psi }}_k(r,A)$ is the matrix-valued spectral density associated with the boundary values of the resolvent. Writing

$${\psi }_j\left(r{e}^{\pm i\pi }\right)=A_j\left(r\right)\pm i{B}_j\left(r\right),$$

where

$$A_j\left(r\right)={\int }_0^1{r}^{\beta }cos\left(\pi \beta \right){\mu }_j\left(d\beta \right),B_j\left(r\right)={\int }_0^1{r}^{\beta }sin\left(\pi \beta \right){\mu }_j\left(d\beta \right),$$

the spectral density is obtained from the jump across the branch cut, and has the form

$${\mathsf{\Psi }}_k(r,A)={\displaystyle \frac{1}{\pi r}}Im\left[{\left(\psi \left(r{e}^{i\pi }\right)\right)}^{-1}{\left(\mathcal{F}\left(A\right)\psi {\left(r{e}^{i\pi }\right)}^{-1}\right)}^k\right].$$

(39)

Thus, the following theorem holds.

Theorem 3.

Let ${\mu }_j,j=1,\dots ,m,$ be bounded Borel measures supported on $(0,1].$ Then the formal solution operator associated with system (37) admits the representation

$$S(t,A)=I+{\displaystyle \sum _{k=1}^{\infty }}{K}_k(t,A),$$

(40)

where operators $K_k(t,A),k=0,\dots ,$ are defined in (38) with kernel operators ${\mathsf{\Psi }}_k(r,A)$ determined in (39).

5. Solution Operators for Fully Coupled Systems

In the preceding sections, we investigated systems where each equation was governed by its own distinct multi-term fractional operator. While those cases allowed for heterogeneous memory across different components, the fractional differentiation was still applied row-wise, meaning the memory of each variable $u_i$ was essentially an internal property of that specific equation. The coupling between components occurred exclusively on the right-hand side through the operator matrix F. This structure preserved a diagonal dominance in the Laplace domain, which allowed the solution operator to be represented as a sum of path-dependent compositions weighted by scalar-valued memory scales.

The system considered in this section represents a further level of complexity, referred to as the fully coupled matrix-order case. In this setting, fractional differentiation is defined via a Hadamard product $A\circ {\mathfrak{D}}^{\mathcal{A}}$, where each interaction between components $u_j$ and $u_k$ is assigned its own fractional order ${\alpha }_{jk}$ and weight $a_{jk}$. This structure implies that the memory of the system is no longer localized within individual equations but is instead globally entangled across components through the left-hand side. We refer to this phenomenon as memory coupling. At the symbolic level, this leads to a matrix-valued fractional symbol $\mathit{B}\left(s\right)={\left(a_{jk}{s}^{{\alpha }_{jk}}\right)}_{j,k=1}^m,$ which is no longer diagonal. As a result, the analysis requires inversion of a fully coupled matrix-valued symbol, in contrast to the diagonal (componentwise) cases considered previously.

Consider the system

$$A\circ {\mathfrak{D}}_{*}^{\mathcal{A}}U\left(t\right)=\mathcal{F}\left(A\right)U\left(t\right),U\left(0\right)=\Phi \in X^m,$$

(41)

where $A=\left(a_{jk}\right)$ is a constant matrix, $\mathcal{A}=\left({\alpha }_{jk}\right)$ with $0<{\alpha }_{jk}\le 1$, and $\mathcal{F}=\left(f_{jk}\right)$ is an $m\times m$ matrix of closed, densely defined linear operators on X. In system (41) the $A\circ {\mathfrak{D}}^{\mathcal{A}}$ is the Hadamard entrywise matrix product of matrices $A=\left\{a_{ij}\right\}$ and ${\mathfrak{D}}^{\mathcal{A}}=\left\{D_{*}^{{\alpha }_{ij}}\right\}.$ The componentwise form of system (41) is

$${\displaystyle \sum _{j=1}^m}{a}_{ij}{D}_{*}^{{\alpha }_{ij}}{u}_j\left(t\right)={\displaystyle \sum _{j=1}^m}{f}_{ij}\left(A\right)u_j\left(t\right),t>0,i=1,\dots ,m.$$

with the initial conditions $u_i\left(0\right)={\phi }_i,i=1,\dots ,m.$

Theorem 4.

Assume that $det\left(A\right)\ne 0$. Then, the formal solution operator $S(t,A)$ corresponding to system (41) admits the operator-valued series representation

$$S(t,A)=I+{\displaystyle \sum _{k=1}^{\infty }}{\mathcal{M}}_k(t,A),$$

(42)

where

$${\mathcal{M}}_k(t,A)={\mathcal{L}}^{-1}\left\{{\displaystyle \frac{1}{s}}{\left(\mathit{B}{\left(s\right)}^{-1}\mathcal{F}\left(A\right)\right)}^k\right\},\mathit{B}\left(s\right)={\left(a_{jk}{s}^{{\alpha }_{jk}}\right)}_{j,k=1}^m.$$

(43)

Proof. 

Applying the Laplace transform $\mathcal{L}\left\{U\left(t\right)\right\}=\widehat{U}\left(s\right)$ to (41), and using the Laplace transform formula (2) for the Caputo–Djrbashian fractional derivative of order $0<\alpha <1$, we obtain, for each $j=1,\dots ,m$,

$${\displaystyle \sum _{k=1}^m}{a}_{jk}\left(s^{{\alpha }_{jk}}{\widehat{u}}_k\left(s\right)-s^{{\alpha }_{jk}-1}{\phi }_k\right)={\displaystyle \sum _{k=1}^m}{f}_{jk}\left(A\right){\widehat{u}}_k\left(s\right).$$

Introducing the matrix-valued Laplace symbol

$$\mathit{B}\left(s\right)={\left(a_{jk}{s}^{{\alpha }_{jk}}\right)}_{j,k=1}^m,$$

the system can be written in vector form on $X^m$ as

$$\mathit{B}\left(s\right)\widehat{U}\left(s\right)-s^{-1}\mathit{B}\left(s\right)\Phi =\mathcal{F}\left(A\right)\widehat{U}\left(s\right),$$

that is,

$$(\mathit{B}\left(s\right)-\mathcal{F}\left(A\right))\widehat{U}\left(s\right)=s^{-1}\mathit{B}\left(s\right)\Phi .$$

We now analyze the invertibility of $\mathit{B}\left(s\right)$. By the Leibniz formula, we have

$$det\left(\mathit{B}\left(s\right)\right)=\sum _{\sigma \in S_m}sgn\left(\sigma \right){\displaystyle \prod _{j=1}^m}{a}_{j,\sigma \left(j\right)}{s}^{{\alpha }_{j,\sigma \left(j\right)}}=\sum _{\sigma \in S_m}{c}_{\sigma }{s}^{E_{\sigma }},$$

where $S_m$ denotes the symmetric group of all permutations of $\{1,\dots ,m\}$, and

$$E_{\sigma }={\displaystyle \sum _{j=1}^m}{\alpha }_{j,\sigma \left(j\right)}.$$

Since $det\left(A\right)\ne 0$, there exists at least one permutation $\sigma \in S_m$ such that $c_{\sigma }\ne 0$. Let

$$E_{max}=\underset{\sigma \in S_m}{max}{E}_{\sigma }.$$

Then, as $\left|s\right|\to \infty$ in a sector avoiding the negative real axis, the terms corresponding to $E_{max}$ dominate asymptotically. Consequently, $det\left(\mathit{B}\right(s\left)\right)\ne 0$ for sufficiently large $\Re \left(s\right)$, and therefore $\mathit{B}\left(s\right)$ is invertible for $\Re \left(s\right)>s_0$, with sufficiently large $s_0>0$. For such s, we rewrite

$$\widehat{U}\left(s\right)=s^{-1}{(\mathit{B}\left(s\right)-\mathcal{F}\left(A\right))}^{-1}\mathit{B}\left(s\right)\Phi =s^{-1}{(I-\mathit{B}{\left(s\right)}^{-1}\mathcal{F}\left(A\right))}^{-1}\Phi .$$

Moreover, for sufficiently large $\Re \left(s\right)$, the operator $\mathit{B}{\left(s\right)}^{-1}\mathcal{F}$ is sufficiently small in the operator topology on $X^m$. Hence on admissible initial datum $\Phi$ the Neumann series expansion converges, and we obtain

$$\widehat{S}(s,A)\Phi =\left(s^{-1}I+{\displaystyle \sum _{k=1}^{\infty }}{s}^{-1}{\left(\mathit{B}{\left(s\right)}^{-1}\mathcal{F}\left(A\right)\right)}^k\right)\Phi .$$

The solution operator $S(t,A)$ is defined as the inverse Laplace transform of $\widehat{S}(s,A)$. Using the Neumann series representation of $\widehat{S}(s,A)$, we perform the inverse Laplace transform termwise and obtain

$$S(t,A)={\mathcal{L}}^{-1}\left\{s^{-1}I\right\}+{\displaystyle \sum _{k=1}^{\infty }}{\mathcal{L}}^{-1}\left\{{\displaystyle \frac{1}{s}}{\left(\mathit{B}{\left(s\right)}^{-1}\mathcal{F}\left(A\right)\right)}^k\right\}.$$

Since ${\mathcal{L}}^{-1}\left\{s^{-1}I\right\}=I,$ it follows that

$$S(t,A)=I+{\displaystyle \sum _{k=1}^{\infty }}{\mathcal{M}}_k(t,A).$$

Each term ${\mathcal{M}}_k(t,A)$ defines a well-defined operator on $X^m$. The convergence of the resulting series in the topology of $X^m$ follows from the convergence of the Neumann series in the Laplace domain, which permits term-by-term inversion of the Laplace transform. The proof is complete. □

The two examples below illustrate that the results obtained for memory decoupled systems considered in the previous sections can be obtained from Theorem 4.

Example 4.

The heterogeneous m-term system, where each component $u_i$ is governed by an m-term fractional operator $P_i\left(D\right)=D^{{\beta }_{i,i}}+{\sum }_{p=1}^{m-1}{a}_i^{\left(p\right)}{D}^{{\alpha }_i^{\left(p\right)}}$, can be recovered by expressing the memory operator as a sum of diagonal Hadamard products.

Consider the system defined by the sum of m entry-wise products

$$\left({\displaystyle \sum _{p=1}^m}{A}^{\left(p\right)}\circ {\mathfrak{D}}^{{\mathcal{A}}^{\left(p\right)}}\right)U\left(t\right)=\mathcal{F}U\left(t\right),$$

(44)

where $A^{\left(p\right)}$ and ${\mathcal{A}}^{\left(p\right)}$ are diagonal matrices. To match the heterogeneous m-term case, we let

$$\begin{array}{cc}\hfill A^{\left(m\right)}& =I,{\mathcal{A}}^{\left(m\right)}=\mathrm{diag}({\beta }_{1,1},\dots ,{\beta }_{m,m})\hfill \\ \hfill A^{\left(p\right)}& =\mathrm{diag}(a_1^{\left(p\right)},\dots ,a_m^{\left(p\right)}),{\mathcal{A}}^{\left(p\right)}=\mathrm{diag}({\alpha }_1^{\left(p\right)},\dots ,{\alpha }_m^{\left(p\right)})\hfill \end{array}$$

for $p=1,\dots ,m-1$. Assume that ${\beta }_{i,i}$ is the dominating order for each i (i.e., ${\beta }_{i,i}>{\alpha }_i^{\left(p\right)}$). The resulting matrix-valued symbol $\mathit{B}\left(s\right)$ is diagonal, that is

$$\mathit{B}\left(s\right)=\mathrm{diag}\left(s^{{\beta }_{1,1}}+{\displaystyle \sum _{p=1}^{m-1}}{a}_1^{\left(p\right)}{s}^{{\alpha }_1^{\left(p\right)}},\dots ,s^{{\beta }_{m,m}}+{\displaystyle \sum _{p=1}^{m-1}}{a}_m^{\left(p\right)}{s}^{{\alpha }_m^{\left(p\right)}}\right).$$

(45)

Identifying $P_i\left(s\right)=s^{{\beta }_{i,i}}+{\sum }_{p=1}^{m-1}{a}_i^{\left(p\right)}{s}^{{\alpha }_i^{\left(p\right)}}$, we obtain $\mathit{B}\left(s\right)=\mathrm{diag}\left(P_i\left(s\right)\right)$. Since the leading term $s^{{\beta }_{i,i}}$ dominates as $\mathrm{Re}\left(s\right)\to \infty$, the symbol $\mathit{B}\left(s\right)$ is invertible for large $s.$ The inverse symbol is $\mathit{B}{\left(s\right)}^{-1}=\mathrm{diag}(P_1{\left(s\right)}^{-1},\dots ,P_m{\left(s\right)}^{-1})$, and the solution operator $S\left(t\right)$ recovers the path-dependent series

$$S(t,A)=I+{\displaystyle \sum _{k=1}^{\infty }}{\mathcal{L}}^{-1}\left\{{\displaystyle \frac{1}{s}}{\left(\mathit{P}{\left(s\right)}^{-1}\mathcal{F}\left(A\right)\right)}^k\right\}.$$

(46)

This demonstrates that the m-term system is a diagonal restriction of a fully coupled framework where the dominating orders ensure well-posedness.

Example 5.

Let $A=I$ and let the matrix of fractional orders be diagonal, $\mathcal{A}=\mathrm{diag}({\beta }_1,\dots ,{\beta }_m),0<{\beta }_j\le 1.$ Then the symbol matrix reduces to $\mathit{B}\left(s\right)=\mathrm{diag}(s^{{\beta }_1},\dots ,s^{{\beta }_m}),$ and the system

$$\mathit{B}\left(s\right)\widehat{U}\left(s\right)-s^{-1}\mathit{B}\left(s\right)\Phi =\mathcal{F}\left(A\right)\widehat{U}\left(s\right)$$

becomes

$$\mathrm{diag}(s^{{\beta }_1},\dots ,s^{{\beta }_m})\widehat{U}\left(s\right)-s^{-1}\mathrm{diag}(s^{{\beta }_1},\dots ,s^{{\beta }_m})\Phi =\mathcal{F}\widehat{U}\left(s\right).$$

In particular, the resolvent structure simplifies to

$$\widehat{U}\left(s\right)=s^{-1}{\left(I-\mathrm{diag}{(s^{{\beta }_1},\dots ,s^{{\beta }_m})}^{-1}\mathcal{F}\left(A\right)\right)}^{-1}\Phi .$$

Since $\mathrm{diag}{(s^{{\beta }_1},\dots ,s^{{\beta }_m})}^{-1}=\mathrm{diag}(s^{-{\beta }_1},\dots ,s^{-{\beta }_m}),$ each component carries an independent fractional scaling in the Laplace domain. Expanding the resolvent yields the Neumann series representation

$$\widehat{S}(s,A)=s^{-1}{\displaystyle \sum _{k=0}^{\infty }}{\left(\mathrm{diag}(s^{-{\beta }_1},\dots ,s^{-{\beta }_m})\mathcal{F}\left(A\right)\right)}^k.$$

In the time domain, this leads to a decoupled family of convolution kernels. In particular, each entry of the solution operator takes the form

$$S_{ij}(t,A)={\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \sum _{i_1,\dots ,i_{k-1}=1}^m}{f}_{i{i}_1}\left(A\right)\dots f_{i_{k-1}j}\left(A\right){\psi }^{(i,i_1,\dots ,i_{k-1},j)}\left(t\right),$$

where

$${\psi }^{(i,i_1,\dots ,i_{k-1},j)}\left(t\right)={\mathcal{L}}^{-1}\left\{{\displaystyle \frac{1}{s{\prod }_{r=1}^{k-1}{s}^{{\beta }_{i_r}}}}\right\}\left(t\right)={\displaystyle \frac{t^{{\beta }_{i_1}+\dots +{\beta }_{i_{k-1}}}}{\mathsf{\Gamma }\left(1+{\beta }_{i_1}+\dots +{\beta }_{i_{k-1}}\right)}}.$$

Thus, in this diagonal fractional-order case, the memory structure is fully encoded by additive fractional exponents along interaction paths, and the solution operator reduces to a weighted operator series with explicit power-law kernels.

In the homogeneous case, where ${\beta }_i=\beta$ for all $i=1,\dots ,m$, the memory structure becomes uniform and the path-dependent exponents simplify. In this situation, each step contributes the same fractional scaling, and the scalar weights reduce to

$${\psi }^{\left(k\right)}\left(t\right)={\mathcal{L}}^{-1}\left\{{\displaystyle \frac{1}{s^{1+k\beta }}}\right\}\left(t\right)={\displaystyle \frac{t^{k\beta }}{\mathsf{\Gamma }(1+k\beta )}}.$$

Consequently, the solution operator simplifies to a standard matrix Mittag–Leffler representation:

$$S(t,A)={\displaystyle \sum _{k=0}^{\infty }}{\left[{\mathcal{F}}^k\left(A\right)\right]}_{ij}{\displaystyle \frac{t^{k\beta }}{\mathsf{\Gamma }(1+k\beta )}},$$

or equivalently in operator form, $S(t,A)=E_{\beta }\left(\mathcal{F}\left(A\right)t^{\beta }\right).$ If $\mathcal{F}\left(A\right)$ is diagonalizable, that is $\mathcal{F}\left(A\right)=P\mathsf{\Lambda }\left(A\right)P^{-1}$, then the spectral representation reads $S(t,A)=P\mathrm{diag}\left(E_{\beta }({\lambda }_j\left(A\right)t^{\beta }\right)P^{-1}.$

If $\mathcal{F}\left(A\right)$ is not diagonalizable, we can use its Jordan canonical decomposition $\mathcal{F}\left(A\right)=P(\mathsf{\Lambda }\left(A\right)+N)P^{-1},$ where $\mathsf{\Lambda }\left(A\right)$ is diagonal, N is nilpotent. In this case, the solution operator admits the representation $S(t,A)=P{E}_{\beta }\left((\mathsf{\Lambda }\left(A\right)+N)t^{\beta }\right)P^{-1}.$ Since N is nilpotent, this expansion is finite in the nilpotent part and reduces to a finite polynomial perturbation of the diagonal Mittag-Leffler dynamics.

6. Discussion

The results obtained in this work highlight the versatility and unifying role of solution operator constructions within fractional-order evolution systems. By extending the classical semigroup framework to accommodate distributed-order, memory-decoupled, and fully coupled structures, we demonstrate that a wide range of fractional dynamics can be embedded into a coherent operator-theoretic setting. The fractional Duhamel principle plays a central role in this unification, providing a systematic mechanism for incorporating both initial data and forcing terms through convolution-based memory kernels. This perspective not only clarifies the structural differences between various classes of fractional systems but also reveals their underlying commonalities in terms of history-dependent evolution. Moreover, the obtained operator representations offer a flexible foundation for further analytical investigations, including stability analysis, regularity theory, and potential numerical approximation schemes.