On Fundamental Solution for Autonomous Linear Retarded Functional Differential Equations

: This document focuses attention on the fundamental solution of an autonomous linear retarded functional differential equation (RFDE) along with its supporting cast of actors: kernel matrix, characteristic matrix, resolvent matrix; and the Laplace transform. The fundamental solution is presented in the form of the convolutional powers of the kernel matrix in the manner of a convolutional exponential matrix function. The fundamental solution combined with a solution representation gives an exact expression in explicit form for the solution of an RFDE. Algebraic graph theory is applied to the RFDE in the form of a weighted loop-digraph to illuminate the system structure and system dynamics and to identify the strong and weak components. Examples are provided in the document to elucidate the behavior of the fundamental solution. The paper introduces fundamental solutions of other functional differential equations.


Introduction
This paper examines and characterizes the fundamental matrix solution for the n-vector autonomous linear retarded functional differential equation (RFDE) with forcing function f (t) and initial function h(θ) where x(t) is an n-vector, the kernel matrix A(ds) is an n × n matrix of Borel measures with support on the interval [0, a], and f and g belong to function spaces that will be described in Section 2.
We shall consider the case that where 0 = θ 0 < θ 1 · · · < θ N = a, A i are n × n matrices over R, δ(s) is the Dirac delta function, and A(s) is an integrable n × n matrix function on the interval [0, a]. Consequently, in terms of the Radon-Nikodym theorem, A(ds) has an absolutely continuous part A(s), a discrete singular part with a finite number of point masses (atoms) corresponding to Dirac delta measures, and a zero continuous singular part. It will be evident from the context whether the symbol A refers to an n × n matrix over Borel measures, an n × n matrix over the real numbers, or a matrix function. The terms ∑ N i=0 A i δ(s − θ i ) and A(s)ds respectively correspond to discrete delays and distributed delay of the kernel matrix A(ds). The fundamental solution Φ(t) satisfies the RFDE Φ (t) = a 0 A(ds)Φ(t − s), t > 0, (4) with initial conditions where I and O are respectively the n × n identity and zero matrices. Note that the fundamental solution Φ(t) also satisfies the Volterra integro-differential equation with initial condition Theorem 1. The fundamental solution Φ(t) is given by where the convolutional powers of the kernel matrix A n (ds) are defined by 1. A 0 (ds) = Iδ(s), 2. A 1 (ds) = A(ds), and 3. A n (ds) = A n−1 (ds) * A 1 (ds) for n ≥ 2, where * denotes the convolution of two matrix measures.
Proof. See the proof of Theorem 2.
This expression for the fundamental solution in the case n = 1 can be found in References [1][2][3]. Note that the computation for the fundamental solution Φ(t) is simplified by the restriction of a finite number of discrete atoms in the discrete singular component, and a zero continuous singular component.
The fundamental solution Φ(t) for the autonomous linear ordinary differential equation (ODE) Φ (t) = AΦ(t) is given (see Reference [4]) by the exponential matrix function Note that the fundamental solution Φ(t) in Equation (8) for the RFDE is analogous to and is a generalization of the fundamental solution in Equation (9) for the ODE, and that it has the form of a convolutional exponential matrix function.
Equation (4) establishes a relationship between two actors in an RFDE: the fundamental solution Φ(t) and the kernel matrix A(ds). Let us consider the other actors and their roles.
The fundamental solution Φ(t) can also be expressed in the form of Laplace inverse as a contour integral in the z-plane where contributions to the integral arise from the characteristic roots {λ j } with multiplicity m r of the characteristic determinant D(z) = det(∆(z)) = 0. We shall show in Section 4 that this gives rise to a spectral representation of the fundamental solution Φ(t) in terms of the exponential solutions of the RFDE In the case that all characteristic roots λ r are simple with multiplicity m r = 1, we have where Ψ r = C(λ r )/D (λ r ) with C(z) = adj(∆(z)). The characteristic roots {λ r } can be arranged in some appropriate order, such as decreasing real part (λ r ) or increasing modulus |λ r |. It is known (References [6][7][8][9]) that Φ ci (t) is convergent for a range of values of t, where −a ≤ t < na. The case t = −a arises in the scalar case in which the exponential solutions are complete and independent (see Reference [10]). Even where convergence is not an issue, we could have Φ ci (t) = Φ(t) for t < t 0 and Φ ci (t) ≡ Φ(t) for t ≥ t 0 , as shown in the example in Section 4.4, where t 0 = 2a. An open question is to characterize the minimum time t 0 for the equality of Φ ci (t) and Φ(t) given by t 0 = inf τ>−a {τ; Φ ci (t) ≡ Φ(t) for t > τ}. Note that Φ ci (t) plays the role of a spectral representation of the fundamental solution Φ(t) for t > t 0 . Figure 1 shows the relationship between the RFDE actors: the kernel matrix A(ds), the fundamental solution Φ(t), the characteristic matrix ∆(z), and the resolvent matrix ∆ −1 (z). Note that Figure 1 shows two different routes to determine the fundamental solution Φ(t) from the kernel matrix A(ds): (i) the direct route from Equation (8), and (ii) the indirect route through the characteristic matrix, the resolvent matrix, and the Laplace inverse. Depending on circumstances, either route may be preferable.
On its own, the fundamental solution is a representative proxy for the solutions of the RFDE. For example, the fundamental solution has exponential bounds that determine the asymptotic behavior of the RFDE solutions. Also the fundamental solution combined with a solution representation gives an exact expression in explicit form for the solution of an RFDE. However, the explicit form of the fundamental solution does not in and of itself convey all the information on its properties. It is the diversity of actors and their interplay that brings a variety of perspectives and the full force to bear on extracting information on the RFDE.

Remark 1.
The right hand side of the linear autonomous RFDE (1) is expressed as where A (ds) = A(−ds), to emphasize the convolution nature of the linear autonomous RFDE (see page 77 of Reference [11]), to highlight the convolutional composition of the fundamental solution, and to reflect the situation that for a linear autonomous RFDE essentially every operation is a convolution. (1) is chosen over the alternatives of (i) distributions and (ii) functions of normalized bounded variation. This choice is mainly a matter of style. Whatever choice is made, the fundamental solution will involve convolutional powers of the kernel matrix.
The remainder of the paper is organized as follows: Section 2 on preliminaries establishes the foundation for the paper. Section 3 applies algebraic graph theory to study RFDEs. It uses a weighted loop-digraph representation of an RFDE to illuminate the system structure, connectivity, and dynamics. The characteristics of the fundamental solution are explored in Section 4. Examples are considered in the main body of the document to elucidate the behavior of the fundamental solution to an RFDE. Section 5 extends the fundamental solution to other functional differential equations. Finally, the conclusions are presented in Section 6.

Ring of Borel Measures
Consider the set B + cs of Borel measures on R + with compact support of the form with finite number of atoms at {θ i }, 0 = θ 0 < θ 1 · · · < θ N , a i ∈ R, and integrable function a(s).
The addition operator + is given by and convolution (aka multiplication) operator * is given by It is straightforward to show that B + cs = B + cs (+, * ) is a commutative ring, with additive identity 0(ds) such that (0 + a)(ds) = a(ds), and the multiplicative identity the Dirac delta measure δ(ds) such that δ(s) * a(ds) = a(ds).

Definition 2.
For a(ds) ∈ B + cs define norm |a(ds)| = sup | f (x)|≤1 a(ds) f (s). Define the support supp(a(ds)) as the complement in R + of the set X = {x : N x a(ds) = 0} for all neighborhoods N x of x that are sufficiently small.

Matrix over Ring of Borel Measures
Definition 3. A n × n matrix over a ring of Borel measures M n (B + cs ) is a n × n array of elements of the ring B + cs .
We have matrix addition and matrix convolution (multiplication) Note that the kernel matrix A(ds) ∈ M n (B + cs ), and well as its its convolutional powers A n (ds) ∈ M n (B + cs ).

Fundamental Solution
We now turn our attention to the Volterra integro-differential equation of convolutional type satisfied by the fundamental solution for an RFDE By integrating we obtain the equivalent Volterra integral equation of convolutional type where the kernel matrix function with It is well known that this Volterra integral equation has a Liouville-Neumann series solution in the form (26) where matrix functions B n (t) are defined by B 1 = B, and B n = B n−1 * B 1 where * denotes a convolution of two matrix functions. This form of the fundamental solution is promulgated in References [14,15,20].
Hence from equality of Laplace transforms 1 The solution Φ(t) is unique and has continuous dependence on the kernel matrix A(ds).
so that the series converges absolutely and uniformly, thereby justifying the formal procedures in the following steps: Intermediate steps are To show that the solution is unique, suppose that we have two solutions Φ 1 (t) and Φ 2 (t) so that for the Volterra integral equation we have Consequently for some constant K > 0, arbitrarily small > 0, and To demonstrate continuous dependence on the kernel matrix A(ds), consider two solutions Φ 1 (t) and Φ 2 (t) of We can write so that so that by Gronwall's inequality we have ||Φ We have already seen an exponential bound ||Φ(t)|| ≤ exp(KT) on the interval [0, T]. A more precise exponential bound is provided by the supremum of the real part of the characteristic roots of the characteristic determinant det(∆(z)) = 0.

Remark 3.
The term fundamental solution Φ(t) has several synonyms in the literature depending on context: A. Cauchy matrix C(t, s) for the nonautonomous RFDE See page 51 of Reference [16]. For the autonomous equation C(t, s) = Φ(t − s). B. Differential resolvent for the Volterra integro-differential equation See page 77 of Reference [11].
C. Green's function or impulse response function for RFDE with h = 0 and f (t) = δ(s). In some parts of the literature, for example Reference [16], a Green's function is associated with a boundary value problem.
so that the kernel matrix and the fundamental solution commute convolutionally.
where the steps involving an infinite sum are justified by uniform convergence.

RFDE Solution
Standard choices for the initial function space are: The standard choice for the forcing function space is L 1 loc ([0, ∞), R n ). Results on the existence, uniqueness and continuous dependence on the initial data can be found in Reference [22] for function space F 1 and Reference [13] for function space F 2 . (1) is

Theorem 4. The representation of the solution x(t) for the RFDE in Equation
Proof. We first consider the case h(0) = 0, h = 0 and f = 0. Taking the Laplace transform we obtain and taking the Laplace inverse, we have We have the equivalent equation with h(0) = 0 and h = 0. This equation has the solution From linearity, we can combine the two solutions to obtain the complete solution.

Remark 6. The integration region R for the double integral is
The alternative expressions for the double integral are: The change in the order of integration in the double integral is justified by the Fubini theorem. The change in the order of the kernel matrix and the fundamental solution is justified since they commute convolutionally.

Application of Algebraic Graph Theory
In this section we apply algebraic graph theory to examine and to provide a pictorial rendition of the system structure and system dynamics of an RFDE. We also use the theory to characterize the system connectivity of the RFDE, and to identify its strong and weak components. The main reference for the application of algebraic graph theory to RFDEs is the excellent book [23].

Weighted Loop-Digraph Representation
We consider a weighted loop-digraph representation for the mathematical objects of the kernel matrix A(ds) and its convolutional powers A n (ds). These mathematical objects are prominent aspects of the RFDE and its fundamental solution Φ(t). In where The component form expression for the RFDE is Likewise Figure 2 depicts the network structure and interactive systems dynamics for a general RFDE in the form of a weighted loop-digraph, albeit for simplicity focusing only on two nodes i and j. A(ds):

Definition 5.
A weighted loop-directed graph or loop-digraph is an ordered triple G = (V, E, W) where V is a set of elements called vertices (or nodes, points), and E is a set of ordered pairs of vertices called directed edges and loops. The directed edges connect two distinct vertices in the given order. Loops connect a specific vertex with itself. W is the weight assigned to each directed edge and loop.
Although for the most part the book [23] focuses on directed graphs without loops and without weights, the applied results in this paper are applicable to directed graphs with loops and weights. Figure 2 shows a weighted loop-digraph for an n-vector RFDE. The weights assigned to directed edge e ij and loop e ii are Borel measures B + cs of a ij (ds) and a ii (ds), respectively. By convention the direction on the directed edge e ij is i to j. However, note that in a physical interpretation the influence of node j on the system dynamics at node i is represented by the Borel measure a ij (ds) in the opposite direction. This consideration should be borne in mind in interpreting the system dynamics from the loop-digraph.
For the components of the second convolutional power A 2 (ds) we have This has the interpretation that the Borel measure a 2 ij (ds) is the sum of convolutions of the Borel measures encountered on the directed edges and loops on a walk of length 2 from vertex i to vertex j. More generally the Borel measure a n ij (ds) is the sum of convolutions of the Borel measures encountered on the directed edges and loops on a walk of length n from vertex i to vertex j. The Borel measures a n ij (ds) thus have a straightforward interpretation in the language of algebraic graph theory.
The following examples illustrate the concepts.
The fundamental solution is where The weighted loop-digraphs for the kernel matrix and its convolutional powers for the RFDE are shown in Figure 3.

Example 2.
The RFDE is The kernel matrix and its convolutional powers are A n (ds) = a n 11 δ(s − n) 0 0 a n 22 δ(s − n) .
The fundamental solution is where The weighted loop-digraphs for the kernel matrix and its convolutional powers for the RFDE are shown in Figure 4.

Strong and Weak Connectivity
Notice that in the first example from the previous subsection that the two components in the RFDE are interconnected and interact, and that there is no interaction in the second example. The two components are interconnected in the first example, and disconnected in the second example. The nature of the interaction between the RFDE components will be obvious in the case that the number of components is small, but less so when the number is large.
In this subsection we apply material from Chapters 2,3,5,14 of Reference [23] to characterize the connectivity between RFDE components and to classify the nature of the connectivity. Notice that this approach focuses on the structure of the RDFE, and is not concerned with the dynamics produced by the Borel measures, i.e., the weights on the directed edges and loops.
We start with some definitions.

Definition 6.
A path from v 1 to v n is an ordered collection of vertices v 1 , v 2 · · · v n interspersed with an ordered collection of directed edges v 1 v 2 , v 2 v 3 , · · · , v n−1 v n . A semipath from v 1 to v n is an ordered collection of vertices v 1 , v 2 · · · v n interspersed with an ordered collection of n − 1 directed edges, one from each pair of directed edges v 1 v 2 or v 2 v 1 , v 2 v 3 or v 3 v 2 , · · · , v n−1 v n or v n v n−1 .
A strict semipath is a semipath that is not a path. A vertex v is reachable from a vertex u if there is a path from u to v. Points u and v are 0-connected if they are not joined by a semipath; 1-connected if they are joined by a semipath but not a path; 2-connected if they are joined by a path in one direction but not the other; 3-connected if they are joined by paths in both directions.
A loop-digraph is strongly connected or strong if every two vertices are mutually reachable. A loop-digraph is unilaterally connected or unilateral if for any two points, at least one is reachable from the other; it is strictly unilateral if it is unilateral but not strong.
A loop-digraph is weakly connected or weak if every two vertices are joined by a semipath; it is strictly weak if it is weak but not unilateral.
A subgraph of a loop-digraph D is a loop-digraph in which the vertices and loops/directed edges are vertices and loops/directed edges of D.
A strong component of a loop-digraph is a maximal strong loop-digraph. A unilateral component of a loop-digraph is a maximal unilateral loop-digraph. A weak component of a loop-digraph is a maximal weak loop-digraph.

Theorem 5 (Theorem 3.2 of Reference [23]). Every vertex and every directed edge is contained in exactly one weak component. Every vertex and every directed edge is contained in at least one unilateral component. Every vertex and every directed edge is contained in exactly one strong component.
We shall now proceed to demonstrate how to decompose a loop-digraph into strong and weak components.
The reachability matrix R = [r ij ] is defined as In terms of the convolutional powers A k (ds) for 0 ≤ k ≤ n − 1 of the kernel matrix A(ds), the reachability matrix R = [r ij ] is given by Other variants of the reachability matrix come into play: The transpose reachability matrix R T = [r T ij ] is given by The symmetric reachability matrix R S = [r S ij ] arises from the symmetrized kernel matrix A S (ds) where a S ij (ds) = a ij (ds) + a ji (ds). It is given by The element-wise product R × R T = [r P ij ] = [r ij r ji ]. It is given by r P ij = 1, vertex i is reachable from vertex j and vertex j is reachable from vertex i, The following theorem says that the strong components of the RFDE are determined by the elementwise product reachability matrix R × R T . Theorem 6 (Theorem 5.8 of Reference [23]). The strong component containing the vertex v i is given by the entries of 1 in the ith row (or column) in the elementwise product reachability matrix R × R T .
The following theorem says that the weak components of the RFDE are determined by the symmetric reachability matrix R S . Theorem 7 (Corollary 5.15a of Reference [23]). The weak component containing the vertex v i is given by the entries of 1 in the ith row (or column) in the symmetric reachability matrix R S .
The connectedness matrix C = [c ij ] is obtained from the reachability matrix R = [r ij ] as follows: If v i and v j are in the same weak component c ij = r ij + r ji + 1.
The connectedness matrix C = [c ij ] takes the values 0, vertices i and j are not connected with a semipath, i.e., they are disconnected, 1, vertices i and j are connected with a strict semipath, 2, vertices i and j are connected with a path in one direction but not the other, 3, vertices i and j are connected with paths in both directions.
(89) Definition 7. A permutation matrix P ∈ M n (B + cs ) has exactly one entry in each row and in each column set to δ(ds). All other entries are 0. A permutation matrix P can be construed as relabelling the vertices in a loop-digraph.
A n × n matrix M ∈ M n (B + cs ) is said to be decomposable if ∃ a permutation matrix P ∈ M n (B + cs ) such that where M 11 and M 22 are square matrices of sizes r × r and n − r × n − r respectively, and O r×n−r and O n−r×r are zero matrices of the indicated sizes. Then the fundamental solution Φ(t) is also decomposable in the form where Φ 11 (t) and Φ 22 (t) are the fundamental solutions corresponding respectively to the kernel matrices A 11 (ds) and A 22 (ds).
In general, if the loop-digraph for an RFDE has p weak components, the kernel matrix A(ds) can be decomposed into the form and the fundamental solution can also be decomposed into a corresponding form where Φ ii (t) is the fundamental solution for the kernel matrix A ii (ds).

Example 3.
To illustrate the results, we consider the RFDE weighted loop-digraph shown in Figure 5. A(ds):

The kernel matrix A(ds) is
where the dummy variable (ds) is suppressed within the matrix to reduce clutter. The symmetric kernel matrix A S (ds) is The reachability matrix R is The transpose reachability matrix R T is The symmetric reachability matrix R S is This yields the weak components {v 1 , v 2 , v 3 , v 6 } and {v 4 , v 5 }. The element-wise product reachability matrix R × R T is The connectedness matrix C is given by By applying the appropriate permutation matrix or simply interchanging the rows of columns of the kernel matrix A(ds) in accordance with the weak components, the kernel matrix is expressed in decomposed form, as follows: The kernel matrix A(ds) is The kernel matrix A(ds) is expressed succinctly in the decomposed form where A 11 (ds) and A 22 (ds) are respectively 4 × 4 and 2 × 2 kernel matrices, and O 2×4 and O 4×2 are respectively 2 × 4 and 4 × 2 zero matrices. The fundamental solution Φ(t) can be expressed in decomposed form where Φ 11 (t) and Φ 22 (t) are respectively the fundamental solutions for the kernel matrices A 11 (ds) and A 22 (ds).

Fundamental Solution is Laplace Inverse of Resolvent Matrix
The fundamental solution satisfies the Volterra integro-differential equation of convolutional type Taking the Laplace transform, we have Consequently, the fundamental solution is the Laplace inverse of the resolvent matrix ∆ −1 (z). Alternatively, the Laplace transform of the fundamental solution Φ(t) is the resolvent matrix ∆ −1 (z).

Tauberian Behavior
From the previous subsection, we havê Assume that all the characteristic roots of the characteristic determinant det(∆(z)) have real part less than zero, so that α 0 < 0. Choose α so that α 0 < α. Then Φ(t) exp(−αt) is integrable and we have (112)

Contour Integral Version of Fundamental Solution
Let us consider the expression for the fundamental solution Φ ci (t) in the form of Laplace inverse as a contour integral in the z-plane For notational convenience, we put C(z) = adj(∆(z)), D(z) = det(∆(z)), and E(z) = 1/ det(∆(z)).
To obtain the contribution to the integral of a characteristic root λ r of det(∆(λ)) = 0 with multiplicity m r , we consider the integral over a small circle centered at λ r containing no other characteristic roots of det(∆(λ)) = 0.
Following the approach in Reference [10], we have We have that E(z), the reciprocal function of D(z), is given by where E 0 (λ r ) = m r !/D ((m r )) (λ r ), and E n (λ r ) for n ≥ 1 is obtained recursively from We have The residue for the characteristic root λ r is the coefficient of (z − λ r ) −1 corresponding to k = k 1 + k 2 + k 3 = m r − 1. It is given by The matrix coefficient Ψ rj is given by For example Ψ r,m r −1 = C(λ r )m r /D (m r ) (λ r ), and for m r = 1 Ψ r,0 = C(λ r )/D (λ r ).
We have arrived at the following theorem: Theorem 10. Let {λ r } be the characteristic roots of the characteristic determinant det(∆(z)) = 0 with multiplicity m r . The version of the fundamental solution Φ ci (t) arising from a contour integral has the form of spectral decomposition of exponential solutions where In the case that all characteristic roots λ r are simple with multiplicity m r = 1 we have where Ψ r = C(λ r )/D (λ r ). For the scalar case n = 1, C(z) ≡ 1 so that Ψ r = 1/D (λ r ).

Example 4.
Consider x 3 (t) = 0 . (127) A n (ds) = 0 for n ≥ 3 so that A(ds) is a nilpotent matrix with index 3. The weighted loop-digraphs for the kernel matrix and its convolutional powers for the RFDE are shown in Figure 6.
The fundamental matrix Φ(t) is The contribution to the contour integral is The contour integral version of the fundamental solution Φ ci (t) is given by Note for this specific example that Φ ci (t) ≡ Φ(t) for t ≥ 2a, but that Φ ci (t) = Φ(t) for 0 ≤ t < 2a. As this example has only a finite number of terms in Φ ci (t) and Φ(t), convergence is not a factor in the discrepancy. The issue is that Φ(t) is compelled to contain truncated powers of t to reflect delay activations, whereas by its nature Φ ci (t) is compelled to contain untruncated powers of t.
An open question is to characterize the minimum time t 0 for the equality of Φ ci (t) and Φ(t) given by t 0 = inf τ>−a {τ; Φ ci (t) ≡ Φ(t) for t > τ}. The matter has been encountered and considered in the literature (see References [6][7][8][9]) with values of t 0 ranging from −a to na.

Semigroup Relations
The fundamental solution Φ(t) for the ODE x (t) = Ax(t) satisfies the semigroup relation Φ(t 1 + t 2 ) = Φ(t 1 )Φ(t 2 ). The generalization for the fundamental solution Φ(t) for an RFDE is as follows: Theorem 11. For t 1 0 and t 2 0 Proof. This follows from the representation of the solution with f = 0 and the consideration that the translate of a solution of an autonomous linear RFDE is also a solution.

Remark 7.
A related result is found in Reference [24].

Remark 8.
It follows from the representation of solutions that the double integral has the following equivalent forms (see Equation (58)):

Semigroup of Solution Operators
In this subsection we consider the role played by the fundamental solution Φ(t) in solution operators mapping the initial function h to the solution x(t; h) in the function spaces M 2 ([−a, 0], R n ) and C([−a, 0], R n ).
The solution x(t; h) is given by Remark 9. Note that a solution operator has to accommodate two aspects of the initial function h: (i) the initial value h(0) acted on by the fundamental solution Φ(t), and (ii) the initial function h as an integrand with kernelΦ(t, α).

Picard Iteration
Definition 8. The n th Picard iteration {Φ n (t)} for the fundamental solution Φ(t) based on its Volterra integral Theorem 12. The n th Picard iteration Φ n (t) is given by

p-th Order RFDE
We consider fundamental solution Φ(t) for the p-th order RFDE with n-vector x(t) By the standard procedure we convert the p-th order RFDE for n-vector x(t) into a first order RFDE for a pn-vector y(t) where the pn × pn companion matrix A(ds) is given by · · · 0 0 0 I n δ(s) · · · 0 · · · · · · · · · · · · · · · A(ds) The pn × pn fundamental solution is where In particular, for p = 2 where See Reference [28].

Infinite Delay
The RDFE with finite delay has a kernel matrix A(ds) ∈ B + cs with compact support and an initial function space defined over a finite interval [−a, 0]. In contrast, the RDFE with infinite (or unbounded) delay has a kernel matrix with non-compact support and an initial function space defined on the infinite interval (−∞, 0]. There is a wide variety of choices for this initial function space and the theory for the RDFE with infinite delay is intricately intertwined with a choice (see Reference [29]). However, the fundamental solution, with its simple initial condition of Φ(0) = I and Φ(t) = O for t < 0 is exempt from these intricacies.
In this expression we assume either (i) we are dealing with the scalar case or (ii) that the matrices A 0 · · · A N commute. Note that for fixed t we have a finite sum if A 0 = 0. The method of steps for computing the fundamental solution Φ(t) will work in this case. (169) is

The fundamental solutions for Equation
It is possible, as demonstrated in the following example, that in the case of an infinite delay that the fundamental solution satisfies an ODE-without a delay.

Example 5.
Consider the scalar case with A(t) = exp(−Kt) ∑ p i=0 a i t i so that Multiplying Equation (172) by exp(Kt), differentiating p + 1 times, and then multiplying by exp(−Kt), we obtain that the fundamental solution Φ(t) satisfies the (p + 2) th order ODE Consequently the fundamental solution satisfies an ODE with no delay.

Nonautonomous RFDE
We consider the fundamental solution Φ(t, s) for the nonautonomous RFDE where A(t, τ) is measurable in t, of bounded variation in τ for s < τ < t, constant for τ ≥ t, and T s dt T s d τ A(t, τ) < ∞ .
The fundamental solution Φ(t, s) (also known as the Cauchy function C(t, s)) satisfies initial conditions Φ(s, s) = I, and Φ(t, s) = O for t < s.
We integrate to obtain the equivalent Volterra integral equation (178)

Conclusions
This paper has extended the expression for and the treatment of the fundamental solution for autonomous linear RFDE from scalar to n-vector, in the following ways: The fundamental solution is presented in the form of a convolutional exponential matrix function. The paper demonstrates how RFDE analysis can be conducted through an interplay of the RFDE actors (fundamental solution, kernel matrix, characteristic matrix, characteristic determinant, resolvent matrix) and tools (Laplace transform and solution representation). Elements of algebraic graph theory are applied in the form of a weighted loop-digraph to illuminate the system structure and dynamics, and to identify the strong and weak components. This is potentially a powerful tool in the case that the RFDE has a special structure, such as coupled oscillators with delay interactions. The paper compares the fundamental solution Φ(t) of an RFDE with its contour integral version Φ ci (t) for a special example, and raises an open question on the minimum time t 0 for equivalence. For t > t 0 , Φ ci (t) plays the role of a spectral representation for the fundamental solution Φ(t). In its explicit form and in tandem with the solution representation, the paper describes the key role played by the fundamental solution in the semigroup of solution operators. The paper considers the fundamental solution for other functional differential equations (p-th order RFDE, RFDE with infinite delay, nonautonomous RFDE), but stops short of Neutral and Partial functional differential equations. Finally, the paper opens the door to future work on fundamental solutions and associated analysis for (i) RFDEs with a special structure and (ii) Neutral and Partial functional differential equations.