Formulation of Time-Fractional Electrodynamics Based on Riemann-Silberstein Vector

In this paper, the formulation of time-fractional (TF) electrodynamics is derived based on the Riemann-Silberstein (RS) vector. With the use of this vector and fractional-order derivatives, one can write TF Maxwell’s equations in a compact form, which allows for modelling of energy dissipation and dynamics of electromagnetic systems with memory. Therefore, we formulate TF Maxwell’s equations using the RS vector and analyse their properties from the point of view of classical electrodynamics, i.e., energy and momentum conservation, reciprocity, causality. Afterwards, we derive classical solutions for wave-propagation problems, assuming helical, spherical, and cylindrical symmetries of solutions. The results are supported by numerical simulations and their analysis. Discussion of relations between the TF Schrödinger equation and TF electrodynamics is included as well.


Introduction
Recently, the problem of electrodynamics of fractional order (FO) has been introduced in the literature. Between various approaches, one can notice the generalizations of Maxwell's equations involving spatial and temporal FO derivatives [1][2][3][4][5], only spatial derivatives [6][7][8], and only temporal FO derivatives [9][10][11][12]. Such generalizations of Maxwell's equations employ FO derivatives in order to describe the dynamics of electromagnetic systems with memory and energy dissipation [13]. Furthermore, the use of FO derivatives in electrodynamics allows one to define the electric potential based on the concept of FO poles [14][15][16]. Hence, a concept of the fractional multipole is formulated, which allows for calculating the electric potential at the distance r from the source, which is proportional to r −α (where α ∈ R is not integer). This approach can be applied to solve some electrostatic-image problems, involving perfectly conducting wedges and cones. In [6], the fractionalization of the duality principle in electrodynamics is proposed, making use of the FO curl operator. This approach can be applied to solve some reflection [17,18] and diffraction [19] problems. Afterwards, in [20], the FO curl operator is proven to be useful for the analysis of the transmission of electromagnetic waves through a slab of reciprocal chiral medium. Hence, the FO curl operator, to some extent, appears to be mathematically equivalent to a chiral medium. As a consequence, some electromagnetic characteristics of media can be modelled using the FO derivatives.
In this paper, we extend the formulation of time-fractional (TF) electrodynamics with the use of the Riemann-Silberstein (RS) vector, i.e., the complex vector which combines simultaneously the electric and magnetic field. This vector was initially introduced by Bernhard Riemann in the lectures on partial-differential equations, which were edited and published in 1901 by Weber [21]. Then, Silberstein studied the RS vector [22,23] and employed it for the quaternionic formulation of the relativity theory [24]. For the historical background of the RS vector, one is referred to [25]. In [25][26][27][28], it is demonstrated that the RS vector provides a uniform and elegant description of electromagnetism, which connects classical and quantum electrodynamics. In [29,30], it is demonstrated that the RS vector is useful for the analysis of antenna and propagation problems in engineering.
To the best of Authors' knowledge, none of formulations of FO electrodynamics has been extended yet using the RS vector. Furthermore, the existing formulations of electrodynamics based on the RS vector do not include the energy dissipation, i.e., they are limited to a homogeneous and static medium, which is usually a vacuum. However, the time-domain fractionalization of Maxwell's equations allows for inclusion of energy dissipation, which is the intrinsic property of media described by fractional-order models (FOMs). Furthermore, it is possible to demonstrate that the energy density of the electromagnetic field and media described by FOM, which is spent to create the field and change the momentum of medium, reflects the entire previous history of the process. Hence, the application of FOM allows for description of memory effects in electrodynamics. Therefore, we investigate this topic focusing on classical-physics solutions. Then, we can demonstrate relations between the TF Schrödinger equation and TF electrodynamics using the proposed formulation based on the RS vector.

Preliminaries
In this section, we introduce the basic notation and terminology, which is used throughout the paper.

Basic Notations
The standard engineering notation for the imaginary unit j = √ −1 is used throughout the paper. We denote real and imaginary parts of the complex number s by (s) and (s), respectively. Then, the conjugate to the complex number s is denoted as s * .
The Fourier transformation is defined for the absolutely integrable function of time f : R → C (i.e., f (t)) by the formula (1) whilst the inverse Fourier transformation is given by the formula For the absolutely integrable function of space g : R 3 → C (i.e., g(r)), we define the spatial Fourier transformation by the formula F (g)(k) = G(k) =ˆR 3 e jk·r g(r)d 3 r whilst the inverse spatial Fourier transformation is given by the formula The convention of the Fourier transformation can be changed into to the one applied in physics and mathematics. For this purpose, the imaginary unit j should be replaced by −i (where i = √ −1) in (1)-(4). Throughout the paper, a branch of logarithm with the domain covering all purely imaginary numbers of the form jω is taken, where ω ∈ R \ {0} ln(s) = ln(|s|e jθ ) = ln |s| + jθ (5) for s = |s|e jθ , s = 0 and θ ∈ (−π, π). From (5), one obtains ln(jω) = ln |ω| + j π 2 sgn(ω).

FO Calculus
The Marchaud derivative of the order β ∈ (0, 1) is defined as Although the Marchaud derivative can be formulated for higher orders β > 1 (see Section 5.6 in [31]), we do not refer to this more general case in our investigations. In (10) and (11), the function f is assumed to be smooth enough, e.g., f ∈ C 1 (R) with | f | bounded by some function not growing too quickly in ±∞ (more about these assumptions later). Applicability of the Marchaud derivative in our investigations stems from its advantages, refer to [32,33]: • The Marchaud derivative of the order β ∈ (0, 1) of the function f exists, if f ∈ C 1 (R) and f (t) = O(|t| β−1−ε ) as |t| → +∞ for some ε > 0 (see Section 5.4 in [31]). If additionally f ∈ L p (R) for p ∈ [1, 1/β), then the Marchaud derivative coincides with the Riemann-Liouville derivative with a base point a = −∞. • For the derivative of the order α of the exponent e st (where s ∈ C is fixed), one obtains for (s) ≥ 0. As discussed in [34] (see (5) therein), this limitation stems from the memory of derivative, which covers the entire interval (−∞, t). Therefore, allowing for (s) < 0 would lead to divergent integrals. In practical terms, it is not a limitation, because one should consider functions bounded at −∞, or even causal (i.e., equal to zero for arguments in (−∞, 0). • The Marchaud derivative, the Riemann-Liouville derivative with a base point in −∞ and the Grünwald-Letnikov derivative [31,35] are equivalent for a very broad class of functions. • When the Marchaud derivative coincides with the Grünwald-Letnikov derivative for the function f , then one obtains where α, β > 0. • In some cases, we would like to assume not only that the Marchaud derivative of the function f : R → R exists, but also that this derivative belongs to an appropriate space L p (R). It is handled by Theorems 6.1 and 6.5 in [31]. The condition that guarantees that and lim t→+∞ f (t) = 0 for λ > max{β, 1/p − β}. The function space H λ (R) denotes the space of Hölder continuous functions with appropriate behaviour in ±∞ (see (1.5) and (1.6) in [31]). One can check that, in order to belong to H λ (R), it is enough to assume that f ∈ C 1 (R) and |t|≥1 | f (t)| q |t| 2(q−1) dt < +∞ (14) for q ≥ 1 1−λ . The last condition is satisfied, e.g., when We employ in our considerations some properties of the null space of the Marchaud derivative. We are not aware of any general result in this area, but we can show some reasonable behaviour assuming that φ is a Hölder continuous function. That is, if one assumes that φ ∈ H λ (R) with exponent λ ∈ (β, 1) and lim t→±∞ φ(t) = 0, then D α φ = 0 implies φ = 0. This is a consequence of Theorems 6.1 and 6.5 from [31]. • The appropriate Fourier transformation identity can be formulated for the Marchaud derivative (see formula (7.4) in [31]) when f ∈ C n (R), α ∈ (n − 1, n), and f along with derivatives up to the order n belongs to L 1 (R). Due to these properties, the Marchaud derivative is preferable in our investigations focused on electrodynamics. It occurs that popular definitions of FO derivative, e.g., Riemann-Liouville and Caputo [31,35], do not satisfy the properties (12) and (13) in general, which are necessary for the analysis of the wave propagation in the media described by FOM (refer to [32,33]).

Vector Fields
It is necessary to formulate precise assumptions related to functions (i.e., vector fields), which we use in our considerations. Suppose now that V ⊂ R 3 is a compact volume with the boundary S being the piecewise smooth surface. For the vector field X : V × R → R 3 (or X : V × R → C 3 ), we always assume that this is C 2 map in the entire domain and the functions X i (x, y, z, ·) (i ∈ {x, y, z}) belong to L 1 (R) for all (x, y, z) ∈ V. Hence, the Fourier transformation (1) can be calculated with respect to the time variablẽ for a fixed (x, y, z) ∈ V. We also assume that the limits as |t| → +∞ of functions X i (x, y, z, ·) (i ∈ {x, y, z}) equal to zero and the time derivative D 1 t X i (x, y, z, ·) = O(t −c ) as |t| → +∞ for c > 0 big enough. It implies that the derivative D β t X(x, y, z, ·) ∈ L 1 (R). One should note that having two functions f , g : R → R, satisfying the assumptions given above for which D which are well defined and differentiable. In (17) and (18), I 1 t denotes the integral of upper integration limit´t −∞ (·) dτ and I In what follows, we often refer to the dot-product notation. For complex vectors A, B ∈ C 3 , the product does not refer to the standard dot product over the field of complex numbers but to the bilinear form 'borrowed' from the real-valued case, i.e., Therefore, to be formal, when we want to refer to the standard dot product of complex vectors A, B ∈ C 3 , we write:

RS Vector
In this section, we propose TF formulation of electrodynamics based on the RS vector. Our FO approach is consistent with the RS vector formulation of classical electrodynamics proposed in [25,29].
Let us consider symmetric Maxwell's equations in free space where E and H are, respectively, the electric-and magnetic-field intensities, D and B are, respectively, the electric-and magnetic-flux densities, J and M are, respectively, the electricand magnetic-current densities, and ρ e and ρ m are, respectively, the electric-and magneticcharge densities. We assume the constitutive relations for the medium described by FOM in the following form [34,36]: In (27) and (28), the parameters β and µ γ denote, respectively, FO equivalents of permittivity and permeability, whose SI units are as follows: In [37], the constitutive relation (27) is described as an 'engineering' model of dielectrics. Such a model describes memory effects and hereditary properties of electromagnetic field in media. It results not only from experimental results of Westerlund and Ekstam published in 1994, but also from the purely empiric Curie's law formulated in 1889 [38]. Assuming analogous form of the constitutive relation for the magnetic field, one can write (28). When β = 1 and γ = 1, the constitutive relations for the linear medium, described by integer-order (IO) model (IOM) with the permittivity 1 and the permeability µ 1 , are obtained For a vacuum, we assume that the permittivity and permeability are respectively denoted by 0 and µ 0 . Comparing to the settings considered before [34], we assume the same FO of all the time derivatives in (23)- (26) in order to formulate TF Maxwell's equations based on the RS vector. In our considerations, we assume that 1 2 ≤ β = γ ≤ 1. Therefore, one can write (23)- (26) as Let us now define the parameter whose the SI unit is Ohm, so may be considered as a form of the wave impedance. Let us formulate the RS vector as follows: Our definition of the RS vector is consistent with the definition applied in engineering [29,30], where the E, H vectors are usually employed to describe electromagnetic systems (e.g., antennas and waveguides). In order to write Maxwell's equations for the RS vector F, we need to introduce an additional notation. Let us define For a vacuum when β = 1, one obtains c 1 = c = (µ 0 0 ) − 1 2 . With the definitions given above, FO Maxwell's Equations (31)-(34) can be written as Without sources (K = 0), (40) can be written as the evolution equation analogous to the TF Schrödinger equation [39][40][41] with the HamiltonianĤ = c β (p · Σ), i.e., where p = −jh∇ and Σ denotes the spin-1 counterparts of the Pauli matrices The TF Schrödinger equationĤψ = jhD β t ψ differs from the ordinary Schrödinger equation due to the TF derivative of the wave function ψ. The ordinary Schrödinger equationĤψ = jh ∂ ∂t ψ has a mathematical form of a diffusion equation and can be obtained for a quantummechanical particle by considering a Gaussian-probability distribution in the space of all possible paths [39,42]. Hence, one can obtain various types of Schrödinger equations for non-Gaussian distributions. For instance, Feynman's path-integral approach is employed in [43] to construct space-fractional quantum mechanics. That is, the Schrödinger equation is obtained with FO space derivatives due to the application of Levy distributions instead of Gaussian distributions to the set of possible paths.
Similarly, the TF Schrödinger equation is obtained if one considers non-Markovian evolution [39,40]. Analogously, the TF Maxwell's Equations (31)- (34), which can be formulated in the form similar to the TF Schrödinger equation, actually stem from the application of non-standard memory function in a transport equation in order to include a temporal dispersion [9].
In general, the FO derivatives stem from relaxation and oscillation processes that exhibit memory and delay. For instance, the FO derivative in the oscillator equation is equivalent, to some extent, to the presence of a dissipative term in the standard oscillator equation [44]. However, this damping represents an intrinsic feature of the motion equation and does not result from additional dissipation like for a damped harmonic oscillator. Similar dumping of dynamics is also visible for the harmonic oscillator described by the TF Schrödinger equation [41]. Analogously, as it is shown later, the application of non-local formalism of FO derivatives to describe complex electromagnetic systems with memory also provides the energy dissipation, which is the intrinsic property of such systems.
The TF Schrödinger equation describes a non-relativistic particle [40]. Analogously, the application of FO constitutive relations (27) and (28) to Maxwell's equations provides non-relativistic solutions obtained for the diffusion-wave equation [45][46][47]. This equation delivers solutions classified as intermidiate cases between the diffusion and the wave processes, which are different in terms of response to a localized disturbance. That is, the diffusion equation models a disturbance which spreads infinitely fast, whereas the wave equation models a disturbance with constant propagation velocity. Hence, for the TF diffusion-wave equation, the fundamental solution possesses a maximum that spreads with a finite velocity and the propagation velocity of a disturbance is infinite [45]. To sum up, the TF diffusion-wave equation is non-relativistic similarly to the standard diffusion equation.
Using the definition of the RS vector (36), one can find that Poynting's vector [48,49] can be written as Then, one can find that It has the same unit as the Poynting vector, i.e., [F · F * ] = Watt/m 2 and can be interpreted as the power density of the plane wave propagating in a medium whose the wave impedance is Z f . Then, using the definition (36), one can formulate inverse formulas for the electromagnetic-field vectors Furthermore, because of the assumptions (19) and (20), one can define the vector G such as: With the use of power-law property (13), (48) can be written as Then, assuming that the inverse operator exists for the FO differentiation, (48) can equivalently be written as For assumptions given in Section 2.3, the relation between F and G is unambiguous. Therefore, one can write Maxwell's Equations (23)-(26) as follows: Hence, based on (50), one can formulate inverse formulas for the electric-and magnetic-flux densities D and B. That is With the use of (53) and (54), one can define the electromagnetic-momentum density based on Minkowski's approach [50,51], i.e.,

Plane-Wave Propagation
Let us consider the plane-wave propagation in a free space without sources (ρ = 0 and K = 0). Let us employ the spatial Fourier transformation (3) to the RS vector, i.e., F(k, t) = F (F(r, t)). For the sake of brevity, the same symbol is employed to denote the Fourier transform, but with the argument k. Application of the spatial Fourier transformation (3) to (40) and (41) gives It implies that k ⊥ F(k, t). Let us introduce the polarisation vectors e ± (k) such as where k = |k|. Eigenvectors of the cross-product (58) differ only by the phase factors. Let us denote e + (k) = e(k) and notice that e − (k) = e(−k) = e * (k).
In [25], the solution to IO Maxwell's Equations (56) and (57) is represented in the RS formalism as follows: Unfortunately, this is not the case for the considered FO generalization, because components of (59) lead to and It implies that f + (k) = f * − (k) = 0, because none of the complex numbers (jω) β nor (−jω) β is purely imaginary for β ∈ (0, 1).
Hence, the general solution to (56) and (57) is written in a more general form The complex conjugate and the minus sign in the second term of (62) stem from our requirement of consistency with the classical theory of the RS vector [25]. Let us consider the first component in (62), i.e., e + (k) f + (k, t). It satisfies (56); hence, using the relation k × e + (k) = −jke + (k), one obtains Similarly, for the component corresponding to the vector e − (k), one obtains By taking the complex conjugate, the last equation can be written as One should note that for both polarization vectors e + (k), e − (k) and the fixed vector k, a FO differential equation is obtained where λ = λ(k) = jc β k. Among the solutions to (66), one may find non-zero functions that are zero on the interval (−∞, a) and which satisfy the Equation (63) on the interval [a, +∞), when D β t denotes the Riemann-Lioville derivative of the order β with a base point a ∈ R. By Theorem 4.1 in [35] (see also (7.2.10) in [52]), the solutions to (66) for λ ∈ R and t ∈ [a, +∞), are given by where A ∈ C is a constant. In (67), E β,β denotes the Mittag-Leffler function of the second type, i.e., it is the entire complex function given by (see, e.g., (4.1.1) in [52]) .
One should note that in [35,52], the problem with a real value of λ is concerned, but the method of the proof works in the same way also for a complex value of λ. Looking for a causal solution, which is bounded in any set [δ, +∞) ⊂ R, one has to take the solution for To keep the notation shorter, when writing E β,β (jc β kt β ), it is assumed to be equal to 0 for also satisfies (56).
Therefore, the general solution to (40) for K = 0 is given by Taking into consideration that´R 3 g(−x) d 3 x =´R 3 g(x) d 3 x, e + (k) = e(k) and e − (k) = e(−k), one can write (71) as follows: Let us consider the propagation of attenuated waves along the z-direction, i.e., k = k i z and Then, for fixed functions f + (k) and f − (k), one can write (71) as follows (remembering that now the spatial Fourier transformation is applied in one dimension only): Taking the special case of f + (k) = 1 , one obtains from (74) for the x coordinate (F(z, t)) x = 1 2πˆR t β−1 (cos(kz)j (E β,β (jc β kt β )) + sin(kz) (E β,β (jc β kt β ))dk.
Let us now consider real and imaginary parts of t β−1 E β,β (jc β kt β ). As one can see, these may be represented as Hence, the real part of (F(z, t)) x , i.e., the electric-field component, can be written as This is Green's function, which is given in [53] for, so called, signalling problem for the diffusion-wave equation with an appropriate boundary condition (cf. formula (17) in [53], where the normalized problem with c β = 1 is concerned). In [53], Green's function is given by and differs from (78) by a multiplicative constant only.
The formula (80) is one of many equivalent representations including the Formula (72) from [54], which is taken from [45,47], and used in our numerical simulations. These two formulations are actually closely related through the integration by parts formula.
Let us also consider the solution (71), which is the boundary case for β = 1 and c β = c.
. Hence, one obtains classical unattenuated solutions in the form: Consistency with solutions for IO Maxwell's equations confirms correctness of our derivations. Now, one can identify the left-and right-handed circularly polarized attenuated waves propagating in the z-direction, which are related with f + and f − functions within the general solution (72), respectively. In Figure 1, the rotating polarisation vectors e(k)t β−1 e −jk·r E β,β (jc β kt β ) and e(k)t β−1 e jk·r E * β,β (−jc β kt β ) (where k = ki z ) are presented for these components, which create helical surfaces. As seen, for both components related to f + and f − functions, the wave amplitude is attenuated and helices rotate in opposite directions.

Fundamental Solutions
Maxwell's Equations (40) and (41) lead to the following equation and eventually Let us consider the TF diffusion-wave equation in free space without sources Let us solve (84) by separating the variables, i.e., assuming that and where i ∈ {x, y, z} for the Cartesian coordinate system. Each of the components Ψ i (r)T i (t) satisfies the scalar equation as follows: Naturally, it is assumed that Ψ i (r) = 0 and T i (t) = 0 in (87). Since both sides of (87) are functions of different variables, then both sides of (87) are equal to the constant q 2 i . Hence, (84) is equivalent to the following set of equations Equation (89) is satisfied by the functions e jω i t when q 2 i = (jω i ) 2β and e −jω i t when q 2 i = (−jω i ) 2β ; hence, the general solution to the i-th coordinate of (84) can be written as To obtain the solution to the vector Equation (84) from (90), the same frequency (i.e., ω = ω x = ω y = ω z ) for all coordinate components has to be assumed. Then, one obtains whereF − =F − (r) andF + =F + (r). Assuming that ω ∈ R, the solution (91) can be reduced to the phasor representation of the RS vector, i.e., F =Fe jωt whereF =F(r) is a complex vector. Then, one obtains the following stationary equation, which depends on the angular frequency ω only: In the case of plane wave propagating in R 3 , one can assume that the propagation takes place along the z-direction. It corresponds to the vector fieldF given bỹ whereF(z) is the scalar complex-valued function and e z is the corresponding unit vector given by (73). Then, inserting (93) into (92), one obtains the corresponding 1-D equation for an unknown complex-valued functionF(z) whose the solution is given byF where ξ = c −1 β (jω) β and F + and F − are certain complex constants. Hence, for a fixed z, one obtains the scalar transfer function in the frequency domain whereF + ∈ C is a certain complex constant. Its value does not influence the causality of the transform, i.e., it can be demonstrated that the transform e −zc −1 β (jω) β is causal iff β ∈ (0, 1) [55]. However, this transform is not relativistically causal, i.e., when t < 0 then the inverse Fourier transform of e −zc −1 β (jω) β is equal to zero, whilst it is different then zero for e −zc −1 β (jω) β e jωzc −1 . One may also consider the radially symmetric solutions to (92). Hence, let us assume that the solutionF is given bỹ where e R is the corresponding unit vector in spherical system. Then, one obtains where R = x 2 + y 2 + z 2 . Then, the 3-D Equation (92) takes the form and the solution R f (R) is given by It gives For a fixed R > 0, one can write the transfer function as where the normalizing factor appears following Section 2.3.1 of [56] and the dependence on ω has the same form as before. Hence, it is a causal transform iff β ∈ (0, 1). One may also consider the solutionF with cylindrical symmetry, i.e., given bỹ where e r is the correspoding unit vector in cylindrical system. With the assumption that the wave propagates in the xy plane (i.e., with z = const), one obtains where r = x 2 + y 2 . Then, the 3-D Equation (92) takes the form which can be written as Equation (106) is known as the Bessel equation of the order α = 0. Its solutions are known as cylinder or Bessel functions (see Section 3.9 in [57] for a recursive definition of cylinder functions and Section 4.3, Equations (3) and (4), for solutions of the Bessel Equation (105) as well as Sections 8.4 and 8.5 in [58]). As one can notice, (105) corresponds to the Equation (3) in Section 4.3 in Watson's book [57] with c = (jω) β /c β and p = −1/2. Hence, the solutions to (105) are functions given by where Z 0 is any Bessel function of the order ν = 0 (as in Sections 8.4 and 8.5 of [58], or the cylinder function C 0 from Watson's book [57]). As a special case of cylinder functions, one can consider (see 8.401 in [58]) the Bessel functions of the first kind J 0 (z), the Bessel functions of the second kind Y 0 (z) and the Hankel functions H 0 (z) (also called the Bessel functions of the third kind). Each of the pairs (J 0 (j(jω) β r/c β ), Y 0 (j(jω) β r/c β )) 0 (j(jω) β r/c β )) forms a fundamental system of solutions to (105). In the considered case, we take a function that tends to 0 as r → +∞. We assume that f (r) = J 0 (j(jω) β r/c β ), where Bessel function J 0 is given by the series (108) Hence, for a fixed r > 0, the transfer function is given by The Bessel function is multiplied by a factor − j 4 , which comes from the normalization as in Section 2.3.2 in [56] (Hankel functions are used in [56], but one should remember that J 0 (z) = 1 2 (H (1) 0 (z) + H (2) 0 (z)) which makes the condition valid for J 0 as well). To illustrate differences between solutions for plane wave, radial, and cylindrical symmetries, we performed simulations showing the electric field observed at certain moments in time at spatial points of varying distance from the source. Three transfer functions (96), (102) and (109) are considered. The algorithm of computations is described in Section 5 of [34] and verified against accurate solutions [59]. The simulation shows the system response to the unit-impulse excitation (which approximates Dirac's delta distribution) in two time moments t 1 = 1 f s and t 2 = 2 f s at the distance varying from 0.1-1 µm. The simulations are performed assuming the sampling period T s = 1 50·750 10 −12 s. It is assumed that the value of c β is equal to c, but expressed in appropriate units, i.e., m s β . The simulations are performed for β ∈ {1, 0.995, 0.99}. Figure 2 corresponds to the plane-wave propagation, Figure 3 corresponds to the radial-symmetry case, while Figure 4 corresponds to the cylindrical symmetry case. In each figure, for β < 1, an inclined wavefront is observed with the maximum of a pulse propagating with finite velocity. For the radial-symmetry case, the values of the plotted function are much larger due to singularity in the formula (102) for R = 0. As one can notice, the electric-field intensity is different than zero before the arrival of the wavefront maximum. This feature means that fundamental solutions significantly differ from Green's functions in IO electrodynamics. The disturbance spreads infinitely fast, which makes the solutions similar to the ones known for the TF diffusion-wave equation being non-relativistic, similarly to the standard diffusion equation. As one can see in all the cases, the decreasing derivative order (i.e., when systems become 'more fractional'), the amplitude of the system response decreases but the pulse width is increasing. Furthermore, the 'tails' of pulses are increasing, which represent the memory effects of FOMs.

Power Conservation
Let us consider the energy-balance equation (i.e., Poynting's theorem) in TF electrodynamics formulated in the time domain (without electric and magnetic currents and Joule's heating, i.e., J = 0 and M = 0, for the sake of brevity) [60] ∇ · (E × H) Then, using (46) and (47), one obtains this theorem formulated based on the RS vector Equation (111) can be rewritten as is the density of accumulated and dissipated power of the field and the medium described by FOM. Hence, denotes the density of electromagnetic-field energy, which is spent to create the field and change the momentum of medium. This integral exists because F is continuous and bounded as well as D β t F belongs to L 1 (R). It depends upon the entire history of the process, which is typical for the TF electrodynamics of media with memory but unnecessary for a vacuum [61]. Important is to notice that we cannot distinguish in (113) and (114) power and energy components, respectively, which are conserved in the field and dissipated in the medium described by FOM.
For a vacuum when β = 1, (114) gives the standard formula for the energy density in classical IO electrodynamics As one can notice, the energy density formula (115) is local in time for β = 1. The case of β ∈ (0, 1) is totally different. The formula (114) is nonlocal in time not only because of the integral I 1 t , but also due to nonlocality of fractional derivatives D β t . Let us demonstrate that nonlocal nature of the derivative D β t introduces appropriate memory effects into the energy-density formula (114). Consider now a local in time RS vector as an 'input signal' to (114), which approximates Dirac's delta distribution where v is a fixed vector such as v · v * = |v| 2 = 1 and T denotes an approximation parameter. One should note that although F T is not a solution to the diffusion-wave Equation (84), it allows to analyse (114) in terms of memory effects. Let us show that F T influences the right-hand side of (114) not only for t ∈ [0, T] but for all times t > 0, when β ∈ (0, 1). It means that memory effects can be observed in this case for t > T. One can see that and It shows that the effect of 'input signal' F T with compact support is accumulated in the density of electromagnetic-field energy up to t → +∞. Memory effects, for T = 0.1 s, are presented in Figure 5.
As one can see, for β = 1, the energy density u is nonzero only for t ∈ [0, T], when the RS vector F T is nonzero as well. However, for β ∈ (0, 1), the energy density u saturates when nonzero F T appears. The amplitude of saturation decreases when β is decreased. One should note that for β = 1 and the discontinuous function F T , we may not use the formula (115), because we may not apply the Leibniz product rule to the non-differentiable functions. Hence, we need to directly apply the formula (114) in a distributional sense.
Suppose now that V ⊂ R 3 is a compact volume with the boundary S being the piecewise smooth surface. Hence, with the use of Gauss's theorem and (44), the theorem (111) can be written in the integral form as Equation (120) describes the balance of the power dissipated and stored in the electromagnetic field in the considered volume V. Since S ∈ R 3 , then one obtains Equation (122) is always valid for S ∈ R 3 , because Hence, the imaginary part of expression under the integral (122) is identically equal to zero. In (121), the volume integral involves dot products between E/H fields and their FO derivatives. On the other hand, one can notice in (122) that the volume integral involves dot products between E/H fields and FO derivatives of H/E fields, respectively. It can be treated, to some extent, as the orthogonality condition between the components of the RS vector and their FO derivatives. With the use of (114), one can calculate the total electromagnetic-field energy (i.e., Hamiltonian of the field and the medium) as follows: (124)

Momentum Conservation
Let us formulate the momentum-conservation theorem [48,49] in TF electrodynamics. Let us assume that M = 0 and ρ m = 0 in (23)- (26). Then, one can formulate the momentum-balance equation without using the definition of Lorentz force and the constitutive relations [62,63]: Lack of symbol between vectors denotes a tensor or dyadic product. The right-hand side of (125) represents the force density which includes, besides the Lorentz force with total fields, other force densities, e.g., Helmholtz's for gases [62]. Above equation can be written using the RS vector as follows: In (126), T denotes the stress tensor, i.e., and f denotes the force density, i.e., For β = 1, (128) reduces to the standard formula for the Lorentz force, i.e., H, then additional contribution to the force density is obtained due to the application of FOM.
Suppose now that V ⊂ R 3 is a compact volume with the boundary S being the piecewise smooth surface. Then, using Gauss's theorem, the formula (126) can be written in the integral form asˆV Equation (130) describes the momentum conservation in the electromagnetic field in the considered volume V, i.e., the total change of the mechanical (i.e., related with an electromagnetic medium) and electromagnetic-field momentum is equal to the momentum flowing into and out of the volume V. Finally, one can calculate the total electromagnetic-field momentum G M and the total angular momentum L M as follows: Although the stress tensor (127) is asymmetric, one can formulate the balance equation of angular momentum using the approach proposed in [64].

Uniqueness of Solutions
Let us consider the boundary S of the volume V ⊂ R 3 , which is the piecewise smooth surface. Then, let us define either the tangential electric-or magnetic-field intensity. Now, suppose that two solutions F 1 and F 2 of TF Maxwell's equations exist in V. Hence, the solution F d = F 1 − F 2 of TF Maxwell's Equations (40) and (41) is obtained assuming that either the tangential electric or magnetic field is equal to zero at the boundary S. The function F d satisfies homogeneous Equations (40) and (41) By applying the Fourier transformation to both sides of (133), one obtains whereF d (r, ω) = F (F d (r, t)). Taking the complex conjugates, one obtains When the operators (F * d ·) and (F d ·) are respectively applied to (135) and (136), and then these equations are subtracted, one obtains Now, using Gauss's theorem, one obtains Since ((jω) β + (−jω) β ) = 2|ω| β cos(β π 2 sgn(ω)), there is Hence, for β ∈ (0, 1), there is´V |F d | 2 dv = 0 and consequentlyF d = 0. Finally, one obtains thatF 1 =F 2 . As one can expect, it means that also for TF Maxwell's equations in the RS representation, the solutions are unique for assumed boundary conditions.

Reciprocity
Reciprocity describes an electromagnetic system, whose a response to a source is unchanged when the source and the measurer are interchanged [65]. We extend this property hereby towards the FO electrodynamics, defined based on the RS vector, by deriving the Lorentz reciprocity theorem [54].
Let us consider a volume V of the medium described by FOM, which includes sources K. Then, taking the Fourier transformation with respect to the variable t, curl Maxwell's Equation (40) can be written as whereF(r, ω) = F (F(r, t)) andK(r, ω) = F (K(r, t)). Let us consider two solutions of FO Maxwell's equationsF i existing for current sourcesK i in the volume V (i = 1, 2). Then, one obtains Taking into consideration that Equation (144) is a differential form of the Lorentz reciprocity theorem in TF electrodynamics, formulated making use of the RS vector. Applying Gauss's theorem to (144), one obtains the integral form of the Lorentz reciprocity theorem Let us assume that the surface integral in (145) vanishes. It is a realistic assumption, because the surface of integration can be extended to infinity from the sources. Then, for β = 1, the electric and magnetic fields are related by the plane-wave formulas Z fH1/2 = n ×Ẽ 1/2 and n ·Ẽ 1/2 = 0, where n is the unit vector normal to the integration surface S pointing outwards. However, for 0 < β < 1, the electromagnetic field is exponentially attenuated towards zero, when the surface of integration is extended to infinity from the sources. Hence, one obtainsˆV This means that volume integrals of the dot products of the electromagnetic fieldsF 1/2 (resulting from the excitationsK 1/2 and measured in their positions) and the current sourcesK 2/1 are the same.

Causality
Now let us quickly review the causality. The complex-valued transfer function/distribution in the frequency domain G : R → C is considered, which is the Fourier transform of a certain time-domain function/distribution g : R → C. One should note that we do not assume that g(t) is real-valued; hence, we are not be able to use any simplification that results from the assumption that g(t) is real valued.
The classical perspective is provided by the Titchmarsh theorem, which works for the function g ∈ L 2 (R). Theorem 1 (see Theorem 1.6.1 in [66] and Theorem 2 in [55]). If a square-integrable function G(ω) fulfills one of the four conditions below, then it fulfills all four of them: (i) The inverse Fourier transform g(t) of G(ω) vanishes for t < 0: (ii) G(v) is, for almost all v, the limit as u → 0 + of an analytic functionG(u + jv) that is holomorphic in the right half-plane and square integrable over any line parallel to the imaginary axis:ˆ∞ −∞ |G(u + jv)| 2 dv < C (u > 0).
(iii) G and G verify the first Plemelj formula: (iv) G and G verify the second Plemelj formula: This theorem suggests two different methods of proving the causality, i.e., one requires searching for an appropriate holomorphic extension to the right-half plane, and the other requires to prove the validity of the Kramers-Krönig (K-K) relations (147) and (148).
One should note that for both functions G z and G R given by (96) and (102), respectively, the appropriate holomorphic extensions exist for β ∈ (0, 1), given by G z (σ + jω) = G z (s) =F + e −zc −1 β s β (149) Since all the assumptions given in point (ii) of the Titchmarsh Theorem 1 are satisfied, it directly proves the causality.
On the other hand, the case of transfer function with the cylindrical symmetry is different. The function G r , given in (109), does not belong to L 2 (R) due to the asymptotics of the Bessel function J 0 given by valid for z ∈ C such that arg z ∈ (−π, π) as |z| → +∞ (see formula 9.2.1 in [67], more details are given in Chapter VII in [57]). The function G r (ω) has a natural holomorphic extension to the right half-plane given by G r (σ + jω) = G r (s) =F + J 0 (js β r/c β ).
However, since we operate outside L 2 (R), we may not directly apply classical methods and we should treat the function G r (ω) as a tempered distribution. Fortunately, there are tools that may be applied here as well: we are going to refer to the theorem proved in [68] and rephrased as Theorem 7 and discussed in [55].
The function (152) satisfies (due to the estimate (151)) all the assumptions of Theorem 2, so it is also a causal transform.

Conclusions
In this paper, the extension of TF electrodynamics towards the formulation based on the RS vector is presented. Up to now, the existing formulations of electrodynamics based on the RS vector do not include the energy dissipation; hence, they are limited to a homogeneous and static medium, which is usually a vacuum. The time-domain fractionalization of Maxwell's equations allows for inclusion of energy dissipation, which is the intrinsic property of media described by FOMs. TF Maxwell's equations are formulated using the RS vector and their properties are analysed from the point of view of classical electrodynamics, i.e., energy and momentum conservation, reciprocity, causality. Therefore, it is possible to demonstrate that the energy density of the electromagnetic field and media described by FOM, which is spent to create the field and change the momentum of medium, reflects the entire previous history of the process. Hence, the application of FOM allows for description of memory effects in electrodynamics. Afterwards, classical solutions are derived for wave-propagation problems, assuming helical, spherical and cylindrical symmetries of solutions. Furthermore, we demonstrate relations between the TF Schrödinger equation and TF electrodynamics using the proposed formulation based on the RS vector.