## 1 Introduction  
Fractional calculus has been used in modeling of electromagnetic phenomena for many years. In particular, time-fractional (TF) derivatives allow for modeling materials with memory effects and energy dissipation [1]. Nonlocality is a crucial property of fractional-order (FO) operators as it allows for including the effects of the electromagnetic-field history on the current response of the media [2]. Notably, the FO models of dielectric media, such as the Cole-Cole or Davidson models, can be used to describe some biological tissues more accurately and over a broader frequency range than the integer-order (IO) ones, e.g., the Debye model [3]. Fractional derivatives can also be applied in modeling of capacitors, which has been experimentally confirmed to be more accurate than the IO classical approach [4, 5]. Spatial fractional derivatives, in the form of a fractional curl operator [6], are applied to the description of chiral media [7]. Another possible application of the FO calculus is the formulation of fractional multipoles, which has proved to be useful in solving some electrostatic problems, especially those including perfectly conducting wedges and cones [8]. The finite-difference time-domain (FDTD) method [9] remains one of the most popular techniques of electromagnetic simulation, which is successfully applied to FO models of media. For the Debye, Drude and Lorentz media, whose relative complex permittivity is described by a rational function involving integer powers of \(j\omega\), the time-domain expression of the constitutive relation can easily be formulated and implemented inside the FDTD algorithm [10]. On the other hand, the dielectric response involving non-integer powers of \(j\omega\) requires implementation of fractional derivatives in the time domain, which introduces difficulties in terms of domain discretization and computational overhead. In the early approach [11; 12], dielectrics with fractional relaxation of the Cole-Cole, Cole-Davidson and Havriliak-Negami models are implemented based on approximation of the dielectric response by the Debye functions. Unfortunately, this technique does not correctly represent the singularity of time-domain susceptibility, which arises from dielectric models with fractional relaxation [13]. Hence, it is not possible, with the use of this technique, to accurately compute the impulse or step response of spatially complex dispersive scatterers. Other approaches employ the Z-transformation for implementation of the Cole-Cole, Cole-Davidson and Havriliak-Negami models of dispersive media. The method based on the Z-transformation and a second-order Taylor approximation of the Cole--Cole formula is proposed in [14]. This approach is then applied to the problem of the penetration of a short electromagnetic pulse into biological matter. The algorithm of the Cole-Cole model implementation in FDTD [15] is developed based on the Z-transformation and the bilinear transformation. The FO differentiators in this model are approximated by polynomials, whose coefficients are found using a least-squares fitting method. In the approach [16], fractional derivatives in the time-domain representation of the dielectric response are circumvented by the use of fast inverse Laplace transformation and [PERSON]'s method. The Laplace transformation is used to transform the frequency-domain response of a medium into the time domain, and [PERSON]'s method is employed to extract the parameters and transform further the time-domain response into the Z-domain in order to incorporate it directly into the FDTD method. The obtained numerical results are found to be in good agreement with those obtained by an analytical method over a broad frequency range, demonstrating the validity of this FDTD scheme. In other approaches [17; 18], the FDTD simulations of the Cole-Cole media employing the auxiliary differential equation (ADE) are proposed.