# The Havriliak-Negami and Jurlewicz-Weron-Stanislavsky relaxation models
revisited: memory functions based study

K. Górska  A. Horzela<EMAIL_ADDRESS><EMAIL_ADDRESS>H. Niewodniczański Institute of Nuclear Physics, Polish Academy of Sciences,
ul.Eljasza-Radzikowskiego 152, PL 31342 Kraków, Poland  K. A. Penson
<EMAIL_ADDRESS>b Sorbonne Université, Campus Pierre et
Marie Curie (Paris VI), CNRS UMR 7600, Laboratoire de Physique Théorique de la
Matière Condensée (LPTMC), Tour 13, 5ième ét., B.C. 121, 4 pl. Jussieu, F
75252 Paris Cedex 05, France

###### Abstract

We provide a review of theoretical results concerning the Havriliak-Negami
(HN) and the Jurlewicz-Weron-Stanislavsky (JWS) dielectric relaxation models.
We derive explicit forms of functions characterizing relaxation phenomena in
the time domain - the relaxation, response and probability distribution
functions. We also explain how to construct and solve relevant evolution
equations within these models. These equations are usually solved by using the
Schwinger parametrization and the integral transforms. Instead, in this work
we replace it by the powerful Efros theorem. That allows one to relate
physically admissible solutions to the memory-dependent evolution equations
with phenomenologically known spectral functions and, from the other side,
with the subordination mechanism emerging from a stochastic analysis of
processes underpinning considered relaxation phenomena. Our approach is based
on a systematic analysis of the memory-dependent evolution equations. It
exploits methods of integral transforms, operational calculus and special
functions theory with the completely monotone and Bernstein functions. Merging
analytic and stochastic methods enables us to give a complete classification
of the standard functions used to describe the large class of the relaxation
phenomena and to explain their properties.

##  I Introduction

Researchers working in the fields of natural and engineering sciences are used
to understand relaxation phenomena as processes which describe effects related
to the delay between the application of an external stress to a system and its
response. Phenomena which exhibit such behavior are observed practically
everywhere in our surroundings: decay of induced dielectric polarization or
magnetization usually does not follow switching off the external
electromagnetic field instantaneously and various luminescence phenomena
persist in the absence of illumination. Similarly, deformed viscoelastic
materials return to their original shapes with shorter or longer delay,
kinetics of chemical reactions in many cases is retarded with respect to
changes of external conditions, and so on. The simplest mathematical model
used to describe relaxation is the exponential, or the Debye F1 , law which
provides us with the time decay rule proportional to the exponential function
$\exp{(-\lambda t})$, $\lambda>0$. The Debye law works well for different
physical systems but suffers more and more difficulties or even fails with the
growing complexity of the system under investigation. The first more
systematic observations of non-exponential relaxations and/or decays were
performed more than one hundred seventy years ago, in the middle of the XIX-th
century, by R. Kohlrausch Kohlrausch1854  who carried out a series of
relaxation experiments. At first he was focused on the mechanical creep but
soon after he became interested in electrical phenomena and successfully
measured the residual discharge current in the Leyden jar. Studying his own
experimental results Kohlrausch found that the data were much better described
by the stretched exponent $\exp{(-\lambda t^{\alpha})}$, $\alpha\in(0,1)$,
which means that the relaxation is slower than predicted by the Debye law.
During the next few decades, deviations from the Debye law were observed also
in relaxation phenomena going far beyond mechanics of viscoelastic media and
effects occurring in simple resistor-capacitor circuits. We mention the
luminescence phenomena BarberanSantos05 ; BarberanSantos08  or the response of
dielectric material to the step input of a direct current voltage as well.
Mathematical models tested to explain experimental results were, as a rule,
based on the inverse power-like decay model. Examples include the Curie-von
Schweidler law of current decay or Nutting’s equation for the stress-strain
relation, both usually treated as purely phenomenological input proposed
without any further justification. This phenomenology dominated period of
relaxation studies ended in the context of physics of dielectrics around the
60s and 70s of the previous century. With the progress of experimental
techniques, new materials were more and more common applied and investigated.
Seminal results were obtained by A. K. Jonscher via analysis of a vast
majority of available dielectric relaxation data. He found regularities hidden
behind them Jonscher77 ; Jonscher83 ; Jonscher96 . Jonscher’s discoveries,
unified in the so-called Universal Relaxation Law (URL), boil down to the
statement that the asymptotics of relaxation functions, both in the frequency
and the time domains, are governed by fractional power-like functions. Also at
that time, more than a hundred years after Kohlrausch’s experiments, the
stretched exponential function got back in the game and found wider interest
among physicists thanks to the work of G. Williams and D. C. Watts who applied
the stretched exponential to fit the dielectric relaxation data in the time
domain WW70 . It should be mentioned that despite its popularity among
experimentalists, the use of the stretched exponent (now called also the
Kohlrausch-Williams-Watts (KWW) function) has remained formal as it has never
found a deeper theoretical explanation. Moreover, its use leads to an
embarrassing situation: the KWW function, if transformed from the time to the
frequency domain, is expressed by rather involved special functions, see Sect.
IV. The latter are defined either by special contour integrals in the complex
plane or through slowly convergent power series. This makes it difficult to
read out correctly the dependence of dielectric permittivities on the
frequency of applied harmonic electric field, even if these relations, nick-
named spectral functions, are not only well-known phenomenologically from the
broadband dielectric spectroscopy experiments, but do appear to be simple
rational functions CJFBottcher78 . For that reason, the KWW model is often
replaced by the Cole-Cole (CC) CJFBottcher78 ; RHilfer02 ; RGarrappa16 , Cole-
Davidson (CD) CJFBottcher78 ; RHilfer02 ; RGarrappa16 , Havriliak-Negami (HN)
CJFBottcher78 ; RHilfer02 ; RGarrappa16 , and Jurlewicz-Weron-Stanislavsky
(JWS) AStanislavsky10 ; KWeron10 ; RGarrappa16  patterns. Their spectral
functions are much simpler than the KWW model but the opposite situation
appears for the time dependence functions.

Contemporary ongoing experiments aimed at collecting results essential for a
better understanding of dielectric relaxation phenomena may be classified in a
twofold way. The data which are measured in a suitable experimental setup
concern either the time dependence of physically meaningful properties
characterizing the system (e.g., the number of previously excited dipoles
which survive some amount of time after sudden switching off the polarizing
field or retreating mechanical deformation of the sample), or the frequency-
dependent response of the system perturbed by a step-like or alternating
(usually harmonic) external electric field which plays the role of time-
dependent external stress. In the first case, the results are traditionally
encoded in the time-dependent relaxation $n(t)$ (or the response
$\phi(t)=-\operatorname{d\\!}n(t)/\operatorname{d\\!}t$) function, while in
the second case, customarily used objects are complex-valued spectral
functions $\widehat{\phi}(\operatorname{i}\\!\omega)$. In dielectric physics,
the spectral functions $\widehat{\phi}(\operatorname{i}\\!\omega)$ mimic
normalized complex dielectric permittivity
$[\widehat{\varepsilon^{\star}}(\omega)-\varepsilon_{\infty}]/[\varepsilon_{0}-\varepsilon_{\infty}]$,
whose real and imaginary parts are responsible for diffractive and absorptive
properties of the medium, and the material constants $\varepsilon_{0}$ and
$\varepsilon_{\infty}$ denote the low (static) and high-frequency values of
the dielectric permittivity. The time and frequency approaches are mutually
connected by the Laplace transform F2
$f(t)\div\widehat{f}(z)=\int_{0}^{\infty}\exp{(-zt)}f(t)\operatorname{d\\!}t={\mathscr{L}}[f(t),z]$
which boils down to the following relations between the spectral, response and
relaxation functions:

$\widehat{\phi}(\operatorname{i}\\!\omega)=\mathscr{L}[\phi(t);\operatorname{i}\\!\omega]=1-\operatorname{i}\\!\omega\,\widehat{n}(\operatorname{i}\\!\omega),$
(I.1)

where we used the fact that

$\phi(t)=-\frac{\operatorname{d\\!}n(t)}{\operatorname{d\\!}t}\quad\text{or,
equivalently,}\quad n(t)=1-\int_{0}^{t}\phi(u)\operatorname{d\\!}u.$ (I.2)

Equation (I.1) is crucial for our further consideration. It gives an
operational rule which constitutes a bridge between the time and frequency
descriptions of the relaxation phenomena and opens possibilities to compare
them. This applies also if the data come from independent measurements
performed in the time and frequency domains on the same, or twin-like
prepared, sample.

We shall now give more details on the phenomenological models specified by CC,
CD, HN, and JWS earlier. In the frequency domain the phenomenological pattern
relevant to the Debye model is equal to

$\widehat{\phi}_{D}(\operatorname{i}\\!\omega)=(1+\operatorname{i}\\!\omega\tau)^{-1}$
(I.3)

and non-Debye relaxations are modelized through:

$\displaystyle\begin{split}\widehat{\phi}_{CC}(\operatorname{i}\\!\omega)&=[1+(\operatorname{i}\\!\omega\tau)^{\alpha}]^{-1},\\\
\widehat{\phi}_{C\\!D}(\operatorname{i}\\!\omega)&=[1+(\operatorname{i}\\!\omega\tau)]^{-\beta},\\\
\widehat{\phi}_{H\\!N}(\operatorname{i}\\!\omega)&=[1+(\operatorname{i}\\!\omega\tau)^{\alpha}]^{-\beta},\\\
\widehat{\phi}_{J\\!W\\!S}(\operatorname{i}\\!\omega)&=1-[1+(\operatorname{i}\\!\omega\tau)^{-\alpha}]^{-\beta}.\end{split}$
(I.4)

Real numbers $\alpha,\beta\in(0,1]$ are called the width and the symmetry
parameters, respectively, and each time they are fitted to the experimental
data together with the material-dependent characteristic time scale $\tau$. It
is easy to see that Eqs. (I.4) extend the Debye rule and that among them the
HN and JWS patterns are the most general involving the CC and CD as their
special cases CJFBottcher78 ; RGarrappa16 . The Debye and CC models emerge
either from $\widehat{\phi}_{H\\!N}(\operatorname{i}\\!\omega)$ and
$\widehat{\phi}_{J\\!W\\!S}(\operatorname{i}\\!\omega)$ for $\alpha=\beta=1$
or for $\beta=1$ and $\alpha\in(0,1)$, respectively. The CD model comes out
from the HN pattern for $\alpha=1$ and $\beta\in(0,1)$. On the other hand,
measurements done in the time domain provide us with the data for the
relaxation function which are usually fitted by the standard KWW function

$n_{KWW}(t)=\exp{[-(t/\tau^{\prime})^{\alpha^{\prime}}]}$ (I.5)

where parameters $\tau^{\prime}$ and $\alpha^{\prime}\in(0,1]$ are unrelated
to those of Eqs. (I.4). Moreover, none of the formulae listed in Eq. (I.4) are
the Laplace transform of Eq. (I.5). Efforts to shed more light on this
troublesome situation have been the subject of both experimental and
theoretical research accompanying the issues of dielectric physics Alvarez91 ;
Alvarez93 ; KGorska21c ; HandH96 ; RHilfer02 ; RHilfer02a .

Focusing attention on the HN pattern we remind that it was introduced in
SHavriliak67  to explain the asymmetry of the Cole-Cole plot, i.e., the
Argand-like diagram in the
$\widehat{\varepsilon^{\prime}}(\omega)=\operatorname{\mathfrak{Re}}[{\widehat{\varepsilon^{\star}}(\omega)}]$,
$\widehat{\varepsilon^{\prime\prime}}(\omega)=\operatorname{\mathfrak{Im}}[{\widehat{\varepsilon^{\star}}(\omega)}]$
plane, observed in dielectric relaxation of polymers. Shortly after Havriliak
and Negami’s discovery, the HN pattern gained popularity because, although
involving only three parameters, it appeared to be so flexible and universal
that within a few years it became a standard model used to describe data
obtained for a large range of relaxation and viscoelastic phenomena.
Consequently, the calculation of its time domain counterpart, leading either
to the response or to the relaxation function, became a serious challenger to
routine parameterizations of their time evolution, based on the KWW functions
or finite sums of exponential decays. This has motivated many authors to
investigate the time dependence of the HN model and explore it using different
mathematical tools. They include the analytical methods like the integral
transforms theory and fractional calculus as well as an the approach worked
out in the probability and stochastic processes theory and known as
subordination methods. We shall show that these seemingly distant approaches
merge entirely into the framework of memory-dependent kinetics.

The general expressions for the time-dependent counterparts of Eqs. (I.4) (as
well as the frequency-dependent analogue of Eq. (I.5)) were found analytically
20 years ago by R. Hilfer RHilfer02 ; RHilfer02a  who, for arbitrary real
values of $\alpha$ and $\beta$, expressed functions relevant for dielectric
relaxation in terms of the Fox H-functions. Note that the applicability of
this class of functions was used 10 years ealier for studies of other
relaxation phenomena, especially in the research devoted to rheological models
of viscoelastic materials Gloeckle91 ; Gloeckle93 ; RMetzler95 ; RMetzler03 ;
TNonnenmacher91 ; HSchiessel95 ; NWTschoegl . In applications essential for
the relaxation phenomena, the Fox H-functions are replaced by much better-
known functions of the Mittag-Leffler family RGarrappa16  and the Meijer
G-function. In fact, its well-known representative, namely the standard
Mittag-Leffler function, naturally enters the time relaxation for the CC
model. Following and exploring this direction of theoretical investigations we
shall show how to express the functions relevant to the HN and JWS relaxation
models, like the response or the relaxation functions and the probability
densities, in terms of the Meijer G-function and the generalized Mittag-
Leffler functions. These families of functions are well-implemented in modern
computer algebra systems (CAS) and may be effectively studied using methods of
classical analysis much more accessible to a larger community. The crucial
step to achieve this purpose is to replace the irrational values of the
parameters sitting in the Fox H-functions with rational numbers. This simple,
but physically well-justified F3  approximation, enables us, as the first
step, to represent results in terms of the Meijer G-functions, and in the
second step in terms of the generalized hypergeometric functions which
specialize to the generalized, or multi-parameter, Mittag-Leffler functions.

##  II The time domain description: general properties of the relaxation
functions

Assume that the object of our investigation is the time-dependent quantity
$\mathcal{J}(t)$, like the relaxation or the response function is, which
reflects the time behavior of some, suitably chosen, property $\mathcal{P}(t)$
of the system. For example, $\mathcal{P}(t)$ can be the number of its
constituents which at an instant of time $t$ are still excited (have survived
the decay), or the depolarization current, or the intensity of luminescence
light, etc. The normalized form of ${\mathcal{J}}(t/\tau)$, with the
characteristic time scale $\tau$ put in, reads:

${\mathcal{J}}(t/\tau)=\frac{\mathcal{P}(t/\tau)-\mathcal{P}_{\infty}}{\mathcal{P}_{0}-\mathcal{P}_{\infty}}$
(II.1)

where $\mathcal{P}_{0}$ and $\mathcal{P}_{\infty}$ denote the values of
$\mathcal{P}(t)$ for $t=0$ and $t=\infty$, respectively. Parameterization
given by Eq. (II.1) implies by construction that ${\mathcal{J}}(0)=1$ and
${\mathcal{J}}(\infty)=0$. We also assume that $\mathcal{J}(t)$ is
experimentally measurable (at least in the sense of Gedankenexperiment) and
that in the system, there exist no internal factors which can revert the
decay, i.e., $\mathcal{J}(t)$ never increases. Such requirements are too vague
to specify $\mathcal{J}(t/\tau)$ satisfactorily and additional assumptions are
needed to make any prospective mathematical modeling of $\mathcal{J}(t/\tau)$.
Two commonly used suggestions are:

######  Proposition 1.

${\mathcal{J}}(t/\tau)=\exp\left(\\!-\int_{0}^{t}r(t^{\prime}/\tau)\operatorname{d\\!}t^{\prime}\\!\right),$
(II.2)

with

$r(t/\tau)=-\tau\frac{\operatorname{d\\!}\,\ln{\mathcal{J}}(\xi)}{\operatorname{d\\!}\xi}|_{\xi=t/\tau},$
(II.3)

which generalizes the standard exponential Debye decay law,
$\operatorname{d\\!}n(t)/\operatorname{d\\!}t=-\lambda n(t)$, by extending it
to the case when the relaxation/decay rate $\lambda$ becomes the time-
dependent quantity, $\lambda\rightarrow r(t)$ instead of being a constant.

######  Proposition 2.

$\displaystyle\begin{split}{\mathcal{J}}(t/\tau)&=\int_{0}^{\infty}p(x)\exp{(-xt/\tau)}\operatorname{d\\!}x\\\
&=\tau\int_{0}^{\infty}p(\tau\xi)\exp{(-\xi
t)}\operatorname{d\\!}\xi\end{split}$ (II.4)

is interpreted as a continuous weighted sum of the Debye decays, provided the
additional condition that $p(x)$ is the probability distribution function
(PDF), i.e., it is non-negative and normalizable. If we focus attention on the
relaxation phenomena and identify ${\mathcal{J}}(t/\tau)$ either with the
relaxation $n(t)$, or the response function $\phi(t)$, then both these
functions may be represented as the Laplace integrals F4

$n(t)=\int_{0}^{\infty}\operatorname{e}^{-\frac{t}{\tau}\xi}g(\xi)\operatorname{d\\!}\xi$
(II.5)

and

$\phi(t)=-\frac{\operatorname{d\\!}n(t)}{\operatorname{d\\!}t}=\int_{0}^{\infty}\operatorname{e}^{-\frac{t}{\tau}\xi}\,\frac{\xi}{\tau}\,g(\xi)\operatorname{d\\!}\xi.$
(II.6)

If we have enough information on the time dependence of $n(t)$ or $\phi(t)$,
e.g., know that their time decay is governed by some KWW-like function, then
Eqs. (II.5) and (II.6) may reduce to the well-defined problems of classical
mathematical analysis whose solutions allow us to determine the function
$g(\xi)$, see e.g. HPollard46 .

The above propositions are strongly supported by simple physical intuition
and, at first glance, look very attractive as a way to solve the problem.
Indeed, imagine that the constituents of the relaxing system are under
consideration, e.g., polarized dipoles, decay according to the Debye law but
each constituent of the system does that with a different decay rate
$\lambda_{k}$; we treat the latter quantity as a non-negative discrete or
continuous random variable whose distribution $g(\lambda)$ we somewhat know.
We also assume that the randomization of elementary decays does not influence
the Debye pattern dominating the process. Under such formulated conditions of
independence the use of a joint probability formula for the fraction of
constituents which survive the decay gives

$\frac{n(t)}{n_{0}}=\int_{0}^{\infty}e^{-\lambda
t}g(\lambda)\operatorname{d\\!}\,\lambda,$ (II.7)

i.e., we have arrived at an expression which has an obtrusive interpretation
to be treated as a weighted sum of exponential decays. Analogous reasoning is
the basis of the majority of routine analyses whose goal is to represent decay
curves describing experimental data as a weighted finite (or infinite) sum of
exponentials RSAnderssen11 ; Bertelsen74 ; Feldman02 ; Johnston06 ;
Nigmatullin12 . However, despite their attractive simplicity and possible
utility, the statements cited above are too naive. Equations (II.3) and (II.4)
suffer from essential physical shortcomings and lack of mathematical
precision. Equation (II.3), if treated as a source of experimental information
concerning $r(t)$, is difficult to be verified both for short and long times.
Another objection is that $r(t)$ calculated from exactly solvable fractional
models of relaxation becomes singular for $t\to 0$ and thus unphysical, e.g.,
KGorska21CNSNCa . The main criticism against Eq. (II.4) comes from the
objection that the interpretation of factors $\exp{(-\xi t)}$ as the Debye
decays, each characterized by its own decay rate $\xi$, is formal and
difficult to be explained. Does it mean that decay rates appearing in such a
way are distributed according to some mysterious function $p(\tau\xi)$? What
is the origin of the time scale $\tau$ introduced ad hoc because of
dimensional reasons but finally acquiring universal meaning? Nevertheless, so
far proposed Eqs. (II.3) and (II.4), if understood in mathematically and
physically correct ways, provide us with useful guideposts which show how to
choose, arrange and push forward theoretical methods whose development will
lead to a better understanding of relaxation phenomena.

The mathematical condition which guarantees that Eq. (II.4) makes sense, comes
from classical mathematical analysis - according to the Bernstein theorem (for
a brief tutorial on it, on completely monotone (CM) and on Bernstein (B)
functions see Sect. IV.1) the sufficient and necessary condition to represent
the function $f(t)$ through the Laplace integral of the measure $g(x)$

$f(t)=\int_{0}^{\infty}\exp{(-xt)}g(x)\operatorname{d\\!}x$ (II.8)

is that $f(t)$ is CM. In Eq. (II.8) we deal with $t\in\mathbb{R}_{+}$. It
means that we work in the realm of real-valued Laplace integral whose
continuation to the complex-valued Laplace transform needs care. Merging Eqs.
(II.3) and (II.4) and requiring complete monotonicity of
${\mathcal{J}}(t/\tau)$ represented by Eq. (II.2), one learns that the minus
exponent, which is in its right-hand side (RHS), i.e.,
$\int_{0}^{t}r(t/\tau)\operatorname{d\\!}t$, must belong to the class of
Bernstein functions, $\int_{0}^{t}r(t/\tau)\operatorname{d\\!}t\in{\rm{B}}$.
Thus, $r(t)$ must be the CM and, consequently, must not be singular. As
mentioned earlier, solvable models, e.g., the CC, show that it is not the
case, which excludes them from further considerations. The distinguished role
of the CM functions in the theory of relaxation phenomena has been noted a
long time ago, starting with early observations concerning viscoelasticity
Day70  and rheology Bagley83 . Additional links between the CM properties and
relaxation phenomena have been provided by the studies of Green’s functions
Anh03  and of the extended family of the Mittag-Leffler-type functions
ECapelas11 ; RGarrappa16 ; KGorska21b ; Hanyga1 ; FMainardi15 ; Tomovski14 ,
describing of the non-Debye relaxation.

Besides their well-established place in the classical mathematical analysis
Berg , the CM and B functions are objects which play an important role also in
the probability and the theory of stochastic processes RLSchilling10 . Thus,
their appearance in the formalism signalizes existence of non-accidental
relations between probabilistic tools and analytical methods of fractional and
operational calculi, if applied to the relaxation theory. Requiring
probabilistic interpretation of Eq. (II.4), besides of the indication that
${\mathcal{J}}(t)$ should be CM, leads to another important result. If we
calculate the Laplace transform of ${\mathcal{J}}(t)$ with respect to $t$,
$\widehat{{\mathcal{J}}}(z)={\mathscr{L}}[{\mathcal{J}}(t),z]$,
$z\in\mathbb{C}$, then, after changing the order of integration, we get

$\widehat{{\cal J}}(z)=\int_{0}^{\infty}\frac{p(t)}{t+z}\operatorname{d\\!}t.$
(II.9)

Under the condition that $p(t)$ is integrable and non-negative for all
$t\in{\mathbb{R}_{+}}$, it means that $\widehat{{\mathcal{J}}}(z)$ belongs to
the class of Stieltjes functions (S) which obeys well-known analytical
properties, play important role in the Laplace transform theory Widder  and,
as we will see, turns out to be of a crucial importance in our further
analysis. Notice that to get Eq. (II.9), we assume
$\operatorname{\mathfrak{Re}}{(z)}\in{\mathbb{R}}_{+}$. From the physical
point of view, this assumption is too restrictive - to make our scheme
compatible with physical data incorporated in Eq. (I.4), we do need to include
also $\operatorname{\mathfrak{Re}}{(z)}=0$. The justification of this
statement comes from the fact that the Laplace transform Eq. (II.9) for
$z=\operatorname{i}\\!\omega$ may be considered as the Fourier transform of a
function vanishing on the negative semiaxis, i.e., the Fourier transform of a
function $\Theta(t)f(t)$ being undefined for $t=0$. The fundamental physical
implication of the Fourier transform, namely the time-frequency
correspondence, implies that
$\widehat{{\mathcal{J}}}(z)|_{z=\operatorname{i}\\!\omega}$ in such a way
understood is closely related to the earlier introduced spectral functions
being objects of direct physical meaning, as reconstructed from experimental
data provided by methods of the broadband dielectric spectroscopy. But the
spectral functions should be (and they are) defined also for $\omega=0$, in
the case of a static field. Mathematically, it means that we have to point out
the value $\widehat{{\mathcal{J}}}(0)$ to redefine the singularity of
$\int_{0}^{\infty}\frac{p(t)}{t}\operatorname{d\\!}t$. Here the Plemelj-
Sokhotski formula F5  enters the game

$\int_{0}^{\infty}\frac{p(t)}{t}\operatorname{d\\!}t=\operatorname{i}\pi
p(0)+{\rm{PV}}\int_{0}^{\infty}\frac{p(t)}{t}\operatorname{d\\!}t$

with the Cauchy principal value

${\rm{PV}}\int_{0}^{\infty}\frac{p(t)}{t}\operatorname{d\\!}t=\lim_{\epsilon\to
0}\int_{\epsilon}^{\infty}\frac{p(t)}{t}\operatorname{d\\!}t$

being introduced. Consequently, we should extend Eq. (II.9) to the form

$\widehat{{\mathcal{J}}}(z)={\rm{const}}+{\rm{PV}}\int_{0}^{\infty}\frac{p(t)}{t+z}\operatorname{d\\!}t,$

which explains different shapes of the spectral functions relevant for the HN
and JWS models.

Summarizing the results obtained so far - using intuitive, and from time to
time also hand-waving arguments, we can expect that functions admittable to
describe non-exponential relaxation phenomena should be searched among
functions belonging to the CM class or eventually its subclasses. Natural
justification of such a conjecture is that the CM functions obey
characteristic properties of relaxation functions, namely they are always non-
negative and non-increasing. Infinite differentiability and sign-alternating
derivatives, which in fact define complete monotonicity, are impossible to be
reliably verified for functions known from phenomenological analyses only.
Intuitively, to require the CM property shared with the exponential function
used to describe the standard Debye decay, can shed more light on the problem.
In Secs. V and VI we will see that it is really the case. Now, leaving aside
the CM functions, we are going to convince the readers that the most promising
factor which will enable us to go forward is to consider memory-dependent
kinetics giving dynamical laws governing the relaxation phenomena.

##  III Memory dependent kinetics

As signalized in the Introduction, the characteristic feature of the
relaxation phenomena is that the relaxing system exhibits delayed response to
the changes of external conditions. A typical example of such a behavior is
extended over time approach to equilibrium after switching off factors
perturbing the system during its previous history. In other words -
approaching equilibrium, even in the case of absent external influence, needs
time and the process under consideration does not satisfy the instantaneous
response principle. This means giving up the time-local evolution rules and
strongly suggests looking for prospective challengers which would be used to
replace the standard time-local evolution equations. The “first choice”
solution of the problem is to take into account the time non-local analogues
of basic equations, in particular to consider integro-differential equations
with kernels suited to mimic memory effects occurring in the system. Thus, to
study the non-Debye relaxation phenomena, as well as other kinetic phenomena
taking place in complex systems, e.g., anomalous diffusion, we should to
reformulate the standard time-local approach. The aim is to construct a new
theoretical scheme which has the time non-locality built-in from the very
beginning, and leads to new evolution laws modified by various memory effects.
The cornerstone of such new formalism is classical works of R. Zwanzig
Zwanzig1 ; Zwanzig2  who, after paying attention to the fact that the
formalism of instantaneous response did not work in many physically
interesting situations, proposed the self-consistent construction of kinetic
equations respecting causality and taking into account delayed response of the
system. The essence of mathematical encoding such delays is to replace point-
like operations, e.g., the usual multiplication of functions, by integral
operators taken in the form of time convolutions of memory functions and
solutions looked for. Thus, it is justified to say that such modified
evolution equations are the "time-smeared" standard ones. Almost 40 years
later, at the very beginning of the current century, it was I. M. Sokolov
Sokolov1  who shed new light on the ideas underpinning Zwanzig’s seminal
works. Sokolov showed that the use of generalized non-Markovian Fokker-Planck
(FP) equations with built-in memory kernels leads to considerable progress in
understanding anomalous transport and non-Debye relaxations. Subsequent years
brought new research techniques and huge amount of earlier unreachable
experimental data. These, urgently needed to be theoretically explained,
pushed forward theoreticians’ efforts and soon led to the development of
effective mathematical tools rooted either in mathematical analysis, e.g.,
provided by the fractional calculus, or coming from the probability/stochastic
processes theory. Taken separately or merged together, these methods have
formed a toolbox whose usage has enabled physicists to solve equations of
anomalous diffusion and non-Debye relaxations for a quite large number of
problems AChechkin21 ; KGorska21 ; TSandev18 ; TSandev17 ; AStanislavsky15 .

###  III.1 The basics

Restricting Zwanzig’s approach to the problems depending on the time only,
e.g., to the simplest analysis of dielectric relaxation F6 , but holding on to
the spirit of seminal paper Zwanzig1  we can write down the evolution equation
as

$\frac{\operatorname{d\\!}n(t)}{\operatorname{d\\!}t}={\mathcal{O}}n(t),$
(III.1)

where $t\in\mathbb{\mathbb{R}_{+}}$. The linear operator ${\mathcal{O}}$ acts
on the time variable but is not simple multiplication by a time-dependent
function, and is the object which represents a more general characterization
of non-locality in time. If ${\mathcal{O}}$ is represented as an integro-
differential operator with the kernel $O(t-\xi)\mathcal{B}$, $\mathcal{B}={\rm
const}$, then Eq. (III.1) becomes

$\frac{\operatorname{d\\!}n(t)}{\operatorname{d\\!}t}=\int_{0}^{t}\\!\\!O(t-\xi)\mathcal{B}n(\xi)\operatorname{d\\!}\xi,$
(III.2)

i.e., has a typical form of the Volterra equation. It is convenient to rewrite
Eq. (III.2) as

$\frac{\operatorname{d\\!}n(t)}{\operatorname{d\\!}t}=-\frac{\operatorname{d\\!}}{\operatorname{d\\!}t}\\!\int_{0}^{t}\\!\\!M(t-\xi)\mathcal{B}n(\xi)\operatorname{d\\!}\xi$
(III.3)

coming directly from the integral master equation AStanislavsky17 ;
AStanislavsky19b

$n(t)=n(0)-\int_{0}^{t}M(t-\xi)\mathcal{B}n(\xi)\operatorname{d\\!}\xi.$
(III.4)

If we take the Laplace transforms (with respect to $t$) of Eqs. (III.2) and
(III.3) or (III.4) and symbolically write down these equations as

$[z-\widehat{O}(z)\mathcal{B}]\,\widehat{n}(z)=z[1-\widehat{M}(z)\mathcal{B}]\,\widehat{n}(z)=n(0)$

then we see that Eqs. (III.2) - (III.4) are equivalent, if in the Laplace
domain $\widehat{O}(z)=z\widehat{M}(z)$, i.e., in the time domain
$O(t)=\operatorname{d\\!}M(t)/\operatorname{d\\!}t+M(0+)$, where
$M(0+)=\lim_{t\rightarrow 0_{+}}M(t)$ is to be specified, cf. AChechkin21 .
Looking for solutions of (III.2) - (III.4) we are interested in the functions
$n(t)$ which are normalizable and non-negative, i.e., fulfill a minimal set of
requirements needed to interpret them as the relaxation functions. The crucial
property to be shown is non-negativity - we will see that requiring this
property allows to find conditions which have to be put on the memory
functions.

The theory of integral equations, more precisely the possibility of a
dichotomic description of the system using either integral or differential
equations, suggests that it would be interesting also for Eq. (III.4) to find
its integro-differential analogue. To achieve this goal let us notice that if
we have given Eq. (III.3) governed by the kernel function $M(t)$ then we may
ask for a function $k(t)$ defined by the convolution

$\int_{0}^{t}k(t-t^{\prime})M(t^{\prime})\operatorname{d\\!}t^{\prime}=1.$
(III.5)

which, if required to be satisfied for any $t>0$, is known in the theory of
integral equations as the Sonine relation AHanyga20 ; AHanyga21 ; YLuchko21 ;
YLuchko21a ; Meerschaert19a . Transformed to the Laplace domain Eq. (III.5)
becomes $\widehat{k}(z)\widehat{M}(z)=z^{-1}$ from which it follows:

$\widehat{k}(z)=[z\widehat{M}(z)]^{-1}.$ (III.6)

General properties required from $M(t)$ in order to satisfy basic conditions
of the Laplace transform theory guarantee that Eq. (III.5) defines the couple
$(M,k)$ uniquely. Moreover, the symmetry of Eq. (III.6) with respect to the
interchange $M(t)\leftrightarrow k(t)$ means that if $M(t)$ is interpreted as
a memory function then $k(t)$ does share this property: if Eq. (III.5) holds
then any suitable $M(t)$ which governs Eq. (III.3) has its partner memory
function $k(t)$ entering the kernel of integro-differential equation

$\int_{0}^{t}k(t-\xi)\frac{\operatorname{d\\!}{n}(\xi)}{\operatorname{d\\!}\xi}\operatorname{d\\!}\xi=-\mathcal{B}n(t),\qquad\mathcal{B}={\rm
const}.$ (III.7)

Equation (III.7) is equivalent to Eq. (III.3) and in mathematics it is known
as the Caputo-Djhrbashyan problem AKochubei11 . The equivalence of Eqs.
(III.3) and (III.7) is seen if we calculate the time derivative in the RHS of
Eq. (III.3) and next integrate the resulting equation with respect to
$\int^{t}_{0}\operatorname{d\\!}t^{\prime}$ using the Sonine relation (III.5)
twice KGorska21a ; KGorska20a . Duality in the description of the same process
using either Eq. (III.3) or Eq. (III.7) was noticed and studied in many papers
KGorska20a ; KGorska21a ; TSandev17 ; TSandev18 . Those were first of all
addressed to the anomalous diffusion (through studying various applications of
the generalized Fokker-Planck equation TSandev17 ; TSandev18 ) but led to
results by no means limited to this class of phenomena KGorska20a ; KGorska21a
. To name different points of view on the problem we list: i.) solvability of
the Cauchy problem for Eq. (III.7) studied either using advanced partial
differential equations methods AKochubei11  or operator calculus tools
developed for the Volterra-type evolution equations JPruess93  and leading to
the appearance of subordination Bazhlekova18 ; Bazhlekova19 , ii.)
investigations devoted to the role played by the Sonine relation and resulting
duality of memories KGorska20a ; KGorska21CNSNCa ; KGorska21a ; KGorska21c ;
AHanyga20 ; AStanislavsky21  and iii.) efforts oriented on understanding
mutual relations between the so-called deterministic and
probability/stochastic approaches to anomalous kinetic phenomena KGorska21 ;
TSandev18 ; TSandev17 ; AStanislavsky17 ; AStanislavsky19b .

Recall that restricting the kinetic problem under consideration to be
depending on the time only means that the action of the FP operator reduces to
multiplication by a constant factor $\mathcal{B}$. Thus we arrive at equations
governing the relaxation phenomena which provide us with other examples of
non-Markovian processes AStanislavsky19b ; AStanislavsky15 . It has to be
noted that Eq. (III.4) enables us to write down the spectral function, i.e.,
the Laplace transform $\widehat{\phi}(z)$ of the response function
$\phi(t)=-\operatorname{d\\!}{n}(t)/\operatorname{d\\!}t$ in terms of
$\widehat{M}(z)$:

$\widehat{\phi}(z)=1-z\widehat{n}(z)=\\{1+[{\mathcal{B}}\widehat{M}(z)]^{-1}\\}^{-1},$
(III.8)

reasulting in

$\widehat{M}(z)={\mathcal{B}}^{-1}\frac{\widehat{\phi}(z)}{1-\widehat{\phi}(z)},$
(III.9)

which gives the mutual relation between the spectral functions and evolution
kernels responsible for the memory effects. This observation has far-going
consequences. From the physicists point of view Eqs. (III.8) and (III.9), if
taken for $z=\operatorname{i}\omega$ with $\omega$ identified as the frequency
of harmonic field used to polarize the sample, couple in a unique way
evolution equations introduced on theoretical background with
phenomenologically known spectral functions given by formulae which classify
the standard non-Debye relaxations as the CC, CD, HN, and JWS patterns (cf.
KGorska21CNSNCa ; KGorska21a ; KGorska21c ; KGorska20a ).

###  III.2 The stochastic approach - a few general remarks

As mentioned, Eqs. (III.3), (III.4), and (III.7) are the Volterra-type
equations GGripenberg90 . They may be considered as equations which, because
of a possible introduction of "atypical" kernels $M(t)$ and $k(t)$, go beyond
the fractional differential equations conventionally used to describe the
anomalous diffusion and non-Debye relaxation phenomena. In turn, Eq. (III.4)
may be identified as the master equation which governs some stochastic
processes underlying mechanisms of anomalous kinetics AStanislavsky19b .
Recall that the memory functions $M(t)$ and $k(t)$ play a dual role in our
scheme - they determine the mathematical structure of equations under study
and provide us with links to observational data fitted in the case of
relaxation experiments by phenomenological spectral functions
$\widehat{\phi}(\operatorname{i}\\!\omega)$. However, if we restrict ourselves
to the experimental data only, then usually we are not able to acquire
knowledge sufficient to determine basic processes underlying physics of
considered phenomena. To proceed further we do need more detailed information,
also coming from mathematics, expected to be obtained from analysis of
properties which the memory functions should obey. Obvious requirements that
the memory functions have to be non-negative and non-increasing are
insufficient to judge the existence and physical applicability of solutions as
well as to find their interpretation. Some extra conditions, like the
Boltzmann fading memory concept, did not satisfy initial expectations due to
the lack of mathematical precision RSAnderssen02 ; RSAnderssen02a .
Clarification and resolution of many doubts appearing in the statistical
description of both anomalous transport and relaxation phenomena came with
using advanced probabilistic methodology, first of all with identification of
just mentioned phenomena with stable stochastic processes.

Stochastic point of view (cf. e.g., Meerschaert01 ; Meerschaert19 ;
AStanislavsky17 ; AStanislavsky15  and (AStanislavsky19b, , Chs. 4.1, 4.3))
led to methods which made it possible to analyse non-Debye relaxations (as
well as other anomalous kinetics processes) using probabilistic concepts like
subordination Bochner1 ; Bochner ; Feller ; RLSchilling98 , infinitely
divisible distributions Lukac  and applications the theory of B and S
functions, see Berg ; RLSchilling10  and Sect. IV.1 for a brief tutorial.
Simultaneously, investigations based on less popular methods of mathematical
analysis (fractional and operational calculi) made the integro-differential
equations (III.3) and (III.7) much better understood if they were governed by
kernels whose the Laplace images belong to the S class KGorska21 ; AKochubei11
; AKochubei19 . Here we would like to turn once more the readers attention to
the importance of the Sonine relation and the special role of the S functions
hidden behind it. Indeed, if one requires that $\widehat{M}(z)$ and
$\widehat{k}(z)$ of Eq. (III.6) belong to the same class of functions then the
most natural and convincing possibility to satisfy this requirement is to put
both of them in the S class. Additionally, belonging to the S class directly
links the memory functions with the complete Bernstein (CB) functions and the
Laplace exponents, objects which characterize PDFs relevant for infinitely
divisible stochastic processes KGorska21CNSNCa ; KGorska21a ; KGorska21c ;
RLSchilling98 ; AStanislavsky19b .

###  III.3 Integral decompositions and subordination

To keep the forthcoming construction as general as possible let us observe
that for the relaxation phenomena $\widehat{n}(z)$ from Eq. (III.8) may be
rewritten in the form

$\widehat{n}(z)=\frac{[z\widehat{M}(z)]^{-1}}{\mathcal{B}+[\widehat{M}(z)]^{-1}}=\int_{0}^{\infty}\operatorname{e}^{-\xi\mathcal{B}}\,\frac{\widehat{\Psi}(z)}{z}\operatorname{e}^{-\xi\widehat{\Psi}(z)}\operatorname{d\\!}\xi.$
(III.10)

To get Eq. (III.10) the Schwinger parametrization was used and
$\widehat{\Psi}(z)=[\widehat{M}(z)]^{-1}$ introduced. $\widehat{\Psi}(z)$ is
called the characteristic or the Laplace-Lévy exponent AStanislavsky17 ;
AStanislavsky19b  which can be connected to the spectral function
$\widehat{\phi}(z)$ via Eq. (III.9) :

$\widehat{\Psi}(z)=\mathcal{B}\,\frac{1-\widehat{\phi}(z)}{\widehat{\phi}(z)}.$
(III.11)

Taking the inverse Laplace transform of Eq. (III.10) we express it in the form
of integral decomposition:

$n(t)=\int_{0}^{\infty}n_{D}(\xi)f_{\widehat{\Psi}}(\xi,t)\operatorname{d\\!}\xi$
(III.12)

the integrand being a product of the Debye relaxation function
$n_{D}(\xi)=\exp(-\mathcal{B}\xi)$ and
$f_{\widehat{\Psi}}(\xi,t)=\mathscr{L}^{-1}[\widehat{\Psi}(z)\exp[-\xi\widehat{\Psi}(z)]/z;t]$
presenting the relation between the time $t$ and integral variable $\xi$. The
functions $n_{D}(\xi)$ and $f_{\widehat{\Psi}}(\xi,t)$ are given by normalized
and non-negative functions which in probability theory means that they are the
PDFs of “parent” F7  and “leading” processes, respectively. If these PDFs are
independent Eq. (III.12) turns out to be the subordination. Within the latter
scheme the integral variable $\xi$ is interpreted as the internal time with
randomize the physical time $t$.

The same results can be obtained by using the Efros theorem (Subsect. IV.2)
giving more flexibility than an application of the Schwinger parametrization.
The use of the Efros theorem allows one to preset $n(t)$ as

$n(t)=\int_{0}^{\infty}h(u)f(u,t)\operatorname{d\\!}u,$ (III.13)

where $h(u)$ denotes the relaxation function of basic model. It is
subordinated by the PDF $f(u,t)$ which expresses the conversion between
physical $t$ and internal $u$ times. The simplest realization of Eq. (III.13)
is for $f(u,t)=\delta(t-u)$ and $h(u)=n(u)$ where $\delta(\cdot)$ denotes the
$\delta$-Dirac distribution. Other examples of Eq. (III.13) are presented in
Subsecs. V.6 and VI.6 where for the HN and JWS models setting separately
$h(u)$ and $f(u,t)$ we obtain dual representations of the same $n(t)$.

##  IV Interlude - elements of the mathematical toolbox

The wide applicability of the relaxation pattern attracts numerous authors to
investigate it from various points of view. The goal of mathematically
oriented research was to find the analytic expressions for the family of
relaxations in the time domain. We mentioned in the Introduction that the
problem was solved by R. Hilfer RHilfer02 ; RHilfer02a  who showed that
solutions looked for may be expressed in terms of the Fox H-functions. This
result, important for theoretical considerations, regrettably did not find
deeper practical interest among experimentalists. It is not difficult to
understand that as the Fox H-functions are special functions which for general
values of parameters are defined symbolically either via the contour integrals
of the Mellin-Barnes type APPrudnikov-v3  or by slowly convergent series. All
together it makes relevant calculations difficult to understand, time-
consuming and prone to mistakes emerging during hand-made manipulations. Thus,
it is justified to emphasize that the use of the Fox H-functions formalism is
not mandatory to invert analytically the Laplace transform in Eq. (I.4) if its
RHS specializes in specifically chosen sets of parameters. As is shown in a
series of research papers and books, e.g., Prabhakar ; book1  in such a case
the suitable Laplace transform essentially simplifies and results in functions
of the Mittag-Leffler family, much better known and easier to be handled both
analytically and numerically. Moreover, relevant Mittag-Leffler functions obey
properties placing them in special classes of analytic functions which for the
positive value of the argument, besides of being non-negative, also boils down
to the completely monotone functions Berg ; RLSchilling10 .

###  IV.1 Non-negatively definite functions

#### The completely monotone (CM) functions

are introduced as non-negative and infinitely differentiable (i.e., belonging
to the $C^{\infty}$ class) functions $\widehat{g}(s)$ on $\mathbb{R}_{+}$
whose all derivatives alternate in sign for any $s\in\mathbb{R}_{+}$:

$(-1)^{n}{\widehat{g}}^{(n)}(s)\geq 0,\quad n=0,1,\ldots.$

Notice that any CM function is a non-increasing and convex function.
According to the Bernstein theorem RLSchilling10 , we can connect in a unique
way CM and non-negative functions:
$s\in[0,\infty)\rightarrow\widehat{g}(s)\in\textrm{CM}$ iff

$\widehat{g}(s)=\int_{0}^{\infty}\exp(-su)g(u)\\!\operatorname{d\\!}u$

with $g(u)\geq 0$ for all $u\in[0,\infty)$. The important property of the CMs
is that the product of two CM functions is CM as well, i.e., $CM\cdot
CM\subset CM$.

#### The Stieltjes (S) functions

are denoted as $\widehat{k}(s)$ and form a subclass of the CM functions. On
the real semiaxis they admit the integral representation

$\widehat{k}(s)=\frac{a}{s}+b+\int_{0}^{\infty}\frac{\sigma(\operatorname{d\\!}u)}{s+u},\qquad
s\in\mathbb{R}_{+},$ (IV.1)

where $a,b\geq 0$ and $\sigma$ is a Borel measure on $[0,\infty)$, such that
$\int_{0}^{\infty}(1+u)^{-1}\sigma(\operatorname{d\\!}u)<\infty$. Equation
(IV.1) can be interpreted as the Laplace transform of the CM function, i.e.,
the Laplace transform of the Laplace transform of the non-negative function
Berg ; RLSchilling10 . This fact was pointed out and commented in the
Introduction.

#### The Bernstein (B) functions

are these non-negative infinitely differentiable functions $\widehat{b}(s)$ on
$\mathbb{R}_{+}$ whose first derivative is CM. It means that

$(-1)^{n-1}\,\widehat{b}^{(n)}(s)\geq 0,\quad n=1,2,\ldots.$

All B functions are non-decreasing and concave. As the main property, we
remark that the composition of CM function with B function gives another CM
function, i.e., $CM(B)\subset CM$. This property is illustrated on the example
of $\exp[-b(s)]$ which for $b(s)\in B$ is CM function.

#### The complete Bernstein (CB) functions,

denoted as $\widehat{c}(s)$, $s\in\mathbb{R}_{+}$, are CB functions such that
$\widehat{c}(s)/s$ is the Laplace transform of the CM function restricted to
the positive semi-axis, or, equivalently, the same way restricted Stieltjes
transform of a positive function.
In the paper we will use the following properties of CB functions:

  * (a1)

the composition of two CB functions is CB function, i.e., $CB(CB)\subset CB$;

  * (a2)

S functions can be obtained by making the composition of another S function
with a CB function or the composition of a CB function with another S
function, i.e., $S(CB)\subset S$ and $CB(S)\subset S$;

  * (a3)

reciprocal (algebraic inverse) of an S function is a CB function. The inverse
property is true as well: $1/\widehat{k}(s)$ is CB function and
$1/\widehat{c}(s)$ is a S function for $\widehat{k}\in S$ and $\widehat{c}\in
CB$.

The illustrative example of differences between CM, S, B, and CB functions is
provided by the power law function. This is presented in Tab. 1.

$f(s)$ |  CMF |  SF |  BF |  CBF
---|---|---|---|---
$s^{\mu}$ | $\mu\leq 0$ | $\mu\in[-1,0]$ | $\mu\in[0,1]$ | $\mu\in(0,1)$
$(s+b)^{\mu}$, $b>0$ | $\mu\leq 0$ | $\mu\in[-1,0]$ | $\mu\in[0,1]$ | $\mu\in(0,1)$
Table 1:  Examples of CMF, SF, BF, and CBF RLSchilling10 ; TSandev18a .

###  IV.2 Two theorems of classical mathematical analysis.

Two theorems of classical mathematical analysis are important for the
mathematical toolbox used in the relaxation theory. First of them, which we
abbreviate as Theorem 1, allows one to make analytic continuation of
completely monotone functions to the complex domain. Thus, if its conditions
are satisfied, we can freely jump between the complex variable $z$, e.g.,
$z=\operatorname{i}\\!\omega$ and real $s>0$. The second one is the Efros
theorem AApelblat21 ; Ditkin51 ; Efros35 ; KGorska12a ; UGraf04 ;
VSMartynenko68 ; LWlodarski52  which generalizes the Borel convolution theorem
for the Laplace transform. We will see that if applied to our purposes, the
Efros theorem justifies integral decompositions and can be interpreted as a
guidepost leading to the subordination approach KGorska21 ; KGorska23 . Just
mentioned theorems state as follows.

Characterization of the Laplace transform of a CM, locally integrable function
(ECapelas11 ; GGripenberg90 ) is given by Theorem 1:
The Laplace transform $\widehat{f}(z)$ of a function $f(s)$, $s>0$, that is
locally integrable on $\mathbb{R}_{+}$ and completely monotone, has the
following properties:

  * (a)

$\widehat{f}(z)$ has an analytical extension to the cut complex plane
$\mathbb{C}\setminus\mathbb{R}_{-}$.

  * (b)

$\widehat{f}(s)=\widehat{f}^{\star}(s)$ for $s\in(0,\infty)$,

  * (c)

$\lim_{s\to\infty}\widehat{f}(s)=0$,

  * (d)

$\operatorname{\mathfrak{Im}}[\widehat{f}(z)]<0$ for
$\operatorname{\mathfrak{Im}}(z)>0$,

  * (e)

$\operatorname{\mathfrak{Im}}[z\widehat{f}(z)]\geq 0$ for
$\operatorname{\mathfrak{Im}}(z)>0$ and $\widehat{f}(s)\geq 0$ for
$s\in(0,\infty)$.

Conversely, every function $\widehat{f}(z)$ that satisfies (a)–(c) together
with (d) or (e), is the Laplace transform of a function $f(s)$, which is
locally integrable on $\mathbb{R}_{+}$ and completely monotone on
$(0,\infty)$.

###### Proof.

The proof of this theorem is presented in Ref. GGripenberg90 . ∎

######  Remark 1.

The Laplace transforms of CM functions satisfying the assumptions of Theorem 1
are defined for $z\in\mathbb{C}\setminus\mathbb{R}_{-}$ and because of
conditions (d) or (e) they can be called the Nevanlinna functions Akhiezier ;
Berg ; RLSchilling10 . For positive argument $s$ they are, according to
RLSchilling10 , the Stieltjes functions.

The Efros theorem Ditkin51 ; Efros35 ; UGraf04 ; VSMartynenko68 ; LWlodarski52
. If $\widehat{G}(z)$ and $\widehat{q}(z)$ are analytic functions, and

${\mathscr{L}}[h(\xi);z]=\widehat{h}(z)$ (IV.2)

as well as

${\mathscr{L}}[f(t);z]=\int_{0}^{\infty}f(t)\operatorname{e}^{-zt}\operatorname{d\\!}t=\widehat{G}(z)\operatorname{e}^{-\xi\widehat{q}(z)},$
(IV.3)

then

$\widehat{G}(z)\,\widehat{h}(\widehat{q}(z))=\int_{0}^{\infty}\\!\Big{[}\\!\int_{0}^{\infty}\\!\\!h(\xi)f(\xi,t)\operatorname{d\\!}\xi\Big{]}\operatorname{e}^{-zt}\operatorname{d\\!}t.$
(IV.4)

###### Proof.

The proof can be found in Ref. Efros35 ; VSMartynenko68 . ∎

######  Remark 2.

From Efros theorem, it follows that

${\mathscr{L}}^{-1}[\widehat{G}(z)\widehat{h}(\widehat{q}(z));t]=\int_{0}^{\infty}\\!\\!h(\xi)f(\xi,t)\operatorname{d\\!}\xi.$
(IV.5)

###  IV.3 Some special functions

#### The generalized hypergeometric functions

${{}_{p}F_{q}}\Big{(}{(c_{p})\atop(d_{q})};x\Big{)}$, $x\in\mathbb{R}$, are
defined via the series NIST ; APPrudnikov-v3 :

${{}_{p}F_{q}}\left({(c_{p})\atop(d_{q})};x\right)\operatorname{{\stackrel{{\scriptstyle{\scriptscriptstyle{\mathsf{def}}}}}{{=}}}}\sum_{r=0}^{\infty}\frac{x^{r}}{r!}\frac{(c_{1})_{r}(c_{2})_{r}\cdots(c_{p})_{r}}{(d_{1})_{r}(d_{2})_{r}\cdots(d_{q})_{r}},$
(IV.6)

where $(c)_{r}=\Gamma(c+r)/\Gamma(c)=c(c+1)\cdots(c+r-1)$ is the Pochhammer
symbol (rising factorial).

#### The Meijer G-functions.

The Meijer G-functions NIST ; APPrudnikov-v3  are defined through the Mellin-
Barnes integrals introduced as inverses in the sense of the Mellin transform
of expressions being ratios of certain $\Gamma$ functions. The definition of
the Meijer G-function reads

$\displaystyle
G^{m,n}_{p,q}\left(x\Big{|}{(a_{p})\atop(b_{q})}\right)\operatorname{{\stackrel{{\scriptstyle{\scriptscriptstyle{\mathsf{def}}}}}{{=}}}}$
(IV.7)
$\displaystyle\int_{\gamma_{L}}\frac{\prod_{i=1}^{m}\Gamma(b_{i}+s)\prod_{i=1}^{n}\Gamma(1-a_{i}-s)}{\prod_{i=n+1}^{p}\Gamma(a_{i}+s)\prod_{i=m+1}^{q}\Gamma(1-b_{i}-s)}\,x^{-s}\frac{\operatorname{d\\!}s}{2\pi\\!\operatorname{i}}.$

Parameters in Eq. (IV.7) are subject to conditions

$\displaystyle\begin{split}&x\neq 0,\quad 0\leq m\leq q,\quad 0\leq n\leq
p;\\\ &a_{i}\in\mathbb{C},\quad i=1,\ldots,p;\qquad b_{i}\in\mathbb{C},\quad
i=1,\ldots,q;\\\
&(a_{p})=a_{1},a_{2},\ldots,a_{p};\quad(b_{q})=b_{1},b_{2},\ldots,b_{q}.\end{split}$

The contour $\gamma_{L}$ lies between the poles of $\Gamma(1-a_{i}-s)$ for
$i=1,\ldots,n$ and the poles of $\Gamma(b_{i}+s)$ for $i=1,\ldots,m$. For
generalization of Eq. (IV.7) to the Fox H-function, see (APPrudnikov-v3, , Eq.
(8.3.2.21)) namely

$G^{m,n}_{p,q}\left(z\Big{|}{(a_{p})\atop(b_{q})}\right)=H^{m,n}_{p,q}\left[z\Big{|}{[a_{p},1]\atop[b_{q},1]}\right].$

Another possibility to express the Meijer G-function in terms of Fox
H-function is provided by (APPrudnikov-v3, , Eq. (8.3.2.22)). It is obtained
by using the Gauss–Legendre formula for the Gamma function in Eq. (IV.7):

$\prod_{j=0}^{n-1}\Gamma\left(z+\frac{j}{n}\right)=(2\pi)^{\frac{n-1}{2}}n^{\frac{1}{2}-nz}\Gamma(nz),$
(IV.8)

for $z\neq 0,-1,-2,\ldots$, and $n=1,2,\ldots$.

The main properties of Meijer G-functions are as follows:

  * (i)

They form a set closed under the Laplace transform. It means that its (direct
and/or inverse) Laplace transform leads to another Meijer G-function:

$\mathscr{L}\left[x^{\alpha-1}G^{m,n}_{p,q}\left(\omega
x^{l}\Big{|}{(a_{p})\atop(b_{q})}\right);z\right]=\\\
\frac{l^{\alpha-\frac{1}{2}}z^{-\alpha}}{(2\pi)^{\frac{l-1}{2}}}G^{m,n+l}_{p+l,q}\left(\frac{\omega
l^{l}}{z^{l}}\Big{|}{\Delta(l,1-\alpha),(a_{p})\atop(b_{q})}\right),$ (IV.9)

and

$\mathscr{L}^{-1}\left[z^{-\alpha}G^{m,n}_{p,q}\left(\frac{\omega}{z^{l}}\Big{|}{(a_{p})\atop(b_{q})}\right);x\right]=\\\
\frac{l^{1/2-\alpha}x^{\alpha-1}}{(2\pi)^{\frac{1-l}{2}}}G^{m,n}_{p,q+l}\left(\frac{\omega
x^{l}}{l^{l}}\Big{|}{(a_{p})\atop(b_{a}),\Delta(l,1-\alpha)}\right).$ (IV.10)

The conditions on parameters under which Eqs. (IV.9) and (IV.10) are held are
given in (APPrudnikov-v3, , Eq. (2.24.3.1)) and (APPrudnikov-v5, , Eq.
(3.38.1.2)), respectively. The symbol $\Delta(n,a)$ denotes the special list
of parameters which reads $a/n,(a+1)/n,\ldots,(a+n-1)/n$.

  * (ii)

Below we quote two identities which will be used in our calculations:

$\displaystyle G^{m,n}_{p,q}\left(z\Big{|}{(a_{p})\atop(b_{q})}\right)$
$\displaystyle=G^{n,m}_{q,p}\left(\frac{1}{z}\Big{|}{(1-b_{q})\atop(1-a_{p})}\right),$
(IV.11) $\displaystyle
z^{\mu}\,G^{m,n}_{p,q}\left(z\Big{|}{(a_{p})\atop(b_{q})}\right)$
$\displaystyle=G^{m,n}_{p,q}\left(z\Big{|}{(a_{p}+\mu)\atop(b_{q}+\mu)}\right),$
(IV.12)

see NIST ; APPrudnikov-v3 .

  * (iii)

Important procedure often applied in calculations is to convert the Meijer
G-function $G^{m,n}_{p,q}$ into the finite sum of the generalized
hypergeometric functions ${{}_{p}F_{q}}$. This formula can be found in, e.g.,
(APPrudnikov-v3, , Eq. (8.2.2.3)) and reads

$G^{m,n}_{p,q}\left(z\Big{|}{(a_{p})\atop(b_{q})}\right)=\sum_{j=1}^{m}\left[\prod_{i=0}^{j-1}\Gamma(b_{i}-b_{j})\prod_{i=j+1}^{m}\Gamma(b_{i}-b_{j})\right]\\\
\times\frac{[\prod_{i=0}^{n}\Gamma(1+b_{j}-a_{i})]}{[\prod_{i=n+1}^{p}\Gamma(a_{i}-b_{j})]\,[\prod_{i=m+1}^{q}\Gamma(1+b_{j}-b_{j})]}z^{b_{j}}\\\
\times{{}_{p+1}F_{q}}\left({1,1+b_{j}-(a_{p})\atop
1+b_{j}-(b_{q})};(-1)^{p-m-n}z\right),$ (IV.13)

where $p\leq q$, $b_{i}-b_{j}\neq 0,\pm 1,\pm 2,\ldots$, $i\neq j$,
$i,j=1,2,\ldots,m$. The case $p\geq q$ can be obtained from Eqs. (IV.11) and
(IV.13). The conditions under which Eq. (IV.13) is held may be found in
(APPrudnikov-v3, , Eqs. (8.2.2.3)).

#### The family of the Mittag-Leffler functions.

The three parameter generalization of the Mittag-Leffler function book1  in
the literature devoted to the fractional calculus and its applications is
known as the Prabhakar function. It was introduced by T. R. Prabhakar
Prabhakar  and reads:

$E_{\alpha,\mu}^{\nu}(z)\operatorname{{\stackrel{{\scriptstyle{\scriptscriptstyle{\mathsf{def}}}}}{{=}}}}\sum_{r\geq
0}\frac{(\nu)_{r}}{r!\,\Gamma(\alpha r+\mu)}z^{r},\qquad
z\in\mathbb{C},\quad\operatorname{\mathfrak{Re}}(\alpha)>0,$ (IV.14)

where the Pochhammer symbol $(\nu)_{r}$ is given below Eq. (IV.6). For
$\mu=\nu=1$ it reduces to the (standard, called also one-parameter) Mittag-
Leffler function $E_{\alpha}(z)$ being the relaxation function of the CC
model, which was first considered in Mittag2 ; Mittag0 ; Mittag1 . For $\nu=1$
it becomes the two-parameters Mittag-Leffler function $E_{\alpha,\mu}(z)$
studied in 1905 by Wiman Wiman .

######  Remark 3.

From the definition (IV.14) it follows that
$E_{\alpha,1}^{1}(-x^{\alpha})\equiv
E_{\alpha}(-x^{\alpha})=1-x^{\alpha}E_{\alpha,1+\alpha}(-x^{\alpha})=\delta(x)-x^{-1}E_{\alpha,0}(-x^{\alpha})$,
where introducing the $\delta$-Dirac function allows us to avoid the pole at
$x=0$.

The comprehensive review of the family of the Mittag-Leffler functions can be
found in RGarra18 ; book1 ; HJHaubold11 . Here, we shall quote a few formulae
employed throughout the paper.

######  Lemma 1.

The asymptotics of $E_{\alpha,\mu}^{\nu}(-x^{\alpha})$ for small and large $x$
is

$E_{\alpha,\mu}^{\nu}(-x^{\alpha})\propto\frac{1}{\Gamma(\mu)}-\frac{\nu\,x^{\alpha}}{\Gamma(\alpha+\mu)}\qquad\text{for}\quad
x\ll 1$ (IV.15)
$E_{\alpha,\mu}^{\nu}(-x^{\alpha})\propto\frac{x^{-\alpha\nu}}{\Gamma(\mu-\alpha\nu)}-\frac{\nu
x^{-\alpha\nu-\alpha}}{\Gamma(\mu-\alpha-\alpha\nu)}\qquad\text{for}\quad x\gg
1.$ (IV.16)

###### Proof.

Asymptotics for $x\ll 1$ comes from the series expansion of three-parameters
Mittag-Leffler function given by Eq. (IV.14). The justification of Eq. (IV.16)
is shown in RGarra18 . ∎

For $x\in\mathbb{R}$ and $\alpha=l/k$ such that $0<l<k$, the three-parameters
Mittag-Leffler function $E_{l/k,\mu}^{\nu}(x)$ can be represented as a finite
sum of $k$ generalized hypergeometric functions KGorska20a :

$E_{l/k,\mu}^{\nu}(x)=\sum_{j=0}^{k-1}\frac{(\nu)_{j}\,x^{j}}{j!\,\Gamma(\mu+\frac{l}{k}j)}\\\
\times{{}_{1+k}F_{l+k}}\left({1,\Delta(k,\nu+j)\atop\Delta(k,1+j),\Delta(l,\mu+\frac{l}{k}j)};\frac{x^{k}}{l^{l}}\right),$
(IV.17)

where $\Delta(n,a)$ is a special list of parameters defined below Eq. (IV.10).

The Meijer G-representation of $E_{\alpha,\mu}^{\nu}(x)$ for $\alpha=l/k$ can
be obtained by employing the Laplace integral of
$E_{\alpha,\mu}^{\nu}(-at^{\alpha})$. The latter is given by (KGorska21b, ,
Eqs. (4), (5), and (8)) and/or (KGorska21c, , Eqs. (A4) and (A5))

$E_{\alpha,\mu}^{\nu}(x)=\alpha^{-1}\int_{0}^{\infty}\operatorname{e}^{x\xi}\xi^{-1-1/\alpha}\varPhi_{\alpha,\mu}^{\nu}(\xi^{-1/\alpha})\operatorname{d\\!}\xi,$
(IV.18)

in which

$\varPhi_{l/k,\mu}^{\nu}(y)=\frac{k^{\nu-1}}{l^{\mu-1}\Gamma(\nu)}\frac{\sqrt{lk}}{(2\pi)^{(k-l)/2}}\frac{1}{y}G^{k,0}_{l,k}\left(\frac{l^{l}}{k^{k}y^{l}}\Big{|}{\Delta(l,\mu)\atop\Delta(k,\nu)}\right).$
(IV.19)

We can distinguish at least two special cases of Eq. (IV.18) involving the
one-sided Lévy stable distribution
$\varPhi_{\alpha}(y)\equiv\varPhi_{\alpha,1}^{1}(y)$ with $\alpha\in(0,1)$ and
$y\in\mathbb{R}_{+}$ related by the Laplace transform to the stretched
exponential function (the KWW relaxation) KAPenson10 ; HPollard46 , see two
examples below.

######  Example 1.

The first of them is obtained from Eq. (IV.19) by taking $\mu=\nu=1$. Due to
the definition of the Meijer G-function in the numerator of (IV.7) the only
non-vanishing terms are equal to $\prod_{i=1}^{m}\Gamma(b_{i}+s)$ whereas in
the denominator survive only the terms $\prod_{i=n+1}^{p}\Gamma(a_{i}+s)$
which read $\Delta(l,1)=\frac{1}{l},\frac{2}{l},\ldots,\frac{l-1}{l},1$ and
$\Delta(k,1)=\frac{1}{k},\frac{2}{k},\ldots,\frac{k-1}{k},1$, respectively.
Delating “$1$” in both these sequences and adding “$0$” allows us to shift the
upper and lower lists to $\Delta(l,0)=0,\frac{1}{l},\ldots,\frac{l-1}{l}$ and
$\Delta(k,0)=0,\frac{1}{k},\ldots,\frac{k-1}{k}$. In consequence, we get the
integral representation of the Mittag-Leffler function first obtained by H.
Pollard in Ref. HPollard48  and independently derived in Refs. KGorska12 ;
KWeron96 :

$E_{\alpha,1}^{1}(-x)\equiv
E_{\alpha}(-x)=\int_{0}^{\infty}\operatorname{e}^{-x\xi}\xi^{-1-1/\alpha}\varPhi_{\alpha}(\xi^{-1/\alpha})\frac{\operatorname{d\\!}\xi}{\alpha}.$

######  Example 2.

As the second example, we take $\mu=\alpha\beta$ and $\nu=\beta$. Then
$\varPhi_{\alpha,\alpha\beta}^{\beta}(y)=\alpha
y^{-\alpha\beta}\varPhi_{\alpha}(y)/\Gamma(\beta)$ for which we have
(KGorska18, , Eq. (11)), this is

$E_{\alpha,\alpha\beta}^{\beta}(-x)=\int_{0}^{\infty}\operatorname{e}^{-x\xi}\frac{\xi^{\beta-(1+\frac{1}{\alpha})}}{\Gamma(\beta)}\varPhi_{\alpha}(\xi^{-\frac{1}{\alpha}})\operatorname{d\\!}\xi.$

The substitution of Eq. (IV.19) into Eq. (IV.18) enables us to calculate the
integral in Eq. (IV.10) and express it in the form

$\displaystyle\begin{split}E_{l/k,\mu}^{\nu}(-x)&=\frac{k^{\nu}}{\Gamma(\nu)}\frac{l^{\frac{l}{k}\nu-\mu+\frac{1}{2}}}{(2\pi)^{k-(1+l)/2}}\,\frac{1}{x^{\nu}}\\\
&\times
G^{k,k}_{k+l,k}\left(\frac{l^{l}}{x^{k}}\Big{|}{\Delta(k,1-\nu),\Delta(l,\mu-\frac{l}{k}\nu)\atop\Delta(k,0)}\right)\\\
&=\frac{k^{\nu}}{\Gamma(\nu)}\frac{l^{\frac{1}{2}-\mu}}{(2\pi)^{k-(1+l)/2}}\\\
&\times
G^{k,k}_{k+l,k}\left(\frac{l^{l}}{x^{k}}\Big{|}{\Delta(k,1),\Delta(l,\mu)\atop\Delta(k,\nu)}\right).\end{split}$
(IV.20)

Employing Eqs. (IV.11) and (IV.13) we can express Eq. (IV.20) in terms of the
generalized hypergeometric functions, vis. Eq. (IV.17).

The Laplace transform of the three-parameters Mittag-Leffler function
(Prabhakar, , Eq. (2.5)):

$\mathscr{L}[t^{\mu-1}E_{\alpha,\mu}^{\nu}(-\lambda
t^{\alpha});z]=\frac{z^{\alpha\nu-\mu}}{(\lambda+z^{\alpha})^{\nu}},\qquad\lambda\in\mathbb{R},$
(IV.21)

is crucially important for solving problems which appear in the relaxation
theory, e.g., it allows us to calculate the response function $\phi(t)$ from
the experimental data encoded in phenomenological spectral functions
$\widehat{\phi}(s)$ whose forms, see Eq. (I.4), match the RHS of Eq. (IV.21).
Both heuristic, and rigorous arguments quoted in Sect. II suggest that looking
for the relaxation functions we should focus our attention on the functions
belonging to the CM class. Here it has to be said that physicists’ interest in
Mittag-Leffler functions as elements of the CM class is by no means limited to
dielectric relaxation phenomena ECapelas11 ; RGarrappa16 ; Hanyga1 ;
FMainardi15 ; Tomovski14  but is also well-grounded in studies of other
phenomena exhibiting the memory-dependent time evolution, to mention
hereditary mechanics and viscoelasticity Hanyga2 ; Mainardibook ; Rabotnov ;
Rossikhin3 ; Rossikhin2 . Mathematicians’ interest to study CM property of the
Mittag-Leffler functions stems from the classical works of H. Pollard and W.
R. Schneider who proved it for $E_{\alpha}(-x)$ HPollard48  and
$E_{\alpha,\beta}(-x)$ WRSchneider96 . Currently the research in the field is
inspired by problems and applications of the special functions theory Samko ;
Shukla ; Simon , fractional calculus and fractional differential equations
book1 ; Mainardibook . Recent results include proofs of the CM character
obeyed by the Prabhakar function FMainardi15 , see the lemma below, and by
$E_{\alpha,\mu}^{\nu}(-at^{\alpha})$ itself if $0<\alpha\leq 1$,
$\mu\geq\alpha\nu$, and $\gamma>0$, as shown in KGorska21b . CM property of
$E_{\alpha,\mu}^{\nu}(-at^{\alpha})$ explains restriction $\beta<1$ put on the
parameter $\beta$ in the Prabhakar function to preserve its CM as a result of
the multiplication of two CM functions.

######  Lemma 2.

$t^{\mu-1}E_{\alpha,\mu}^{\nu}(-at^{\alpha})$ is CM function for $0<\alpha\leq
1$, $0<\mu\leq\alpha\nu<1$, and $\nu>0$.

###### Proof.

Different proofs of this lemma were presented in KGorska21b ; FMainardi15 ;
Tomovski14 . ∎

######  Lemma 3.

$x^{-\alpha\beta-1}E^{\beta}_{-\alpha,-\alpha\beta}(-\lambda x^{-\alpha})$ can
be replaced by
$\lambda^{-\beta}x^{-1}E_{\alpha,0}^{\beta}(-x^{\alpha}/\lambda)$,
$\lambda\in\mathbb{R}$.

Proof. Lemma 3 comes from direct calculations:

$\displaystyle\mathscr{L}[x^{-\alpha\beta-1}$ $\displaystyle
E_{-\alpha,-\alpha\beta}^{\beta}(-\lambda x^{-\alpha});z]$
$\displaystyle=\frac{1}{(z^{-\alpha}+\lambda)^{\beta}}=\lambda^{-\beta}\frac{z^{\alpha\beta}}{(\lambda^{-1}+z^{\alpha})^{\beta}}$
$\displaystyle=\mathscr{L}^{-1}[\lambda^{-\beta}x^{-1}E_{\alpha,0}^{\beta}(-x^{\alpha}/\lambda);z].\qquad\qed$

######  Lemma 4.

$E^{\beta}_{\alpha,1}(-\lambda
x^{\alpha})=\lambda^{-\beta}\int_{0}^{x}u^{-\alpha\beta-1}E_{-\alpha,-\alpha\beta}^{\beta}(-u^{-\alpha}/\lambda)\operatorname{d\\!}u,$

$\lambda\in\mathbb{R}$.

Proof. Lemma 4 can be shown by making the direct calculations analogical like
in the proof of Lemma 3 and using
$\mathscr{L}^{-1}[z^{-1}f(z);x]=\int_{0}^{x}\mathscr{L}^{-1}[f(z);u]\operatorname{d\\!}u$.
∎

We also mention that the fractional derivative in Riemann-Liouville sense of
$x^{\mu-1}E_{\alpha,\mu}^{\nu}(\lambda x^{\alpha})$ yields (book1, , Eq.
(5.1.34)), namely

$\left(D^{\beta}[y^{\mu-1}E_{\alpha,\mu}^{\nu}(\lambda
y^{\alpha})]\right)\\!(x)=x^{\mu-\beta-1}E_{\alpha,\mu-\beta}^{\nu}(\lambda
x^{\alpha}).$ (IV.22)

##  V The Havriliak-Negami model and its physical content

In this section, we shall illustrate the general theoretical approach
presented in Secs. I-III on the example of the HN relaxation model and its
special cases, like the Debye, CC and CD models. The HN relaxation pattern was
introduced in SHavriliak67  to parametrize experimental data describing the
frequency dependence of the complex dielectric permittivity
$\widehat{\varepsilon^{\star}}(\omega)$ measured in polymers. Despite its
purely phenomenological origin and apparent simplicity, the applicability of
the HN model went far beyond its initial implementation. The model has
appeared well-working for a much larger plethora of dielectric phenomena and
has become the “first choice” method to analyse experimental relaxation data
for various relaxation phenomena taking place in different complex systems, by
no mere limited to those of condensed matter physics and materials science.
Unexpected examples of its utility include the use of the HN function for
monitoring the contamination in sandstone VSaltas07 , or investigations of
complex systems representing plant tissues of fresh fruits and vegetables, for
which the HN relaxation in the frequency range $10^{7}-1.8\times 10^{9}$ Hz
was shown to be an useful tool of analysis RRNigmatulin06 .

###  V.1 The spectral function
$\widehat{\phi}_{H\\!N}(\alpha,\beta;\operatorname{i}\\!\omega)$

We start our consideration by recalling the spectral function for the HN model

$\widehat{\phi}_{H\\!N}(\alpha,\beta;\operatorname{i}\\!\omega)=\left[\frac{\widehat{\varepsilon^{\star}}(\omega)-\varepsilon_{\infty}}{\varepsilon-\varepsilon_{\infty}}\right]_{H\\!N}=[1+(\operatorname{i}\\!\omega\tau)^{\alpha}]^{-\beta},$
(V.1)

$\alpha,\beta\in(0,1]$, given by Eqs. (I.4). The range of $\alpha$ and $\beta$
parameters (called respectively the width and asymmetry) are indicated by the
experiment which implies that they belong to $(0,1]$. We recall that for
$\alpha=\beta=1$ we have the Debye spectral function, for $\alpha\in(0,1)$ and
$\beta=1$ it reduces to the CC case, and for $\alpha=1$ and $\beta\in(0,1)$
becomes the CD model.

Separating real and imaginary parts of
$\widehat{\varepsilon^{\star}}_{H\\!N}(\omega)=\widehat{\varepsilon^{\prime}}_{H\\!N}(\omega)-\operatorname{i}\widehat{\varepsilon^{\prime\prime}}_{H\\!N}(\omega)$
we get after some complex algebra (see Ref. CJFBottcher78 ):

$\widehat{\varepsilon^{\prime}}_{H\\!N}(\omega)=\varepsilon_{\infty}+\frac{(\varepsilon-\varepsilon_{\infty})\cos\big{[}\beta\,\theta_{\alpha/2}\big{(}(\omega
t)^{2}\big{)}\big{]}}{\big{[}1+2(\omega\tau)^{\alpha}\cos(\pi\alpha/2)+(\omega\tau)^{2\alpha}\big{]}^{\beta/2}},$

and

$\widehat{\varepsilon^{\prime\prime}}_{H\\!N}(\omega)=\frac{(\varepsilon-\varepsilon_{\infty})\sin\big{[}\beta\,\theta_{\alpha/2}\big{(}(\omega
t)^{2}\big{)}\big{]}}{\big{[}1+2(\omega\tau)^{\alpha}\cos(\pi\alpha/2)+(\omega\tau)^{2\alpha}\big{]}^{\beta/2}},$

where

$\theta_{\alpha}(y)=\arctan\left(\frac{\sin(\pi\alpha)}{y^{-\alpha}+\cos(\pi\alpha)}\right)$
(V.2)

provides restrictions on the parameters of the model FMainardi15 .

###  V.2 The response function $\phi_{H\\!N}(\alpha,\beta;t)$

To find the response function $\phi_{H\\!N}(\alpha,\beta;t)$ firstly we set
$z=\operatorname{i}\\!\omega$ in Eq. (V.1) and calculate the inverse Laplace
transform of $\widehat{\phi}_{H\\!N}(\alpha,\beta;z)$. Then, with the help of
Eq. (IV.21) we get

$\displaystyle\begin{split}\phi_{H\\!N}(\alpha,\beta;t)&=\tau^{-\alpha\beta}\mathscr{L}^{-1}[(\tau^{-\alpha}+s^{\alpha})^{-\beta};t]\\\
&=\frac{1}{\tau}\,\big{(}\mbox{$\frac{t}{\tau}$}\big{)}^{\alpha\beta-1}E^{\beta}_{\alpha,\alpha\beta}\big{(}\\!-\big{(}\mbox{$\frac{t}{\tau}$}\big{)}^{\alpha}\big{)},\end{split}$
(V.3)

where $E_{\alpha,\alpha\beta}^{\beta}(-(t/\tau)^{\alpha})$ is the three-
parameter Mittag-Leffler function described in Sect. IV.3. For
$\alpha=\beta=1$ it reduces to the Debye response function
$\phi_{H\\!N}(1,1;t)\equiv\phi_{D}(t)=\exp(-t/\tau)/\tau$ while for $\alpha=1$
and $\beta\in(0,1)$ we get the CD response function

$\displaystyle\begin{split}\phi_{H\\!N}(1,\beta;t)\equiv\phi_{C\\!D}(\beta;t)&=\tau^{-1}\big{(}\mbox{$\frac{t}{\tau}$}\big{)}^{\beta-1}E_{1,\beta}^{\beta}\big{(}\\!-\\!\big{(}\mbox{$\frac{t}{\tau}$})\big{)}\\\
&=\frac{(t/\tau)^{\beta-1}}{\tau\Gamma(\beta)}\operatorname{e}^{-t/\tau},\end{split}$

and for $\alpha\in(0,1)$ and $\beta=1$ the CC response function

$\phi_{H\\!N}(\alpha,1;t)\equiv\phi_{CC}(\alpha;t)=\tau^{-1}\big{(}\mbox{$\frac{t}{\tau}$}\big{)}^{\alpha-1}E_{\alpha,\alpha}\big{(}\\!-\\!\big{(}\mbox{$\frac{t}{\tau}$}\big{)}^{\alpha}\big{)}.$

For rational $\alpha=l/k$ the above formulae can be expressed in terms of
special functions directly implemented in the CAS, like the Meijer G-function
and/or the finite sum of the generalized hypergeometric functions
${{}_{p}F_{q}}$:

$\displaystyle\phi_{H\\!N}$
$\displaystyle(l/k,\beta;t)=(2\pi)^{\mbox{$\frac{1+l}{2}$}-k}\frac{\sqrt{l}k^{\beta}}{\Gamma(\beta)}\,\frac{1}{t}$
(V.4) $\displaystyle\times
G^{k,k}_{l+k,k}\left(\frac{l^{l}\tau^{l}}{t^{l}}\Big{|}{\Delta(k,1-\beta),\Delta(l,0)\atop\Delta(k,0)}\right)$
$\displaystyle=\frac{1}{\tau}\sum_{j=0}^{k-1}\frac{(-1)^{j}(\beta)_{j}}{j!\,\Gamma(\frac{l}{k}(\beta+j))}\left(\frac{t}{\tau}\right)^{\\!\\!\frac{l}{k}(\beta+j)-1}$
(V.5)
$\displaystyle\times{{}_{1+l}F_{l+k}}\left({1,\Delta(k,\beta+j)\atop\Delta(k,1+j),\Delta(l,\frac{l}{k}(\beta+j))};\frac{(-1)^{k}t^{l}}{l^{l}\tau^{l}}\right).$

In Fig. 1 we present $\phi_{H\\!N}(\alpha,\beta;t)$ for $\alpha=1/2$ and
$\beta=1/2,1,2$ as well as for $\beta=3$.

Figure 1: (Color online) Plot of $\phi_{H\\!N}(1/2,\beta;t)$ given by Eq.
(V.4) for $\tau=1$ and $\beta=1/2$ (the red continuous curve no. I), $\beta=1$
(the blue dashed curve no. II), $\beta=2$ (the black dash-dotted curve no.
III), and $\beta=3$ (the green dotted curve no. IV).

It is seen that for $\beta=3$, $\phi_{H\\!N}(1/2,3;t)$ is unimodal with the
maximum at $t=t_{\rm max}$. Thus it can not be CM function. Looking at plots
in Fig. 1, one may see the candidates to be a CM function are
$\phi_{H\\!N}(1/2,\beta;t)$ with $\beta\leq 2=1/\alpha$. This is confirmed by
Lemma 3 applied for Eq. (V.3) but the case III has to be excluded as
experimental data restrict values of $\beta$ to the interval $(0,1]$.

For reader’s convenience we quote the asymptotic behavior of response function
$\phi_{H\\!N}(\alpha,\beta;t)$ for a small and large time presented also in
(RGarrappa16, , Eq. (3.29)):

$\phi_{H\\!N}(\alpha,\beta;t)\propto\frac{1}{\tau\Gamma(\alpha\beta)}\left(\frac{t}{\tau}\right)^{\alpha\beta-1}\quad\text{for}\quad\frac{t}{\tau}\ll
1\quad\text{and}\\\
\phi_{H\\!N}(\alpha,\beta;t)\propto-\frac{\beta}{\tau\Gamma(-\alpha)}\left(\frac{t}{\tau}\right)^{-1-\alpha}\quad\text{for}\quad\frac{t}{\tau}\gg
1.$ (V.6)

The asymptotics at $t\to 0$ depends on the value of $\beta$ namely
$\phi_{H\\!N}(\alpha,\beta;t)$ tends to infinity for $\beta<1/\alpha$ (note
that in such a case the singularity is integrable), goes to $\tau^{-1}$ for
$\beta=1/\alpha$, and approaches 0 for $\beta>1/\alpha$.

###  V.3 The probability density $g_{H\\!N}(\alpha,\beta;\xi)$

Thanks to Theorem 1 (Sect. IV.2) it turns out that the spectral function has
the integral representation as in Eq. (II.9):

$\widehat{\phi}_{H\\!N}(\alpha,\beta;z)=\int_{0}^{\infty}\frac{p_{H\\!N}(\alpha,\beta;\xi)}{\xi+z}\operatorname{d\\!}\xi,$
(V.7)

where $p_{H\\!N}(\alpha,\beta;\xi)$ is PDF. Restricting the complex $z$ onto
$s\in\mathbb{R}_{+}$ we observe that $\widehat{\phi}_{H\\!N}(\alpha,\beta;z)$
is a S function on the real axis and, thus, we confirm the information
presented in Ref. KGorska21c . From the Schwinger parametrization it follows
that Eq. (V.7) is the Laplace transform of the Laplace transform, namely
$\widehat{\phi}_{H\\!N}(\alpha,\beta;z)=\mathscr{L}[\phi_{H\\!N}(\alpha,\beta;t);z]$
where
$\phi_{H\\!N}(\alpha,\beta;t)=\mathscr{L}[p_{H\\!N}(\alpha,\beta;\xi);t]$.
Because $\phi_{H\\!N}(\alpha,\beta;t)$ is a CM function for $\alpha\in(0,1]$
and $\beta\leq 1/\alpha$, then from the Bernstein theorem it emerges that
$p_{H\\!N}(\alpha,\beta;\xi)$ is a non-negative function for

$0<\alpha<1\quad\text{and}\quad 0<\beta\leq 1/\alpha.$ (V.8)

Furthermore, from Eq. (II.6) and
$\phi_{H\\!N}(\alpha,\beta;t)={\mathcal{L}}[p_{H\\!N}(\alpha,\beta;\xi);t]$ we
get

$\displaystyle\begin{split}g_{H\\!N}(\alpha,\beta;\xi)&=\frac{\tau}{\xi}\,p_{H\\!N}(\alpha,\beta;\xi)\\\
&=\frac{\tau}{\xi}\,\mathscr{L}^{-1}[\phi_{H\\!N}(\alpha,\beta;t);\xi],\end{split}$
(V.9)

being also a non-negative function, since the product of the non-negative
functions is also non-negative. That statement agrees with one of the results
presented in (KGorska18, , Sect. 4).

For the rational $\alpha$ such that $0<\alpha=l/k<1$, we substitute
$\phi_{H\\!N}(\alpha,\beta;t)$ in Eq. (V.4) and employ Eq. (IV.10). Hence, we
conclude that $g_{H\\!N}(l/k,\beta;\xi)$ can be written as

$\displaystyle\begin{split}g_{H\\!N}(&l/k,\beta;\xi)=(2\pi)^{\frac{1+l}{2}-k}\frac{\sqrt{l}k^{\beta}}{\Gamma(\beta)}\,\frac{1}{\xi}\\\
&\times\mathscr{L}^{-1}\left[\frac{\tau}{t}G^{k,k}_{l+k,k}\left(\frac{l^{l}\tau^{l}}{t^{l}}\Big{|}{\Delta(k,1-\beta),\Delta(l,0)\atop\Delta(k,0)}\right);\xi\right]\\\
&=(2\pi)^{l-k}\frac{k^{\beta}}{\Gamma(\beta)}\,\frac{1}{\xi}\,G^{k,k}_{l+k,l+k}\left(\xi^{l}\Big{|}{\Delta(k,1-\beta),\Delta(l,0)\atop\Delta(k,0),\Delta(l,0)}\right),\end{split}$
(V.10)

with the symbols $\Delta(n,a)$ defined below Eq. (IV.17). Next, using Eqs.
(IV.11) and (IV.13) as well as the Gauss-Legendre multiplication formula for
$\Gamma$ functions in Eq. (V.10), we can express $g_{H\\!N}(l/k,\beta;\xi)$ as
a finite sum of $k$ generalized hypergeometric functions

$\displaystyle\begin{split}g_{H\\!N}(l/k,\beta;\xi)&=\frac{1}{\pi}\sum_{j=0}^{k-1}\frac{(-1)^{j}}{j!}\frac{(\beta)_{j}}{\xi^{1+\frac{l}{k}n_{j}}}\sin(\mbox{$\frac{l}{k}$}n_{j}\pi)\\\
&\times{{}_{k+1}F_{k}}\left({1,\Delta(k,n_{j})\atop\Delta(k,1+j)};\frac{(-1)^{l-k}}{\xi^{l}}\right)\end{split}$
(V.11)

with $n_{j}=\beta+j$. Equation (V.11) gives the form of
$g_{H\\!N}(l/k,\beta;\xi)$ which is convenient and efficiently applicable in
calculations using the standard computer algebra packages. For example, for
the CC relaxation ($\beta=1$) it is seen that, due to appropriate
cancellations, ${{}_{k+1}F_{k}}$’s reduce to ${{}_{1}F_{0}}\left({1\atop
0},u\right)$, which can be further simplified to $(1-u)^{-1}$ where
$u=(-1)^{l-k}\xi^{-l/k}$. Employing it and Eqs. (1.353.1) and (1.353.3) on p.
38 of Gradshteyn07  to the sum over $j$, we get

$\displaystyle\begin{split}g_{H\\!N}(\alpha,1;\xi)&\equiv
g_{CC}(\alpha;\xi)\\\
&=\frac{\xi^{\alpha-1}\sin(\alpha\pi)}{\pi(\xi^{2\alpha}+2\xi^{\alpha}\cos(\alpha\pi)+1)},\end{split}$
(V.12)

with $0<\alpha=l/k<1$. We point out that Eq. (V.12) was obtained in
ECapelasDeOliveira14 ; BDybiec10 ; RGorenflo08 ; KWeron96  using different
methods, see Eq. (3.24) on p. 245 in RGorenflo08 , Eqs. (22) and (39) in
ECapelasDeOliveira14  or Eq. (26) in BDybiec10 . Distributions
$g_{CC}(\alpha;\xi)$ are non-negative functions for $\alpha\in(0,1)$ and they
share the following properties: (i) $g_{CC}(\alpha;\xi)$ are non-negative and
non-increasing for $0<\alpha<1$ and $\xi\geq 0$; and (ii) $g_{CC}(\alpha;\xi)$
go to infinity at $\xi=0$, and vanish for $\xi\to\infty$. For the CD
relaxation ($\alpha=1$) in Eq. (V.11) we take $l=k=1$ which leads to
(RGarrappa16, , Eq. (3.19)):

$g_{H\\!N}(1,\beta;\xi)\equiv
g_{C\\!D}(\beta;\xi)=\frac{\sin(\beta\pi)}{\pi\xi(\xi-1)^{\beta}}\,\Theta(\xi-1),$

with the Heaviside step function $\Theta(\cdot)$ which guarantees that
$g_{C\\!D}(\beta;\xi)$ is real for $\xi\leq 1$.

Our last task in this subsection is to find the series representation of
$g_{H\\!N}(\alpha,\beta;\xi)$ for $\alpha=l/k$. To achieve this goal, we use
the series representation of the generalized hypergeometric functions
${{}_{p}F_{q}}$, see Eq. (IV.6). Following such a way we get Eq. (V.11) as the
double sum: one over $r$ ($r=0,1,2,\ldots$) which comes from the series
representation of ${{}_{p}F_{q}}$, and the second one over $j$
($j=0,1,\ldots,k-1$) which appears in the Eq. (V.11) itself. Changing the
summation index $kr+j\to n$ we arrive at the expression

$g_{H\\!N}(\alpha,\beta;\xi)=\frac{1}{\pi}\,\sum_{n=0}^{\infty}\frac{(-1)^{n}}{n!}\frac{(\beta)_{n}}{\xi^{1+\alpha(\beta+n)}}\sin[\alpha(\beta+n)\pi],$
(V.13)

which, after representing the sine function as an imaginary part of
$\exp[\operatorname{i}\\!\pi\alpha(\beta+n)]$, using the integral
representation of the $\Gamma$ function, and applying Eq. (IV.14), leads us to
the integral form of Eq. (V.13):

$\displaystyle\begin{split}g_{H\\!N}(\alpha,\beta;\xi)&=\frac{1}{\pi}\operatorname{\mathfrak{Im}}\left\\{\operatorname{e}^{\operatorname{i}\\!\alpha\beta\pi}\int_{0}^{\infty}e^{-u\xi}u^{\alpha\beta}\right.\\\
&\times\left.E_{\alpha,\alpha\beta}^{\beta}(-u^{\alpha}\operatorname{e}^{\operatorname{i}\\!\alpha\pi})\operatorname{d\\!}u\right\\}.\end{split}$
(V.14)

Applying Eq. (IV.21) to Eq. (V.14) and employing de Moivre’s formula to
calculate the imaginary part, we rederive the function
$g_{H\\!N}(\alpha,\beta;\xi)$ in the form obtained and extensively studied in
the just quoted Ref. FMainardi15 . Namely, for $0<\alpha<1$ and $0<\beta\leq
1/\alpha$, we get two solutions

$g_{H\\!N}(\alpha,\beta;\xi)=\pm\frac{1}{\pi\xi}\frac{\sin\left[\beta\,\theta_{\alpha}(\xi)\right]}{[\xi^{2\alpha}+2\xi^{\alpha}\cos(\pi\alpha)+1]^{\beta/2}}$
(V.15)

where $\theta_{\alpha}(y)$ is defined in Eq. (V.2) and the sign in Eq. (V.15)
depends on the choice of the branch of the arctan function in Eq. (V.2) having
the essential singularity for $\xi=[\cos(\pi\alpha+\pi)]^{-1/\alpha}$.
Equation (V.15) for $\beta=1$ is identically equal to Eq. (V.12). In addition,
Eq. (V.8) implies the non-negativity of Eq. (V.15). The denominator in Eq.
(V.15) is always positive so $g_{H\\!N}(\alpha,\beta;\xi)$ remains non-
negative if $\beta\theta_{\alpha}(1)\in[0,\pi]$ mod $2n\pi$. Since Eq. (V.2)
gives $\theta_{\alpha}(1)\in[0,\pi\alpha]$, this leads to $\beta\leq
1/\alpha$. This fully agrees with the observations made in (RGarrappa16, , Eq.
(3.47)) and (FMainardi15, , Eq. (8)). However, $0<\alpha,\beta<1$ so that
$\beta\leq 1/\alpha$ is always satisfied for these range of parameters
$\alpha$ and $\beta$.

Fig. 2a illustrates $g_{H\\!N}(\alpha,\beta;\xi)$ as the function of $\xi$ for
fixed $\alpha$ and different values of $\beta$. The opposite situation is
presented in Fig. 2b where $\beta$ is fixed while $\alpha$ is changing. To get
both plots we can use either Eq. (V.15) which contains the more friendly
elementary functions, or the Meijer G-representation of
$g_{H\\!N}(\alpha,\beta;\xi)$ given by Eq. (V.10). In the first case we should
separate the range of $\xi$ into two sectors: $\xi\in[0,\xi_{p})$ and
$\xi\in[\xi_{p},\infty)$. At the point $\xi_{p}$ $g_{H\\!N}(\alpha,\beta;\xi)$
changes the sign such that we take $g_{H\\!N}(\alpha,\beta;\xi)$ with a plus
sign in the range $\xi\in[0,\xi_{p})$ and $g_{H\\!N}(\alpha,\beta;\xi)$ with a
minus sign for $\xi\in[\xi_{p},\infty)$. In the Meijer G-case, the required
separation is done automatically via the use of CAS.

Figure 2: (Color online) Plot of $g_{H\\!N}(\alpha,\beta;\xi)$ as a function
of $\xi$ given by Eq. (V.10) for $\alpha,\beta={\rm const}$. In Fig. 2a we
take $\alpha=3/4$ and $\beta=1/3$ (the red continuous curve), $\beta=4/3$ (the
blue dashed curve), and $\beta=7/3$ (the green dashed-dotted curve). It should
be observed that $g_{H\\!N}(\alpha,\beta;\xi)$ contains a negative part for
$\beta>1/\alpha$. In Fig. 2b we have $\beta=1/3$ and $\alpha=1/4$ (the red
continuous curve), $\alpha=1/2$ (the blue dashed curve), and $\beta=3/4$ (the
green dashed-dotted curve).

###  V.4 The relaxation function $n_{H\\!N}(\alpha,\beta;t)$

The HN relaxation function may be obtained in at least in two ways: (i) by
inserting $g_{H\\!N}(\alpha,\beta;t)$ into Eq. (II.5) or (ii) by calculating
the inverse Laplace transform of $\widehat{n}_{H\\!N}(\alpha,\beta;z)$ where
$\widehat{n}_{H\\!N}(\alpha,\beta;z)$ is derived from Eq. (I.1) adjusted for
the HN model. Both these ways lead to the same results, so we choose (ii) as
more convenient for us. Here, $n_{H\\!N}(\alpha,\beta;t)$ reads

$\displaystyle\begin{split}n_{H\\!N}(\alpha,\beta;t)&=\mathscr{L}^{-1}[\widehat{n}_{H\\!N}(\alpha,\beta;z);t]=\\\
&=\mathscr{L}^{-1}[z^{-1};t]-\tau^{-\alpha\beta}\mathscr{L}^{-1}\left[\frac{z^{-1}}{(\tau^{-\alpha}+z^{\alpha})^{\beta}};t\right]\\\
&=1-\big{(}\mbox{$\frac{t}{\tau}$}\big{)}^{\alpha\beta}E_{\alpha,1+\alpha\beta}^{\beta}\big{(}\\!-(\mbox{$\frac{t}{\tau}$})^{\alpha}\big{)},\end{split}$
(V.16)

where $z=\operatorname{i}\\!\omega$ and we used the Laplace transform given by
Eq. (IV.21). To express $n_{H\\!N}(\alpha,\beta;t)$ in the language of the
Meijer G-function we use Eq. (IV.20), which for $\alpha=l/k$ allows us to
write

$\displaystyle\begin{split}n_{H\\!N}(\alpha,\beta;t)&=1-(2\pi)^{\frac{1+l}{2}-k}\,\frac{k^{\beta}l^{-1/2}}{\Gamma(\beta)}\\\
&\times
G_{l+k,k}^{k,k}\left(\frac{l^{l}\tau^{l}}{t^{l}}\Big{|}{\Delta(k,1-\beta),\Delta(l,1)\atop\Delta(k,0)}\right).\end{split}$
(V.17)

The relaxation function $n_{H\\!N}(\alpha,\beta;t)$ given by Eq. (V.17) is
plotted in Fig. 3 for $\alpha=1/2$ and $\beta=1/4,1/2$, and $3/4$.

Figure 3: (Color online) Plot of $n_{H\\!N}(1/2,\beta;t)$ as the function of
$t$ for $\tau=1$ and $\beta=1/2$ (the red continuous curve no. I), $\beta=1$
(the blue dashed curve no. II), and $\beta=2$ (the green dash-dotted curve no.
III).

Using appropriate representation of the Mittag-Leffler function given by Eq.
(IV.20) we can express $n_{H\\!N}(\alpha,\beta;t)$ in terms of these
functions. So, for $\alpha\in(0,1)$ and $\beta=1$, Eq. (V.16) leads to

$\displaystyle\begin{split}n_{H\\!N}(\alpha,1;t)\equiv
n_{CC}(\alpha;t)&=1-E^{1}_{\alpha,1+\alpha}\big{(}\\!-\\!(\mbox{$\frac{t}{\tau}$})^{\alpha}\big{)}\\\
&=E_{\alpha}\big{(}\\!-\\!(\mbox{$\frac{t}{\tau}$})^{\alpha}\big{)}.\end{split}$
(V.18)

For $\alpha=1$ and $\beta\in(0,1)$ it gives

$\displaystyle\begin{split}n_{H\\!N}(1,\beta;t)\equiv
n_{CD}(\beta;t)&=1-E^{\beta}_{1,1+\beta}\big{(}\\!-\\!\mbox{$\frac{t}{\tau}$}\big{)}\\\
&=\Gamma\big{(}\beta,\mbox{$\frac{t}{\tau}$}\big{)}/\Gamma(\beta),\end{split}$
(V.19)

where
$\Gamma(a,y)=\int_{y}^{\infty}u^{\alpha-1}\operatorname{e}^{-u}\operatorname{d\\!}u$
is the upper incomplete gamma function NIST .

Applying Lemma 1 to Eq. (V.16) we find the asymptotics of
$n_{H\\!N}(\alpha,\beta;t)$ which is proportional to

$n_{H\\!N}(\alpha,\beta;t)\propto
1-\frac{1}{\Gamma(1+\alpha\beta)}\left(\frac{t}{\tau}\right)^{\alpha\beta}\quad\text{for}\quad\frac{t}{\tau}\ll
1$

and

$n_{H\\!N}(\alpha,\beta;t)\propto\frac{\beta}{\Gamma(1-\alpha)}\left(\frac{t}{\tau}\right)^{-\alpha}\quad\text{for}\quad\frac{t}{\tau}\gg
1,$

as is given in RGarrappa16 ; KGorska21c ; RHilfer02a .

###  V.5 Evolution equation for $n_{H\\!N}(\alpha,\beta;t)$

According to the general formalism developed in Subsect. III.1, the relaxation
function $n_{H\\!N}(\alpha,\beta;t)$ satisfies the integro-differential
equations (III.4) and (III.7). They contain the memory functions
$M_{H\\!N}(\alpha,\beta;t)$ and $k_{H\\!N}(\alpha,\beta;t)$ connected by the
Sonine equation (III.5). The memory functions $M_{H\\!N}$ and $k_{H\\!N}$ in
the Laplace domain can be determined algebraically using Eqs. (III.9) and
(III.6):

$\widehat{M}_{H\\!N}(\alpha,\beta;z)={\cal
B}^{-1}\\{[1+(z\tau)^{\alpha}]^{\beta}-1\\}^{-1}\quad\text{and}\\\
\widehat{k}_{H\\!N}(\alpha,\beta;z)=\frac{{\cal
B}}{z}\\{[1+(z\tau)^{\alpha}]^{\beta}-1\\}.$ (V.20)

The Laplace inversion with respect to time $t$ domain gives

$M_{H\\!N}(\alpha,\beta;t)=({\cal B}t)^{-1}\sum_{r\geq
0}\big{(}\mbox{$\frac{t}{\tau}$}\big{)}^{\alpha\beta(1+r)}E_{\alpha,\alpha\beta(1+r)}^{\beta(1+r)}\big{[}\\!-\big{(}\mbox{$\frac{t}{\tau}$}\big{)}^{\alpha}\big{]}$

and

$k_{H\\!N}(\alpha,\beta;t)={\cal
B}\big{(}\mbox{$\frac{\tau}{t}$}\big{)}^{\alpha\beta}E_{\alpha,1-\alpha\beta}^{-\beta}\big{[}\\!-\big{(}\mbox{$\frac{t}{\tau}$}\big{)}^{\alpha}\big{]}-{\cal
B},$

and for details see Refs. KGorska21CNSNCa ; KGorska21a ; AAKhamzin14 . As is
shown in Subsect. III.1 the solutions of Eqs. (III.4) and/or (III.7) are
equivalent. Thus, we can choose the equation more convenient for us. Equation
(III.7) contains less complicated memory kernel $k_{H\\!N}(\alpha,\beta;t)$
and immediately leads to the evolution equation proposed either in
(RGarrappa16, , Eq. (3.40)) or in (KGorska21CNSNCa, , Eq. (33)):

${{}^{C\\!}(D^{\alpha}+\tau^{-\alpha})^{\beta}}n_{H\\!N}(\alpha,\beta;t)=-\tau^{-\alpha\beta}.$
(V.21)

The pseudo-operator on the LHS above belongs to the class of Prabhakar-like
integral operators, and for the case under consideration is defined as
(RGarrappa16, , Eq. (B.23))

${{}^{C\\!}(D^{\alpha}+\tau^{-\alpha})^{\beta}}n_{H\\!N}(\alpha,\beta;t)\\\
=\int_{0}^{t}(t-\xi)^{-\alpha\beta}E_{\alpha,1-\alpha\beta}^{-\beta}\big{[}\\!-\\!\big{(}\mbox{$\frac{t-\xi}{\tau}$}\big{)}^{\alpha}\big{]}\dot{n}_{H\\!N}(\alpha,\beta;\xi)\operatorname{d\\!}\xi,$
(V.22)

where we keep the notation of (RGarrappa16, , Appendix B) and denote
$\dot{n}_{H\\!N}(\alpha,\beta;t)=\operatorname{d\\!}{n}_{H\\!N}(\alpha,\beta;t)/\\!\operatorname{d\\!}t$.
It has to be pointed out that the pseudo–operator (V.22) must be distinguished
from the pseudo–operator ${(^{\rm c\\!}D^{\alpha}+\tau^{-\alpha})^{\beta}}$
considered in RGarrappa16  as an alternative to the LHS of Eq. (V.22), with
the Caputo derivative ${{}^{\rm c\\!}D}^{\alpha}$ (see IPodlubny99 ) sitting
in. The difference is evidently seen if we take $\beta=1$, which simplifies HN
relaxation to the CC model. Then, it can be concluded that
${{}^{C\\!}(D^{\alpha}+\tau^{-\alpha})}n_{CC}(\alpha;t)=({{}^{\rm
c\\!}D}^{\alpha}+\tau^{-\alpha})n_{CC}(\alpha;t)-\tau^{-\alpha}$, where
$n_{CC}(\alpha;t)=n_{H\\!N}(\alpha,1;t)$.

###  V.6 Two variants of the subordination approach to
$n_{H\\!N}(\alpha,\beta;t)$

Let us start with the first type of subordination, namely Eq. (III.12),
according to which we have

$n_{H\\!N}(\alpha,\beta;t)=\int_{0}^{\infty}n_{D}(\xi)f_{H\\!N}(\alpha,\beta;\xi,t)\operatorname{d\\!}\xi,$
(V.23)

where

$f_{H\\!N}(\alpha,\beta;\xi,t)={\cal B}\operatorname{e}^{{\cal
B}\xi}\mathscr{L}^{-1}\Big{\\{}\Big{[}\frac{\tau^{\alpha\beta}}{z}(\tau^{-\alpha}+z^{\alpha})^{\beta}-\frac{1}{z}\Big{]}\\\
\times\operatorname{e}^{-\xi{\cal
B}\tau^{\alpha\beta}(\tau^{-\alpha}+z^{\alpha})^{\beta}};t\Big{\\}}.$ (V.24)

represents the probability density distribution of (random) internal time
$\xi$ with respect to the laboratory time $t$. It is found by using Eq.
(III.11). Inserting of Eq. (V.24) into (V.23) results in

$n_{H\\!N}(\alpha,\beta;t)={\cal
B}\int_{0}^{\infty}\\!\mathscr{L}^{-1}\Big{\\{}\Big{[}\frac{\tau^{\alpha\beta}}{z}(\tau^{-\alpha}+z^{\alpha})^{\beta}-\frac{1}{z}\Big{]}\\\
\times\operatorname{e}^{-\xi{\cal
B}\tau^{\alpha\beta}(\tau^{-\alpha}+z^{\alpha})^{\beta}};t\Big{\\}}\operatorname{d\\!}\xi,$
(V.25)

where the term $\exp({\cal B}\xi)$ in Eq. (V.24) canceled
$n_{D}(\xi)=\exp(-{\cal B}\xi)$.

The inverse Laplace transform $\mathscr{L}^{-1}[-,\cdot]$ in the integrand of
(V.25) can be calculated by employing the Efros theorem give the Eq. (IV.5) of
it with $\widehat{G}_{2}(z)=z^{-1}$ and
$\widehat{q}_{2}(z)=\tau^{-\alpha}+z^{\alpha}$. That gives

$\displaystyle\mathscr{L}^{-1}\Big{\\{}\frac{1}{z}\Big{[}\tau^{\alpha\beta}$
$\displaystyle(\tau^{-\alpha}+z^{\alpha})^{\beta}-1\Big{]}\operatorname{e}^{-\xi\mathcal{B}\tau^{\alpha\beta}(\tau^{-\alpha}+z^{\alpha})^{\beta}};t\Big{\\}}\operatorname{d\\!}\xi$
$\displaystyle=\int_{0}^{\infty}\mathscr{L}^{-1}[(\tau^{\alpha\beta}z^{\beta}-1)\operatorname{e}^{-\xi\mathcal{B}\tau^{\alpha\beta}z^{\beta}};u]$
$\displaystyle\times\mathscr{L}^{-1}[z^{-1}\operatorname{e}^{-u(z^{\alpha}+\tau^{-\alpha})};t]\operatorname{d\\!}u,$
(V.26)

where $z\div t$ and $z\div u$ constitute the Laplace pairs for the first and
second inverse Laplace transforms, respectively. The only term which depends
on $\xi$ is equal to $\exp(-\xi\mathcal{B}\tau^{\alpha\beta}z^{\beta})$. Thus,
inserting Eq. (V.6) into Eq. (V.25) and changing the order of integration we
get

$\displaystyle n_{H\\!N}(\alpha,\beta;\tau,t)$
$\displaystyle=\mathcal{B}\int_{0}^{\infty}\\!\mathscr{L}^{-1}\Big{[}(\tau^{\alpha\beta}z^{\beta}-1)\int_{0}^{\infty}\\!\operatorname{e}^{-\xi\mathcal{B}\tau^{\alpha\beta}z^{\beta}}\operatorname{d\\!}\xi;u\Big{]}$
$\displaystyle\times\mathscr{L}^{-1}[z^{-1}\operatorname{e}^{-u(z^{\alpha}+\tau^{-\alpha})};t]\operatorname{d\\!}u$
$\displaystyle=\int_{0}^{\infty}\\!\operatorname{e}^{-u\tau^{-\alpha}}\mathscr{L}^{-1}\Big{[}1-\frac{1}{\tau^{\alpha\beta}z^{\beta}};u\Big{]}$
$\displaystyle\times\mathscr{L}^{-1}[z^{-1}\operatorname{e}^{-uz^{\alpha}};t]\operatorname{d\\!}u.$
(V.27)

Because $\exp(-u\tau^{-\alpha})$ does not depend on $z$ we can extract this
term from the second inverse Laplace transform in Eq. (V.6) and, next, merge
it to the first inverse Laplace transform in Eq. (V.6). Using the property
$\operatorname{e}^{-at}\mathscr{L}^{-1}[\widehat{g}(z);t]=\mathscr{L}^{-1}[\widehat{g}(z+a);t]$
we get

$\displaystyle\operatorname{e}^{-u\tau^{-\alpha}}\\!\mathscr{L}^{-1}\Big{[}1-\frac{1}{\tau^{\alpha\beta}z^{\beta}};u\Big{]}=\mathscr{L}^{-1}\Big{[}1-\frac{1}{(1+\tau^{\alpha}z)^{\beta}};u\Big{]}$
$\displaystyle=\delta(u)-\phi_{CD}(\beta;u).$

The real function $\phi_{C\\!D}(\beta;u)\equiv\phi_{H\\!N}(1,\beta;u)$ is the
response function for the CD model which determines
$\phi_{C\\!D}(\beta;u)=\mathscr{L}^{-1}[(1+Ts)^{-\beta};u]$ with
$T=\tau^{\alpha}$. Having this in mind we present Eq. (V.6) as

$\displaystyle n_{H\\!N}($
$\displaystyle\alpha,\beta;t)=\mathscr{L}^{-1}[z^{-1}\int_{0}^{\infty}\delta(u)\operatorname{e}^{-uz^{\alpha}}\operatorname{d\\!}u;t]$
$\displaystyle-\int_{0}^{\infty}\phi_{C\\!D}(\beta;u)\mathscr{L}^{-1}[z^{-1}\operatorname{e}^{-uz^{\alpha}};t]\operatorname{d\\!}u$
$\displaystyle=1-\int_{0}^{\infty}\\!\\!\phi_{C\\!D}(\beta;u)\mathscr{L}^{-1}[z^{-1}\operatorname{e}^{-uz^{\alpha}};t]\operatorname{d\\!}u$
$\displaystyle=1+\int_{0}^{\infty}\\!\\!\dot{n}_{C\\!D}(\beta;u)\mathscr{L}^{-1}[z^{-1}\operatorname{e}^{-uz^{\alpha}};t]\operatorname{d\\!}u$
(V.28)

where we used taken with minus sign the standard relation
$\phi_{C\\!D}(\beta;u)=-\dot{n}_{C\\!D}(\beta;u)$ between the response
function and the time derivative of relaxation function. Then, we employ the
Leibniz formula which allowed us to shift the derivative over $u$ from the
relaxation function $n_{C\\!D}(\beta;u)$ onto
$\mathscr{L}^{-1}[z^{-1}\exp(-uz^{\alpha});t]$. Thus Eq. (V.6) becomes

$\displaystyle n_{H\\!N}(\alpha,\beta;t)$
$\displaystyle=-\int_{0}^{\infty}\\!\\!n_{C\\!D}(\beta;u)\mathscr{L}^{-1}\Big{[}z^{-1}\frac{\operatorname{d\\!}}{\operatorname{d\\!}u}\operatorname{e}^{-uz^{\alpha}};t\Big{]}\operatorname{d\\!}u$
$\displaystyle=\int_{0}^{\infty}\\!\\!n_{C\\!D}(\beta;u)f(\alpha;u,t)\operatorname{d\\!}u,$
(V.29)

where $f(\alpha;u,t)$ can be expressed in terms of the one-sided Lévy stable
distribution $\varPhi_{\alpha}(\sigma)$ with $\alpha\in(0,1)$ and
$\sigma\in\mathbb{R}_{+}$ KAPenson10

$\displaystyle\begin{split}f(\alpha;u,t)&=\mathscr{L}^{-1}[z^{\alpha-1}\operatorname{e}^{-uz^{\alpha}};t]\\\
&=\frac{t}{\alpha u^{1+1/\alpha}}\varPhi_{\alpha}(tu^{-1/\alpha}).\end{split}$
(V.30)

That means that we arrived at the second type of subordination, in which the
CD relaxation appears, which plays the role of the parent process and is
subordinated by the Lévy-like process described by the PDF $f(\alpha;u,t)$.
This result explicitly confirms our previous claim KGorska21  that the
subordination description of the same process may be realized in different
ways which eventually may be understood as an effect of “nested” processes.

##  VI The Jurlewicz-Weron-Stanislavski model and its physical content

The next illustration which shows of applicability of the so far developed
universal scheme presented in Secs. I-III is the Jurlewicz-Weron-Stanislavsky
(JWS) relaxation model. The JWS relaxation model (for its brief description
see RGarrappa16 ; AStanislavsky17 ) complements and modifies the HN model
leaving unchanged its general structure. The model was introduced in
AJurlewicz10 ; AStanislavsky10 ; KWeron10  to explain discrepancies emerging
in some experiments between the results coming out from the Jonscher URL and
those described by the HN pattern, see (AStanislavsky17, , Fig. 1). (According
to statistical enumeration measuring applicability of different models, the
JWS pattern fits and reproduces the data for approximately 20$\%$ of
relaxation experiments Jonscher83 ; KStanislavski16 .) Therefore it is
legitimate to analyse this model applying theoretical framework and methods
developed in the current work.

###  VI.1 The spectral function
$\widehat{\phi}_{J\\!W\\!S}(\alpha,\beta;\operatorname{i}\\!\omega)$

We shall begin by recalling Eq. (I.4):

$\displaystyle\begin{split}\widehat{\phi}_{J\\!W\\!S}(\alpha,\beta;\operatorname{i}\\!\omega)&=\left[\frac{\widehat{\varepsilon^{\star}}(\omega)-\varepsilon_{\infty}}{\varepsilon-\varepsilon_{\infty}}\right]_{J\\!W\\!S}\\\
&=1-(\operatorname{i}\\!\omega\tau)^{\alpha\beta}\widehat{\phi}_{H\\!N}(\alpha,\beta;\operatorname{i}\\!\omega)\\\
&=1-\widehat{\phi}_{H\\!N}(-\alpha,\beta;\operatorname{i}\\!\omega)\end{split}$
(VI.1)

for $\alpha,\beta\in(0,1]$. For $\alpha=\beta=1$ it approaches the Debye model
whereas for $\alpha\in(0,1)$ and $\beta=1$ tends to the CC pattern. The
relaxation model with $\alpha=1$ and $\beta\in(0,1)$ is called the mirror CD
relaxation (MCD) (AStanislavsky17, , Fig. 4). Note that
$\widehat{\phi}_{J\\!W\\!S}(\alpha,\beta;\operatorname{i}\\!\omega)$ tends to
$1$ for $\omega\to 0$ and it vanishes in the limit $\omega\rightarrow\infty$.
The real and imaginary parts of
$\widehat{\varepsilon^{\star}}_{J\\!W\\!S}(\omega)$ read

$\widehat{\varepsilon^{\prime}}_{J\\!W\\!S}(\alpha,\beta;\omega)=\varepsilon-\frac{(\varepsilon-\varepsilon_{\infty})(\omega\tau)^{\alpha\beta}\cos\big{[}\beta\,\theta_{\alpha/2}\big{(}(\omega
t)^{-2}\big{)}\big{]}}{\big{[}1+2(\omega\tau)^{\alpha}\cos(\pi\alpha/2)+(\omega\tau)^{2\alpha}\big{]}^{\beta/2}}$

and

$\widehat{\varepsilon^{\prime\prime}}_{J\\!W\\!S}(\omega)=\frac{(\varepsilon-\varepsilon_{\infty})(\omega\tau)^{\alpha\beta}\sin\big{[}\beta\,\theta_{\alpha/2}\big{(}(\omega
t)^{-2}\big{)}\big{]}}{\big{[}1+2(\omega\tau)^{\alpha}\cos(\pi\alpha/2)+(\omega\tau)^{2\alpha}\big{]}^{\beta/2}},$

respectively, with $\theta_{\alpha}(y)$ introduced in Eq. (V.2).

###  VI.2 The response function $\phi_{J\\!W\\!S}(\alpha,\beta;t)$

The response function $\phi_{J\\!W\\!S}(\alpha,\beta;t)$ is the inverse
Laplace transform of $\widehat{\phi}_{J\\!W\\!S}(\alpha,\beta;z)$,

$\displaystyle\phi_{J\\!W\\!S}(\alpha,\beta;t)$
$\displaystyle=\mathscr{L}^{-1}[1;t]-\mathscr{L}^{-1}\Big{[}\frac{z^{\alpha\beta}}{(\tau^{-\alpha}+z^{\alpha})^{\beta}};t\Big{]}$
$\displaystyle=\delta(t)-\tau^{-1}\big{(}\mbox{$\frac{t}{\tau}$}\big{)}^{-1}E_{\alpha,0}^{\beta}\big{(}-\big{(}\mbox{$\frac{t}{\tau}$}\big{)}^{\alpha}\big{)},$
(VI.2)

where we employ the Laplace transform of the Mittag-Leffler function given by
Eq. (IV.21). Notice that the $\delta$-Dirac distribution corresponds to the
integrable singularity at $t=0$. We remark that Eq. (VI.2) is also presented
in (RGarrappa16, , Eq. (3.43)). For $\alpha\in(0,1)$ and $\beta=1$ we obtain
the CC model, see Remark 3, while the MCD pattern appears for $\alpha=1$ and
$\beta\in(0,1)$. From Lemma 3 we can formally rewrite the JWS response
function in terms of the HN response function, namely

$\phi_{J\\!W\\!S}(\alpha,\beta;t)=\delta(t)-\phi_{H\\!N}(-\alpha,\beta;t).$

To rephrase it in terms of functions implemented in the CAS we use Eq. (VI.2)
and the Mittag-Leffler representation of the response function
$\phi_{J\\!W\\!S}(\alpha,\beta;t)$. This enables one to express
$\phi_{J\\!W\\!S}(\alpha,\beta;t)$ for rational $\alpha=l/k$ in terms of
Meijer G-function and the finite sums of generalized hypergeometric functions.
With the help of Eqs. (IV.20) and (V.3) we have

$\displaystyle\phi_{J\\!W\\!S}($ $\displaystyle
l/k,\beta;t)=\delta(t)-(2\pi)^{\frac{1+l}{2}-k}\,\frac{k^{\beta}l^{1/2}}{\tau\Gamma(\beta)}\left(\frac{t}{\tau}\right)^{\\!-1}$
$\displaystyle\times
G^{k,k}_{k+l,k}\left(\frac{l^{l}\tau^{l}}{t^{l}}\Big{|}{\Delta(k,1),\Delta(l,0)\atop\Delta(k,\beta)}\right)$
(VI.3)
$\displaystyle=\delta(t)-\frac{1}{\tau}\sum_{j=0}^{k-1}\frac{(-1)^{j}(\beta)_{j}}{j!\,\Gamma(\frac{l}{k}j)}\Big{(}\frac{t}{\tau}\Big{)}^{\frac{l}{k}j-1}$
(VI.4)
$\displaystyle\times{{}_{1+k}F_{l+k}}\left({1,\Delta(k,\beta+j)\atop\Delta(k,1+j),\Delta(l,\frac{l}{k}j)};\frac{(-1)^{k}t^{l}}{\tau^{l}l^{l}}\right)\,.$

The asymptotic behavior of the Mittag-Leffler functions presented by Lemma 1
result in (RGarrappa16, , Eq. (3.45)), that is

$\phi_{J\\!W\\!S}(\alpha,\beta;t)\propto\frac{\beta}{\tau\Gamma(\alpha)}\left(\frac{t}{\tau}\right)^{\alpha-1}\quad\text{for}\quad\frac{t}{\tau}\ll
1\quad\text{and}\\\
\phi_{J\\!W\\!S}(\alpha,\beta;t)\propto-\frac{1}{\tau\Gamma(-\alpha\beta)}\left(\frac{t}{\tau}\right)^{-\alpha\beta-1}\quad\text{for}\quad\frac{t}{\tau}\gg
1.$ (VI.5)

Because $\widehat{\phi}_{J\\!W\\!S}(\alpha,\beta;s)$, $s>0$, has the CM
character for $\beta\leq 1/\alpha$, then thanks to the Bernstein theorem
(Sect. IV.1) we know that $\phi_{J\\!W\\!S}(\alpha,\beta;t)$ is a non-negative
function in the same range of parameters $\alpha$ and $\beta$. That is
confirmed in Fig. 4.

Figure 4: (Color online) Plot of $\phi_{J\\!W\\!S}(1/2,\beta;t)$ as a function
of $t$ given by Eq. (VI.2) for $\tau=1$ and $\beta=1/2$ (the red continuous
curve no. I), $\beta=1$ (the blue dashed curve no. II), $\beta=2$ (the black
dash-dotted curve no. III), and $\beta=3$ (the green dotted curve no. IV).

###  VI.3 The probability density $g_{J\\!W\\!S}(\alpha,\beta;\xi)$

By analogy with to Subsect. V.3 we can represent the spectral function
$\widehat{\phi}_{J\\!W\\!S}(\alpha,\beta;z)$, $z=\operatorname{i}\\!\omega$,
as follows

$\widehat{\phi}_{J\\!W\\!S}(\alpha,\beta;z)=\int_{0}^{\infty}\frac{p_{J\\!W\\!S}(\alpha,\beta;\xi)}{\xi+z}\operatorname{d\\!}\xi,$
(VI.6)

which for $z=s\in\mathbb{R}_{+}$ is a S function KGorska21c . The Schwinger
parametrization enables us to write
$\widehat{\phi}_{J\\!W\\!S}(\alpha,\beta;z)=\mathscr{L}[\phi_{J\\!W\\!S}(\alpha,\beta;t);z]$,
where from Eq. (II.6) it follows that
$\phi_{J\\!W\\!S}(\alpha,\beta;t)=\mathscr{L}[p_{J\\!W\\!S}(\alpha,\beta;\xi);t]$
with
$p_{J\\!W\\!S}(\alpha,\beta;\xi)=\xi\,g_{J\\!W\\!S}(\alpha,\beta;\xi)/\tau$.
Hence,

$g_{J\\!W\\!S}(\alpha,\beta;\xi)=\frac{\tau}{\xi}\,\mathscr{L}^{-1}[\phi_{J\\!W\\!S}(\alpha,\beta;t);\xi].$
(VI.7)

Derivation of the exact form of $g_{J\\!W\\!S}(\alpha,\beta;\xi)$ for
$\alpha=l/k$ relies on calculating the inverse Laplace transform sitting in
Eq. (VI.7). To do this we employ the Meijer G-representation of the JWS
response function. That gives

$\mathscr{L}^{-1}[\phi_{J\\!W\\!S}(l/k,\beta;t);\xi]=\mathscr{L}^{-1}[\delta(t);\xi]\\\
-(2\pi)^{\frac{1+l}{2}-k}\frac{\sqrt{l}k^{\beta}}{\tau\Gamma(\beta)}\,\mathscr{L}^{-1}\left[\Big{(}\frac{t}{\tau}\Big{)}^{-1}\right.\\\
\left.\times
G^{k,k}_{k+l,k}\left(\frac{l^{l}\tau^{l}}{t^{l}}\Big{|}{\Delta(k,1),\Delta(l,0)\atop\Delta(k,\beta)}\right);\xi\right].$
(VI.8)

The first inverse Laplace transform in Eq. (VI.8),
$\mathscr{L}^{-1}[\delta(t);\xi]$, can be obtained by using the so-called
limits representation of $\delta$-distribution, namely
$\delta(t)=\lim_{\epsilon\to+0}\epsilon\,t^{\epsilon-1}/2$ for
$t\in\mathbb{R}_{+}$ online1 , and next changing the order of the inverse
Laplace transform and the limit. That leads to

$\displaystyle\begin{split}\mathscr{L}^{-1}[\delta(t);\xi]&=\mathscr{L}^{-1}\Big{[}\\!\lim_{\epsilon\to+0}\frac{1}{2}\epsilon\,t^{\epsilon-1};\xi\Big{]}\\\
&=\frac{1}{2}\lim_{\epsilon\to+0}\epsilon\,\mathscr{L}^{-1}[t^{\epsilon-1};\xi]\\\
&=\frac{1}{2}\lim_{\epsilon\to+0}\frac{\epsilon\,\xi^{-\epsilon}}{\Gamma(1-\epsilon)}=0,\qquad\text{for}\quad\xi\in\mathbb{R}_{+}.\end{split}$

To calculate the second inverse Laplace transform we use Eq. (IV.10):

$\displaystyle g_{J\\!W\\!S}($ $\displaystyle
l/k,\beta;\xi)=-(2\pi)^{\frac{1+l}{2}-k}\frac{\sqrt{l}k^{\beta}}{\Gamma(\beta)}\,\frac{\tau}{\xi}$
(VI.9)
$\displaystyle\times\mathscr{L}^{-1}\left[t^{-1}G^{k,k}_{k+l,k}\left(\frac{l^{l}\tau^{l}}{t^{l}}\Big{|}{\Delta(k,1),\Delta(l,0)\atop\Delta(k,\beta)}\right);\xi\right]$
$\displaystyle=-(2\pi)^{l-k}\frac{k^{\beta}}{\Gamma(\beta)}\,\frac{1}{\xi}\,G^{k,k}_{k+l,k+l}\left(\xi^{l}\Big{|}{\Delta(k,1),\Delta(l,0)\atop\Delta(k,\beta),\Delta(l,0)}\right).$

The last formula becomes more readable if the Meijer G-function is expressed
in terms of the finite sum of the hypergeometric functions ${{}_{p}F_{q}}$
according to Eq. (IV.13). That implies

$\displaystyle\begin{split}g_{J\\!W\\!S}(&l/k,\beta;\xi)=\frac{1}{\pi}\sum_{j=0}^{k-1}\frac{(-1)^{j}}{j!}\frac{(\beta)_{j}}{\xi^{1-\frac{l}{k}n_{j}}}\sin(\mbox{$\frac{l}{k}$}n_{j}\pi)\\\
&\times{{}_{k+1}F_{k}}\left({1,\Delta(k,n_{j})\atop\Delta(k,1+j)};(-1)^{l-k}\xi^{l}\right),\end{split}$
(VI.10)

with $n_{j}=\beta+j$. This is similar to $g_{H\\!N}(l/k,\beta;\xi)$ given by
Eq. (V.11) except of terms involving dependence on $\xi$. For the CC model
with $\beta=1$ and $\alpha\in(0,1)$ we obtain
$g_{J\\!W\\!S}(\alpha,1;\xi)\equiv g_{CC}(\alpha;\xi)$ given by Eq. (V.12).
For the MCD model, i.e., $\alpha=l/k=1$, Eq. (VI.10) becomes

$g_{J\\!W\\!S}(1,\beta;\xi)\equiv
g_{M\\!C\\!D}(\beta;\xi)=\frac{\sin(\beta\pi)}{\pi\xi(\xi^{-1}-1)^{\beta}}\Theta(1-\xi),$

in which the Heaviside step function $\Theta(\cdot)$ ensures that
$g_{M\\!C\\!D}(\beta;\xi)$ is a real function for $0<\xi<1$.

Following the same procedure as in the case of the HN model we insert in Eq.
(VI.10) the series definition of the generalized hypergeometric function
${{}_{k+1}F_{k}}$, i.e., Eq. (IV.6). That gives the double sum: one is over
$j=0,\ldots,k-1$ and another $r=0,1,\ldots$. They can be reduced to the single
sum due to relation $\sum_{r\geq 0}\sum_{j=0}^{k-1}a_{j+kr}=\sum_{r\geq
0}a_{r}$. Thus, we get an analogue of Eq. (V.13)

$\displaystyle g_{J\\!W\\!S}($
$\displaystyle\alpha,\beta;\xi)=\frac{1}{\pi}\,\sum_{r=0}^{\infty}\frac{(-1)^{r}}{r!}\frac{(\beta)_{r}}{\xi^{1-\alpha(\beta+r)}}\sin[\alpha(\beta+r)\pi]$
$\displaystyle=\frac{1}{\pi\xi}\,\operatorname{\mathfrak{Im}}\left\\{\sum_{r=0}^{\infty}\frac{(-1)^{r}}{r!}(\beta)_{r}\,(\xi\operatorname{e}^{\operatorname{i}\\!\pi})^{\alpha(\beta+r)}\right\\},$
(VI.11)

where $(\xi\operatorname{e}^{\operatorname{i}\\!\pi})^{\alpha(\beta+r)}$ can
be defined by the integral representation of $\Gamma$ function, i.e
$(\xi\operatorname{e}^{\operatorname{i}\\!\pi})^{\lambda}=-\int_{0}^{\infty}\operatorname{e}^{-\xi
u}(u\operatorname{e}^{\operatorname{i}\\!\pi})^{\lambda-1}\operatorname{d\\!}u/\Gamma(-\lambda)$.
Thereafter, we substitute it into Eq. (VI.11) and change the order of sum and
integral. Employing the series-form definition of the three-parameter Mittag-
Leffler function (IV.14) leads to

$\displaystyle g_{J\\!W\\!S}($
$\displaystyle\alpha,\beta;\xi)=\frac{1}{\pi\xi}\operatorname{\mathfrak{Im}}\left\\{\operatorname{e}^{-\operatorname{i}\alpha\beta\pi}\int_{0}^{\infty}\operatorname{e}^{-\xi
u}u^{-\alpha\beta-1}\right.$ $\displaystyle\left.\times
E_{-\alpha,-\alpha\beta}^{\beta}(-u^{-\alpha}\operatorname{e}^{-i\alpha\pi})\operatorname{d\\!}u\right\\}$
$\displaystyle=\frac{1}{\pi\xi}\operatorname{\mathfrak{Im}}\left\\{\int_{0}^{\infty}\operatorname{e}^{-\xi
u}u^{-1}E_{\alpha,0}^{\beta}(-u^{\alpha}\operatorname{e}^{\operatorname{i}\\!\alpha\pi})\operatorname{d\\!}u\right\\},$
(VI.12)

where the relation between two Prabhakar functions comes from the Lemma 3.
Inserting the Laplace transform Eq. (IV.21) into Eq. (VI.3) and employing de
Moivre’s formula to calculate the imaginary part we have

$g_{J\\!W\\!S}(\alpha,\beta;\xi)=\pm\frac{1}{\pi}\,\frac{\xi^{\alpha\beta-1}\sin[\beta\,\theta_{\alpha}(\xi^{-1})]}{[\xi^{2\alpha}+2\xi^{\alpha}\cos(\pi\alpha)+1]^{\beta/2}},$
(VI.13)

where $\theta_{\alpha}(y)$ is defined in Eq. (V.2). Repeating the
considerations quoted below Eq. (V.15) we can deduce that
$g_{J\\!W\\!S}(\alpha,\beta;\xi)$ is non-negative only for $0<\alpha\leq 1$
and $0<\beta\leq 1/\alpha$ \- that is seen in Figs. 5. In order to avoid
dividing the interval $\xi$ into two intervals like in
$g_{H\\!N}(\alpha,\beta;\xi)$ we illustrate $g_{J\\!W\\!S}(\alpha,\beta;\xi)$
by employing its Meijer G-representation given by Eq. (VI.9). We have two
plots, Fig. 5a and Fig. 5b, where $g_{J\\!W\\!S}$ are presented as the
function of $\xi$ for given $\alpha$ and $\beta$. In Fig. 5a we keep the value
of $\alpha={\rm const}$ and change $\beta$ whereas in Fig. 5b we consider
opposite situation.

Figure 5: (Color online) Plot of $g_{J\\!W\\!S}(\alpha,\beta;\xi)$ given by
Eq. (VI.9) for given parameters $\alpha$ and $\beta$ as the function of $\xi$.
Descriptions are the same as in Figs. 2.

###  VI.4 The relaxation function $n_{J\\!W\\!S}(\alpha,\beta;t)$

Analogically to the HN case the relaxation function for the JWS model can be
obtained in least in two ways: either inserting
$g_{J\\!W\\!S}(\alpha,\beta;\xi)$ found in Subsect. VI.3 into Eq. (II.5) or
using Eq. (I.1) which connects the Laplace form of the relaxation function
$\widehat{n}_{J\\!W\\!S}(\alpha,\beta;\operatorname{i}\\!\omega)$ with the
spectral function
$\widehat{\phi}_{J\\!W\\!S}(\alpha,\beta;\operatorname{i}\\!\omega)$.
Following the second approach we get (RGarrappa16, , Eq. (3.44)):

$n_{J\\!W\\!S}(\alpha,\beta;t)=E_{\alpha,1}^{\beta}\big{(}-(t/\tau)^{\alpha}\big{)},$
(VI.14)

which after using Lemma 4 and Eq. (I.2) can be rewritten as

$\displaystyle\begin{split}n_{J\\!W\\!S}(\alpha,\beta;t)&=\int_{0}^{t}\phi_{H\\!N}(-\alpha,\beta;u)\operatorname{d\\!}u\\\
&=1-n_{H\\!N}(-\alpha,\beta;t).\end{split}$ (VI.15)

Moreover, with the help of Eqs. (IV.20) we can express Eq. (VI.14) for
$\alpha=l/k$ in the language of the Meijer G-representation:

$\displaystyle n_{J\\!W\\!S}($
$\displaystyle\alpha,\beta;t)=(2\pi)^{\frac{1+l}{2}-k}\,\frac{k^{\beta}l^{-1/2}}{\Gamma(\beta)}$
$\displaystyle\times
G^{k,k}_{k+l,k}\left(\frac{l^{l}\tau^{l}}{t^{l}}\Big{|}{\Delta(k,1),\Delta(l,1)\atop\Delta(k,\beta)}\right)$
(VI.16)

plotted for $\tau=1$, $\alpha=1/2$ and $\beta=1/4,1/2,3/4$ in Fig. 6.

Figure 6: (Color online) Plot of $n_{J\\!W\\!S}(1/2,\beta;t)$ given by Eq.
(VI.4). The curves have analogical characterizations as in Fig. 3.

For $\alpha\in(0,1)$ and $\beta=1$ the relaxation function
$n_{J\\!W\\!S}(\alpha,\beta;t)$ goes to the CC relaxation function (V.18)
whereas for $\alpha=1$ and $\beta\in(0,1)$ it tends to

$n_{J\\!W\\!S}(1,\beta;t)\equiv
n_{M\\!C\\!D}(\beta;t)=E_{1,1}^{\beta}(-t/\tau)={{}_{1}F_{1}}\left({\beta\atop
1};-\frac{t}{\tau}\right).$

The asymptotics of $n_{J\\!W\\!S}(\alpha,\beta;t)$ obtained from Eqs. (IV.15)
and (IV.16) reads

$n_{J\\!W\\!S}(\alpha,\beta;t)\propto
1-\frac{\beta}{\Gamma(1+\alpha)}\left(\frac{t}{\tau}\right)^{\alpha}\quad\text{for}\quad\frac{t}{\tau}\ll
1\quad\text{and}\\\
n_{J\\!W\\!S}(\alpha,\beta;t)\propto\frac{1}{\Gamma(1-\alpha\beta)}\left(\frac{t}{\tau}\right)^{-\alpha\beta}\quad\text{for}\quad\frac{t}{\tau}\gg
1.$ (VI.17)

like in (RGarrappa16, , Eq. (3.46)).

###  VI.5 Evolution equation for $n_{J\\!W\\!S}(\alpha,\beta;t)$

The memory kernels $M_{J\\!W\\!S}(\alpha,\beta;t)$ and
$k_{J\\!W\\!S}(\alpha,\beta;t)$ obtained with the help of Eqs. (III.9) and
(III.6) furnish

$M_{J\\!W\\!S}(\alpha,\beta;t)={\cal
B}^{-1}\\{t^{-1}E_{\alpha,0}^{-\beta}[-(t/\tau)^{\alpha}]-\delta(t)\\}\quad\text{and}\\\
k_{J\\!W\\!S}(\alpha,\beta;t)={\cal B}\sum_{r\geq
0}E^{\beta(1+r)}_{\alpha,1}[-(t/\tau)^{\alpha}],$ (VI.18)

see KGorska21CNSNCa ; KGorska21b , and in the Laplace domain they read

$\widehat{M}_{J\\!W\\!S}(\alpha,\beta;z)={\cal
B}^{-1}[z^{-\alpha\beta}(\tau^{-\alpha}+z^{\alpha})^{\beta}-1]$

and

$\widehat{k}_{J\\!W\\!S}(\alpha,\beta;z)={\cal B}\\{z[1+(\tau
z)^{-\alpha}]^{\beta}-z\\}^{-1}.$ (VI.19)

The expression $M_{J\\!W\\!S}(\alpha,\beta;t)$ is simpler so to get the
evolution equation for the JWS model we take Eq. (III.4) instead of Eq.
(III.7). That gives:

$\int_{0}^{t}(t-\xi)^{-1}E_{\alpha,0}^{-\beta}[-\tau^{-\alpha}(t-\xi)^{\alpha}]\,n_{J\\!W\\!S}(\alpha,\beta;\xi)\operatorname{d\\!}\xi=1,$
(VI.20)

where we have used the sampling property of the Dirac $\delta$–function,
$\int_{0}^{t}\delta(t-\xi)n(\xi)\operatorname{d\\!}\xi=n(t)$ for
$t\in\mathbb{R}_{+}$. Proceeding as in Ref. KGorska21CNSNCa  we can state that
Eq. (VI.20) coincides with the equation (RGarrappa16, , Eq. (3.50)):

$(D^{\alpha}+\tau^{-\alpha})^{\beta}n_{J\\!W\\!S}(\alpha,\beta;t)=\frac{t^{-\alpha\beta}}{\Gamma(1-\alpha\beta)}.$

(completed with a suitable initial condition). The pseudo-operator in the
above equation reads

$(D^{\alpha}+\tau^{-\alpha})^{\beta}n_{J\\!W\\!S}(\alpha,\beta;t)\\\
=\frac{\operatorname{d\\!}}{\operatorname{d\\!}t}\int_{0}^{t}(t-\xi)^{-\alpha\beta}E_{\alpha,1-\alpha\beta}^{-\beta}\Big{[}-\frac{(t-\xi)^{\alpha}}{\tau^{\alpha}}\Big{]}\,n_{J\\!W\\!S}(\alpha,\beta;\xi)\,\operatorname{d\\!}\xi.$
(VI.21)

It was introduced in AStanislavsky15 ; KStanislavski16  and discussed in
Appendix B of RGarrappa16 . For $\beta=1$ and $\alpha\in(0,1)$ it gives the
evolution equation written in terms of fractional derivative in the sense of
Riemann-Liouville derivative $D^{\alpha}$, $\alpha\in(0,1)$ IPodlubny99 , used
instead of the fractional derivative in the Caputo sense. Using the link
between these derivatives, i.e.,
$({{}^{c\\!}D^{\alpha}}f)(t)=(D^{\alpha}f)(t)-t^{-\alpha}f(0)/\Gamma(1-\alpha)$,
we get the same formula.

###  VI.6 Two kinds of subordination approaches of
$n_{J\\!W\\!S}(\alpha,\beta;t)$

As the first type of subordination we take the Debye subordination, i.e., Eq.
(III.12) with $f_{\hat{\Psi}}(\xi,t)$ given now by

$f_{J\\!W\\!S}(\alpha,\beta;\xi,t)={\cal B}\operatorname{e}^{{\cal
B}\xi}\mathscr{L}^{-1}\Big{[}\frac{z^{\alpha\beta-1}}{(z^{\alpha}+\tau^{-\alpha})^{\beta}-z^{\alpha\beta}}\\\
\times\operatorname{e}^{-{\cal
B}\xi\frac{(z^{\alpha}+\tau^{-\alpha})^{\beta}}{(z^{\alpha}+\tau^{-\alpha})^{\beta}-z^{\alpha\beta}}};t\Big{]}.$
(VI.22)

After the cancelation of the Debye relaxation function by $\exp({\cal B}\xi)$
coming from Eq. (VI.22) we get

$n_{J\\!W\\!S}(\alpha,\beta;x,t)={\cal
B}\int_{0}^{\infty}\\!\mathscr{L}^{-1}\left[\frac{z^{\alpha\beta-1}}{(z^{\alpha}+\tau^{-\alpha})^{\beta}-z^{\alpha\beta}}\right.\\\
\left.\times\operatorname{e}^{-\xi\frac{{\cal
B}(z^{\alpha}+\tau^{-\alpha})^{\beta}}{(z^{\alpha}+\tau^{-\alpha})^{\beta}-z^{\alpha\beta}}};t\right]\operatorname{d\\!}\xi.$
(VI.23)

To calculate the inverse Laplace transform in the integrand of Eq. (VI.23) we
apply once more the Efros theorem but this time with $\widehat{G}_{2}(z)$ and
$\widehat{q}_{2}(z)$ given below Eq. (V.25). Thus, the RHS of Eq. (VI.23)
becomes

$\displaystyle\mathscr{L}^{-1}$
$\displaystyle\left[\frac{z^{\alpha\beta-1}}{(z^{\alpha}+\tau^{-\alpha})^{\beta}-z^{\alpha\beta}}\operatorname{e}^{-\xi\frac{{\cal
B}(z^{\alpha}+\tau^{-\alpha})^{\beta}}{(z^{\alpha}+\tau^{-\alpha})^{\beta}-z^{\alpha\beta}}};t\right]$
$\displaystyle=\int_{0}^{\infty}\\!\mathscr{L}^{-1}\left[\frac{(z-\tau^{-\alpha})^{\beta}}{z^{\beta}-(z-\tau^{-\alpha})^{\beta}}\operatorname{e}^{-\xi\frac{{\cal
B}z^{\beta}}{z^{\beta}-(z-\tau^{-\alpha})^{\beta}}};u\right]$
$\displaystyle\times\mathscr{L}^{-1}[z^{-1}\operatorname{e}^{-u(z^{\alpha}+\tau^{-\alpha})};t]\operatorname{d\\!}u$

Inserting it into Eq. (VI.23) and changing the order of integration we can
simplify the calculations. It results in

$\displaystyle n_{J\\!W\\!S}($ $\displaystyle\alpha,\beta;x,t)={\cal
B}\int_{0}^{\infty}\\!\operatorname{e}^{-u\tau^{-\alpha}}\mathscr{L}^{-1}\left[\frac{(z-\tau^{-\alpha})^{\beta}}{z^{\beta}-(z-\tau^{-\alpha})^{\beta}}\right.$
$\displaystyle\left.\times\int_{0}^{\infty}\\!\operatorname{e}^{-\xi\frac{{\cal
B}z^{\beta}}{z^{\beta}-(z-\tau^{-1\alpha})}}\operatorname{d\\!}\xi;u\right]\,\mathscr{L}^{-1}[z^{-1}\operatorname{e}^{-uz^{\alpha}};t]\operatorname{d\\!}u$
$\displaystyle=\int_{0}^{\infty}\\!\mathscr{L}^{-1}\left[\frac{(z-\tau^{-\alpha})^{\beta}}{z^{\beta}};u\right]\operatorname{e}^{-u\tau^{-\alpha}}$
$\displaystyle\times\mathscr{L}^{-1}[z^{-1}\operatorname{e}^{-uz^{\alpha}};t]\operatorname{d\\!}u.$
(VI.24)

The next observation is

$\displaystyle\mathscr{L}^{-1}\Big{[}$
$\displaystyle\frac{(z-\tau^{-\alpha})^{\beta}}{z^{\beta}};u\Big{]}\operatorname{e}^{-u\tau^{-\alpha}}=\mathscr{L}^{-1}[(z\tau^{\alpha})^{\beta}\widehat{\phi}_{C\\!D}(\beta;z);u]$
$\displaystyle=\mathscr{L}^{-1}[1-\widehat{\phi}_{M\\!C\\!D}(\beta;z);u],$
(VI.25)

where we applied Eq. (VI.1) for $\alpha=1$. Then, it turns out from that
$\widehat{\phi}_{M\\!C\\!D}(\beta;z)=1-(z\tau^{\alpha})^{\beta}\widehat{\phi}_{C\\!D}(\beta;z)$.
Since $\mathscr{L}^{-1}[1;u]=\delta(u)$ and the inverse Laplace transform of
the spectral function is equal to the response function, then Eq. (VI.6) is
expressed as $\delta(u)-\phi_{M\\!C\\!D}(\beta;u)$. Consequently,
$n_{J\\!W\\!S}$ comes out as

$\displaystyle n_{J\\!W\\!S}($
$\displaystyle\alpha,\beta;x,t)=\mathscr{L}^{-1}\left[z^{-1}\int_{0}^{\infty}\\!\delta(u)\operatorname{e}^{-uz^{\alpha}}\operatorname{d\\!}u;t\right]$
$\displaystyle-\int_{0}^{\infty}\phi_{M\\!C\\!D}(\beta;u)\mathscr{L}^{-1}\left[z^{-1}\operatorname{e}^{-uz^{\alpha}};t\right]\operatorname{d\\!}u.$
(VI.26)

The last steps to complete the calculation are: the use
$\phi_{M\\!C\\!D}(\beta;u)=-\dot{n}_{M\\!C\\!D}(\beta;u)$, the Leibniz
formula, and, finally, rewrite the RHS of Eq. (VI.6) as

$\displaystyle n_{J\\!W\\!S}($
$\displaystyle\alpha,\beta;x,t)=\int_{0}^{\infty}n_{M\\!C\\!D}(\beta;u)\mathscr{L}^{-1}\left[z^{-1}\frac{d}{du}\operatorname{e}^{-uz^{\alpha}};t\right]\operatorname{d\\!}u$
$\displaystyle=\int_{0}^{\infty}n_{M\\!C\\!D}(\beta;u)\mathscr{L}^{-1}\left[z^{\alpha-1}\operatorname{e}^{-uz^{\alpha}};t\right]\operatorname{d\\!}u$
$\displaystyle=\int_{0}^{\infty}n_{M\\!C\\!D}(\beta;u)f(\alpha;u,t)\operatorname{d\\!}u$
(VI.27)

with $f(\alpha;u,t)$ given by Eq. (V.30). Equation (VI.6) can be interpreted
as the second type of subordination in which $n_{M\\!C\\!D}(\beta;u)$ is
subordinated by $f(\alpha;u,t)$.

We see that the HN and the JWS relaxation models lead to at least two types of
subordinations: one described by Eqs. (V.23) and (V.6) for the HN model, as
well as second describe by Eqs. (VI.23) and (VI.6) for the JWS model. Looking
for physical interpretation of this fact enables us to suspect that with the
growing complexity of the relaxing system, the simple partition of the process
into two components, namely the parent and leading processes is not enough to
reflect and understand all properties of the system, in particular it may be
necessary to take into account the fact that the parent and leading processes
can have a complex structure on their own.

##  VII Summary

The broadband dielectric spectroscopy allows us to obtain experimental data
which, if fitted to the spectral functions
$\widehat{\phi}_{(\cdot)}(\operatorname{i}\\!\omega)$, enable us to classify
the latter as corresponding to one among of the standard relaxation models.
Simpler of them, the CC, CD, and MCD models, emerge as reductions of three
parameter ones, namely the HN and JWS patterns. HN and JWS models reduce to
the CC relaxation for $\alpha\in(0,1)$ and $\beta=1$. For $\alpha=1$ and
$\beta\in(0,1)$ the HN model goes to CD pattern, whereas the JWS pattern leads
to the MCD model. All these spectral functions were tabulated in Tab. 4 from
which one sees that

$\widehat{\phi}_{J\\!W\\!S}(\alpha,\beta;\operatorname{i}\\!\omega)+\widehat{\phi}_{H\\!N}(-\alpha,\beta;\operatorname{i}\\!\omega)=1.$
(VII.1)

The knowledge of phenomenologically found spectral functions is crucial for
our investigations. From one side, the methods of dielectric relaxation theory
and the Laplace transform enable us to recover the response $\phi(t)$ and
relaxation $n(t)$ functions, see Eqs. (I.1) and (I.2), from the knowledge of
spectral function. Obviously, the reverse procedure is also justified - we can
transform $\phi(t)$ and/or $n(t)$ into
$\widehat{\phi}(\operatorname{i}\\!\omega)$. For the readers convenience we
itemized $\phi(t)$ and $n(t)$ in Tabs. 4 and 4 which give

$\phi_{J\\!W\\!S}(\alpha,\beta;t)+\phi_{H\\!N}(-\alpha,\beta;t)=\delta(t)$
(VII.2)

and

$n_{J\\!W\\!S}(\alpha,\beta;t)+n_{H\\!N}(-\alpha,\beta;t)=1.$ (VII.3)

From another side the importance of spectral function
$\widehat{\phi}(\operatorname{i}\\!\omega)$ comes from the fact that its
knowledge enables us to find the memory functions $k(t)$ and $M(t)$. The
latter are basic objects used to determine the time evolution of the system
under study, see Eq. (III.9). To develop this idea we took into account that
the memory functions $M(t)$ and $k(t)$ are mutually related by the Sonine
equation. This implied that the relevant evolution equations, see Eqs.
(III.3)/(III.4) and (III.7), led to the same results. Thus we concluded that
the search for principles which govern the evolution may be done in two free
chosen ways. The simplest is to consider integro-differential equations
interpreted as memory dependent. For standard relaxation function these
equations are presented in Tab. 5. Here, we would like to pay the readers
attention that this approach prefers the evolution equations for the HN and
JWS models in the form which involves the pseudo-operators
${{}^{C\\!}(D^{\alpha}+\tau^{-\alpha})^{\beta}}$ and
${(D^{\alpha}+\tau^{-\alpha})^{\beta}}$. These pseudo-operators are defined by
Eq. (V.22) and Eq. (VI.21) and differ only by the position of the time
derivative, like it is in the fractional derivatives of the Caputo and
Riemann-Liouville senses. The relation between them is given though
(RGarrappa16, , Eqs. (B.24) or (B.25)), this is

${{}^{C\\!}(D^{\alpha}+\tau^{-\alpha})^{\beta}}n(t)\\\
={(D^{\alpha}+\tau^{-\alpha})^{\beta}}n(t)-t^{-\alpha\beta}E_{\alpha,1-\alpha\beta}^{-\beta}\big{(}-(\mbox{$\frac{t}{\tau}$})^{\alpha}\big{)}n(0+).$
(VII.4)

The second observation concerning the spectral functions is that they are
involved in the definition of characteristic functions $\widehat{\Psi}(z)$,
named also the Laplace-Lévy exponents. This correspondence goes through using
the memory function $M(t)$ and appears to be essential if one links the memory
dependent evolution schemes and subordination approach, see Eq. (III.11). Our
approach introduces an essential, previously unnoticed, novelty - in the
dielectric relaxation theory we can distinguish, at least two, subordination
patterns. Besides of the subordination which comes from the Schwinger
parametrization (called by us the basic one), see Eq. (III.12), we found its
alternative coming from the Efros theorem. Both these subordinations are
connected to the various choice of “internal” timing, see Tab. 6. For instant,
the HN relaxation function $n_{HN}(t)$ can be build for from the Debye
relaxation in which we introduce internal time $\xi$. The relation between the
physical time $t$ and internal time $\xi$ is given by the PDF of leading
process, here $f_{\widehat{\Psi}(z)}(\xi,t)$. Another possibility to get
$n_{HN}(t)$ is from the CD relaxation function in which the timing
characterized by the PDF $f(\alpha;\xi,t)$.

$\widehat{\phi}_{HN}$ | $\widehat{\phi}_{JWS}$ | $\widehat{\phi}_{CC}$ | $\widehat{\phi}_{CD}$ | $\widehat{\phi}_{MCD}$
---|---|---|---|---
$[1+(\operatorname{i}\\!\omega\tau)^{\alpha}]^{-\beta}$ | $1-[1+(\operatorname{i}\\!\omega\tau)^{-\alpha}]^{-\beta}$ | $[1+(\operatorname{i}\\!\omega\tau)^{\alpha}]^{-1}$ | $(1+\operatorname{i}\\!\omega\tau)^{-\beta}$ | $1-[1+(\operatorname{i}\\!\omega\tau)^{-1}]^{-\beta}$

Table 2:  The spectral functions $\widehat{\phi}(\operatorname{i}\\!\omega)$
of the various type of relaxation models, namely the HN, JWS, CC, CD, and
postulated MCD models.

| Mittag-Leffler function | Meijer G-function
---|---|---
$\phi_{H\\!N}$ | $\tau^{-1}\big{(}\frac{t}{\tau}\big{)}^{\alpha\beta-1}E^{\beta}_{\alpha,\alpha\beta}\big{(}\\!-\big{(}\frac{t}{\tau}\big{)}^{\alpha}\big{)}$ | $(2\pi)^{\mbox{$\frac{1+l}{2}$}-k}\mbox{$\frac{\sqrt{l}k^{\beta}}{t\Gamma(\beta)}$}G^{k,k}_{l+k,k}\left(\frac{l^{l}\tau^{l}}{t^{l}}\Big{|}{\Delta(k,1-\beta),\Delta(l,0)\atop\Delta(k,0)}\right)$
$\phi_{J\\!W\\!S}$ | $\delta(t)-\tau^{-1}\big{(}\mbox{$\frac{t}{\tau}$}\big{)}^{-1}E_{\alpha,0}^{\beta}\big{(}-\big{(}\mbox{$\frac{t}{\tau}$}\big{)}^{\alpha}\big{)}$ | $\delta(t)-(2\pi)^{\mbox{$\frac{1+l}{2}$}-k}\mbox{$\frac{\sqrt{l}k^{\beta}}{t\Gamma(\beta)}$}G^{k,k}_{l+k,k}\left(\frac{l^{l}\tau^{l}}{t^{l}}\Big{|}{\Delta(k,1),\Delta(l,0)\atop\Delta(k,\beta)}\right)$
$\phi_{CC}$ | $\tau^{-1}(t/\tau)^{\alpha-1}E_{\alpha,\alpha}[-(t/\tau)^{\alpha}]$ | $(2\pi)^{\mbox{$\frac{1+l}{2}$}-k}\frac{\sqrt{l}k}{t}G^{k,k}_{l+k,k}\left(\frac{l^{l}\tau^{l}}{t^{l}}\Big{|}{\Delta(k,0),\Delta(l,0)\atop\Delta(k,0)}\right)$
$\phi_{C\\!D}$ | $\tau^{-1}(t/\tau)^{\beta-1}E_{1,\beta}^{\beta}[-(t/\tau)]$ | $(t/\tau)^{\beta-1}\operatorname{e}^{-t/\tau}/[\tau\Gamma(\beta)]$
$\phi_{M\\!C\\!D}$ | $\delta(t)-\tau^{-1}\big{(}\mbox{$\frac{t}{\tau}$}\big{)}^{-1}E_{\alpha,0}\big{(}-\big{(}\mbox{$\frac{t}{\tau}$}\big{)}^{\alpha}\big{)}$ | $\delta(t)+\beta\tau^{\beta}(\tau+t)^{-1-\beta}$

Table 3:  The response functions $\phi(t)$ presented for the various types of
relaxation models, namely HN, JWS, CC, CD, and MCD, in the language of Mittag-
Leffler function and Meijer G-function.

| Mittag-Leffler function | Meijer G-function
---|---|---
$n_{H\\!N}$ | $1-\big{(}\mbox{$\frac{t}{\tau}$}\big{)}^{\alpha\beta}E_{\alpha,1+\alpha\beta}^{\beta}\big{(}\\!-(\mbox{$\frac{t}{\tau}$})^{\alpha}\big{)}$ | $1-(2\pi)^{\frac{1+l}{2}-k}\,\frac{k^{\beta}l^{-1/2}}{\Gamma(\beta)}G_{l+k,k}^{k,k}\left(\frac{l^{l}\tau^{l}}{t^{l}}\Big{|}{\Delta(k,1-\beta),\Delta(l,1)\atop\Delta(k,0)}\right)$
$n_{J\\!W\\!S}$ | $E_{\alpha,1}^{\beta}\big{(}-(t/\tau)^{\alpha}\big{)}$ | $(2\pi)^{\frac{1+l}{2}-k}\,\frac{k^{\beta}l^{-1/2}}{\Gamma(\beta)}G_{l+k,k}^{k,k}\left(\frac{l^{l}\tau^{l}}{t^{l}}\Big{|}{\Delta(k,1),\Delta(l,1)\atop\Delta(k,\beta)}\right)$
$n_{CC}$ | $E_{\alpha}[-(t/\tau)^{\alpha}]$ | $(2\pi)^{\mbox{$\frac{1+l}{2}$}-k}\frac{\sqrt{l}k}{t}G^{k,k}_{l+k,k}\left(\frac{l^{l}\tau^{l}}{t^{l}}\Big{|}{\Delta(k,0),\Delta(l,0)\atop\Delta(k,0)}\right)$
$n_{C\\!D}$ | $1-E^{\beta}_{1,1+\beta}(-t/\tau)$ | $\Gamma(\beta,t/\tau)/\Gamma(\beta)$
$n_{M\\!C\\!D}$ | $E_{1,1}^{\beta}(-t/\tau)$ | ${{}_{1}F_{1}}\left({\beta\atop 1};-\frac{t}{\tau}\right)$

Table 4:  It presents the relaxation functions $n(t)$ for the various types of
relaxation models, namely HN, JWS, CC, CD, and MCD, in the language of Mittag-
Leffler function and Meijer G-function.

model | evolution equation
---|---
HN | ${{}^{C\\!}(D^{\alpha}+\tau^{-\alpha})^{\beta}}n_{H\\!N}(\alpha,\beta;t)=-\tau^{-\alpha\beta}$
JWS | $(D^{\alpha}+\tau^{-\alpha})^{\beta}n_{J\\!W\\!S}(\alpha,\beta;t)=t^{-\alpha\beta}/\Gamma(1-\alpha\beta)$
CC | ${{}^{c}D^{\alpha}}n_{CC}(\alpha;t)=-\tau^{-\alpha}n_{CC}(\alpha;t)$
CD | ${{}^{C\\!}\big{(}\frac{\operatorname{d\\!}}{\operatorname{d\\!}t}+\frac{1}{\tau}\big{)}^{\\!\beta}}n_{C\\!D}(\beta;t)=-\tau^{-\beta}$
MCD | $\big{(}\frac{\operatorname{d\\!}}{\operatorname{d\\!}t}+\frac{1}{\tau}\big{)}^{\\!\beta}n_{M\\!C\\!D}(\beta;t)=t^{-\beta}/\Gamma(1-\beta)$
Table 5:  Table 5 presents the evolution equation for the various types of
relaxation models, namely HN, JWS, CC, CD, and MCD. The relaxation function
$n(t)$ at the initial time $t_{0}=0$ is equal to 1. The pseudo-operators
${{}^{C\\!}(D^{\alpha}+\tau^{-\alpha})^{\beta}}$ as well as
$(D^{\alpha}+\tau^{-\alpha})^{\beta}$ are, respectively, defined in Eqs.
(V.22) and (VI.21).
model | $h(\xi)$ | $f(\xi,t)$
---|---|---
HN | $n_{D}(\xi)$ |  $f_{H\\!N}(\alpha,\beta;\xi,t)$ given by Eq. (V.24)
| $n_{C\\!D}(\beta;\xi)$ |  $f(\alpha;\xi,t)$ given by Eq. (V.30)
JWS | $n_{D}(\xi)$ |  $f_{J\\!W\\!S}(\alpha,\beta;\xi,t)$ given by Eq. (VI.23)
| $n_{M\\!C\\!D}(\beta;\xi)$ | $f(\alpha;\xi,t)$
CC | $n_{D}(\xi)$ | $f(\alpha;\xi,t)$
| $n_{CC}(\xi)$ | $\delta(t-\xi)$
CD | $n_{D}(\xi)$ | $f_{C\\!D}(\beta;\xi,t)=f_{H\\!N}(1,\beta;\xi,t)$
| $n_{C\\!D}(\beta;\xi)$ | $\delta(t-\xi)$
MDC | $n_{D}(\xi)$ | $f_{M\\!C\\!D}(\beta;\xi,t)=f_{J\\!W\\!S}(1,\beta;\xi,t)$
| $n_{M\\!C\\!D}(\beta;\xi)$ | $\delta(t-\xi)$
Table 6:  The $h(\xi)$ and $f(\xi,t)$ appeared in Eq. (III.12) expressing
subordination approach. Observe that for HN and JWS model we present two non-
trivial types of subordinations.

##  VIII Outlook

The standard relaxation models exemplified by the HN and JWS patterns depend
on three material dependent parameters $\alpha$, $\beta$, and $\tau$, each
time adjusted to the experimental data. These models satisfactorily describe
relaxation phenomena characterized by one-peak (unimodal) behavior in the
frequency domain. However there exist materials which for frequencies $\gtrsim
10^{5}$ Hz exhibit either slower decay or multi-peak behavior of
polarizability. This means that so far studied simple relaxation models do not
cover all possibilities thoroughly enough. As prospective challenger models
used to describe experimentally more complex phenomena we mention the excess
wing (EW) model KGorska21CNSNCb ; RHilfer17 ; RRNigmatulin16  or models
involving sums of standards relaxation patterns taken with different
parameters Liu20 . We remark that the EW model preserves complete monotonicity
in the time domain while the second approach may lead to the lack of this
property and consequently to abandon mathematical methods which result in the
complete monotonicity concepts and reflect seemingly imposing, but far from
complete, physical interpretation of relaxation in terms of summing up the
Debye decays.