Fractional calculus (FC) was introduced more than 300 years ago as a generalization of classical derivative and integral definitions. It is receiving increasing attention for a growing number of applications in different sciences such as physics, biology, chemistry, engineering, finance, mechanics, optics and, in particular, for modeling physical phenomena related to non-Markovian processes, signal and image processing, dielectric relaxation, viscoelasticity, electromagnetism, control theory, pharmacokinetics, fluids, heat transfer, and so on [1,2,3,4,5,6,7].

It is well known that the description of real-world problems using mathematical models based on the differentiation and integration provides a huge contribution to science. However, the classical PDE models fall short of describing many complex anomalous systems and phenomena featuring persistent memory effects and long-range interaction in the medium. The integer-order differential and integral operators are local, i.e., the interactions between two domains happen only through the contact. On the other hand, the space and time-fractional integro-differential operators are not local. This means that interactions among components extend to a neighborhood of each component (space non-locality) as well as being able to model systems in which the reaction to an external excitation is not instantaneous but depends on the history of the system (time non-locality). For this reason, the integral and differential operators of non-integer order are more effective tools for describing media and systems with non-local and hereditary properties of power-law type, long-range memory and/or fractal properties. Moreover, the unusual properties of the fractional derivatives and integrals, including the violation of the standard product, chain, and semigroup rules, can be used to describe unusual properties of complex systems and media [8].

As it is well known, Maxwell’s equations are the theoretical background for understanding the classical electromagnetic phenomena. In their differential form, Maxwell’s equations involve integer-order calculus applied to vector differential operators of the first order. Therefore, the extension of the fractional paradigm in electromagnetism needs a consistent FVC theory. During the last few decades, different fractional generalizations of vector operators such as gradient, divergence, curl, and Laplacian have been suggested. The potential applications of fractional calculus in electromagnetic theory were explored in [9,10] with focus on multipole expansion, theory of images in electrostatics, scalar Helmholtz equation. The definition of the FO curl operator in electromagnetic theory was introduced in [11] to derive solutions for Maxwell’s equations by fractionalizing the principle of duality. This approach was applied to solve some reflection and diffraction problems [12]. The calculation of the source distributions corresponding to fractional dual solutions and intermediate fractional dual solutions was illustrated in [13]. The behavior of fractional dual solutions in metamaterials with negative permittivity and permeability, as well as homogeneous and non-homogeneous chiral media was investigated in [14,15,16] using a higher-order fractional curl operator. The definition of FO gradient and FO divergence were proposed in 1998 [17] and 2006 [18], respectively, as well as a consistent formulation of FVC based on Riemann–Liouville integration and the Caputo differentiation was suggested in [19]. Using this approach, the fractional differential and integral vector operators were defined and the fractional Green, Stokes, and Gauss theorems were proved. After 2008, other articles illustrating special aspects of FVC began to appear in the scientific literature. Using Grünwald–Letnikov fractional calculus, the classical vectorial operators of gradient, divergence, curl, and Laplacian could be generalized to fractional orders and an extension of the Helmholtz decomposition theorem for both fractional time and space was proposed [20]. Based on these results, the fractional versions of the classical Green, Stokes, and Ostrogradski–Gauss theorems were introduced in [21]. Moreover, higher-order fractional Green and Gauss formulas were derived in [22]. However, just in 2021, a generalization of FVC was proposed to include a general form of non-locality in kernels of fractional vector differential and integral operators [23].

Thanks to the aforementioned developments in FVC, different approaches have been suggested to fractionalize the electromagnetic field equations. The idea of such a generalization relies on the application of material constitutive relations extended in terms of FO derivatives. Using the fractional differential forms, the non-local Maxwell’s equations were introduced in [24]. The authors (i) defined fractional vector and scalar potentials, (ii) developed fractional conservative law for the electric charge, (iii) derived fractional wave equations and fractional Poynting theorem. A solution to the source-free fractional wave equation in isotropic and homogeneous dielectric media was presented in [25] with the aim of maintaining the physical units of the system for any value of the fractional operator exponent. Moreover, it was highlighted that the resulting fractional space–time waves are independent of the physical structure, chemical composition, and polarization of the material. The fractional wave equation in space-time and in a conducting media was introduced and analyzed in [26]. The authors showed that fractional electromagnetic fields have different characteristics compared to the classical ones and they exhibit anomalous behavior that is known in the literature as centrovelocity. Solutions for fractional generalization of the Laplacian were derived in [27] with the aim of describing the electric fields in non-local media with power-law spatial dispersion. In particular, a generalizations of Debye’s permittivity and a weak spatial dispersion of power-law type in plasma-like media were suggested. Using the formulation of TF electrodynamics with the Riemann–Silberstein vector, a compact form of FO Maxwell equations was presented in [28]. The proposed formulation allows for inclusion of energy dissipation and establishes a relation between the TF Schrödinger equation and TF electrodynamics. Using the two-sided Ortigueira–Machado derivative, the plane wave signal propagation in media described by FO model was analyzed in terms of system causality [29]. In particular, causality of the transfer function is proven for certain values of the asymmetry parameter, corresponding to the left-sided Grünwald–Letnikov derivative. In addition, an exact solution of the spherical and cylindrical wave equations in fractional dimensional space was presented in [30,31].

The modeling of electromagnetic wave propagation, radiation, and scattering in complex fractal structures has attracted growing attention and seen an increasing number of applications. Moreover, fractal models are becoming popular as they allow describing media of great complexity and rich structure using just a few parameters. Starting from the idea that isotropic and anisotropic fractal media in an Euclidean space can be replaced by some fractional continuous mathematical model [32,33], the electromagnetic wave propagation in fractal media can be studied through an ordinary continuum model with non-integer dimensional space. Using the concepts of power-law density of states and fractional-order integration, fractal electrodynamics of fractal distribution of charges, currents, and fields was described in [4,34]. Electromagnetic equations in fractional space were presented in [35,36]. In particular, a modified vector differential operator in a D-dimensional fractional space was defined with the aim of introducing a fractional space generalization of differential electromagnetic equations in a far-field region. In [37], the wave propagation in fractal media was studied by solving a generalized vector Helmholtz equation in terms of plane wave solutions. In [38], the authors derived electromagnetic equations in generally anisotropic fractal media using a dimensional regularization approach. Such an approach allows to define fractal gradient, divergence, and curl operators, satisfying the four basic identities of vector calculus as well as the Helmholtz decomposition theorem. Furthermore, it readily leads to the generalization of Green–Gauss and Stokes theorems, as well as of the charge conservation equation in anisotropic fractals. Again, the study of transmission and reflection of electromagnetic waves at dielectric–fractal interfaces was presented in [39]. The authors assumed that the fractionality exists only along the z-axis and that the permeability of the fractional medium is approximately the same as in the integer space. Therefore, they derived the expressions for transmission and reflection coefficients for parallel and perpendicular polarizations, as well as the general solution for plane waves.

Dielectric materials have an important role in science and technology. Hence, their relaxation properties are the subject of numerous new research activities. For instance, the dielectric response knowledge of an unknown material inside a lossy and dispersive structures is of paramount importance in GPR techniques and material characterization [40,41,42]. The accurate knowledge of frequency-dependent dielectric susceptibility is also relevant in microwave computational dosimetry, dielectric spectroscopy and imaging for early-stage cancer diagnostics and treatment, as well as in the study of the materials involved in the production and processing of crops and food of agricultural origin [43,44]. Furthermore, dielectric properties play an important role in the analysis of body area networks, nano-networks in living biological tissues, and for in-body electromagnetic communications [45,46,47,48].

Great efforts have been made to model the frequency dispersive behavior of the dielectric susceptibility for a variety of materials in different frequency ranges. Simple models based on exponential decay in time of the relaxation function were introduced in [49,50,51]. Numerous experimental studies have demonstrated that anomalous relaxation and diffusion processes can exhibit a time behavior that is different from the exponential law. These dynamic dielectric properties seem to be a feature of disordered media as spin or dipole glasses, polymers, biopolymers, emulsions, and microemulsions, disordered ferroelectrics, biological cells and tissues, porous materials [52]. Thanks to the relevant advances in measurement science, a significant amount of experimental data has been collected to validate the existence of non-exponential relaxation processes. Meanwhile, the attention of researchers has been focused on modeling dielectric response over broader frequency ranges [51,52,53]. The resulting models allow evaluating performance indexes as well as optimizing the dielectric response so as to meet given design specifications or to generate unusual electromagnetic properties [54,55]. The right knowledge of the dielectric response over broad frequency ranges is of utmost importance for the development of accurate theoretical models and computational techniques aimed at determining the electromagnetic field propagation properties inside complex and disordered dielectric materials. In particular, numerical simulations could be very useful for identifying the main parameters involved in non-invasive diagnosis and medical sensors, to evaluate the characteristics of in-body communications among nano devices in terms of data rate and transmission range [48], as well as to provide guidelines for computational dosimetry and specific therapeutic approaches.

Important empirical relaxation laws or relationships to describe non-exponential relaxation processes causing broadness, asymmetry, and dielectric dispersion excess were introduced in [56,57,58,59,60,61,62,63]. Even though these relationships have been moderately successful in describing relaxation phenomena in different condensed matter systems, they do not entirely clarify the anomalies in the dielectric response of disordered systems from the physical point of view. With the help of a model based on a self-similar relaxation process, in [64], the authors attempted to understand the physical process leading to the appearance of non-exponential relaxation in the Cole–Davidson expression as well as the meaning of the non-integer parameter which appears in it. They stressed that contrary to the TRD concept, the self-similar one results in a relaxation with discontinuous behavior, such that the times when the interaction with an external field exists are distributed over some fractal set. In [65], the authors derived different relaxation functions assuming relaxation times bounded between the upper and lower limits of self-similarity. Moreover, they predicted that at times shorter than that corresponding to the lowest self-similarity level, the relaxation should be of Debye-like type, whatever the pattern of non-classical relaxation is at longer times. In [66], the author attempted to provide a unified interpretation of general relaxation phenomena by exploring the physical mechanism underlying dielectric relaxation phenomena. In particular, it was demonstrated that dielectric materials possess the normal structure and the nematic phase in which the slowly fluctuating dielectric response is influenced and constrained by the collective one.

Even if there is no universal law capable of describing any dielectric response in any frequency band, the proposed empirical relationships share the common feature of power law decay for large times. This behavior allows deriving integro-differential equations with time derivatives and integrals of non-integer order that can be incorporated into Maxwell’s equations. The use of fractional derivative techniques to describe the anomalous dielectric response, in the time domain, of an inhomogeneous medium is demonstrated in [65]. Here, the authors established the relationship between anomalous relaxation and dimensionality of a temporal fractal ensemble characterizing the non-equilibrium state of a medium. An equation for the relaxation function containing FO integral and differential operators was derived and solved in [64]. A model based on viscoelastic analysis techniques employing FO operators was proposed in [67] to evaluate the dielectric response of a material in both radio frequency and terahertz bandwidths. Moreover, a survey on the main dielectric models and their characterization in terms of FO differential operators can also be found in [7]. In [68], it was proved that electromagnetic propagation in a wide class of dielectric media can be described by FO differential equations in the time domain. On the basis of recent non-local effects observed in EMMs, a fractional Drude model was introduced to describe electromagnetic characteristics of metamaterials [69]. A scientific discussion on the analysis of perfect cloaking based on FO formulation of electromagnetics was initiated in [70]. In particular, a two-dimensional cloak was numerically demonstrated in media described by FO models. Moreover, a tuning of the signal arrival time was obtained by means of a small perturbation of the time-derivative orders in Maxwell’s equations [71].

Other key issues which should not be neglected pertain to the development and implementation of numerical methods for finding approximate solutions of FO Maxwell’s equations, especially when fractional problems involving large spaces and long-time iterations need to be solved. In fact, due to the non-locality of the fractional operators, the computational cost concerning the numerical implementation of fractional differential equations is much heavier than that relevant to classical integer-order equations. For example, the relative complex permittivity functions corresponding to non-integer powers of jω