Discriminant and Root Trajectories of Characteristic Equation of Fractional Vibration Equation and Their Effects on Solution Components

: The impulsive response of the fractional vibration equation z (cid:48)(cid:48) ( t ) + bD α t z ( t ) + cz ( t ) = F ( t ) , b > 0, c > 0,0 ≤ α ≤ 2, is investigated by using the complex path-integral formula of the inverse Laplace transform. Similar to the integer-order case, the roots of the characteristic equation s 2 + bs α + c = 0 must be considered. It is proved that for any b > 0, c > 0 and α ∈ ( 0,1 ) ∪ ( 1,2 ) , the characteristic equation always has a pair of conjugated simple complex roots with a negative real part on the principal Riemann surface. Particular attention is paid to the problem as to how the couple conjugated complex roots approach the two roots of the integer case α = 1, especially to the two different real roots in the case of b 2 − 4 c > 0. On the upper-half complex plane, the root s ( α ) is investigated as a function of order α and with parameters b and c , and so are the argument θ ( α ) , modulus r ( α ) , real part λ ( α ) and imaginary part ω ( α ) of the root s ( α ) . For the three cases of the discriminant b 2 − 4 c : > 0, = 0 and < 0, variations of the argument and modulus of the roots according to α are clariﬁed, and the trajectories of the roots are simulated. For the case of b 2 − 4 c < 0, the trajectories of the roots are further clariﬁed according to the change rates of the argument, real part and imaginary part of root s ( α ) at α = 1. The solution components, i.e., the residue contribution and the Hankel integral contribution to the impulsive response, are distinguished for the three cases of the discriminant.

Let f (t) be piecewise continuous on (0, +∞) and integrable on any finite subinterval of (0, +∞). Then, the Riemann-Liouville fractional integral of f (t) of order β is defined as for β > 0, and I 0 t f (t) = f (t) for β = 0, where Γ(·) is the gamma function. Let α be a positive real number, m − 1 < α < m, m ∈ N + . Then, the Riemann-Liouville and Caputo fractional derivatives of f (t) of order α are defined as Theoretical analyses and experimental simulations indicated that the stress-strain relationship of a viscoelastic body can be better characterized by using fractional constitutive equations, such as the Scott-Blair model [4,9], the Kelvin-Voigt, Maxwell and Zener models [10,[13][14][15] and others [16][17][18]. Thus, dynamics involving a viscoelastic medium with a fractional stress-strain relationship lead to fractional vibration equations or fractional oscillator equations. Different types of fractional vibration equations were presented and analyzed, such as in [7,8,[19][20][21][22][23][24][25][26]. We mention the more general forms known as the multi-term fractional Bessel equations [27] and the multi-term fractional quasi-Bessel equations [28], where the existence theory of solutions was constructed in the class of fractional series solutions.
We consider the fractional vibration equation in the form where the zero initial values z(0) = z (0) = 0 are equipped with, when α is the integer 0, 1 or 2, D α t z(t) denotes the integer-order derivative z (α) (t), otherwise D α t denotes the Riemann-Liouville or Caputo fractional derivative operator. The impulsive response to the Dirac delta driving function F(t) = δ(t) is expressed by the inverse Laplace transform as z(t; α) = L −1 1 s 2 + bs α + c . (2) We note that under the zero initial conditions z(0) = z (0) = 0, the Riemann-Liouville and Caputo fractional derivatives of order α are consistent [2]. Otherwise for nonzero initial values, which the fractional derivative must be indicated. The impulsive response is also known as the Green's function or the fundamental solution [2, 23], from which the solution of Equation (1) with nonzero initial values can be expressed through integration and convolution.
By the complex path-integral formula of the inverse Laplace transform, Equation (2) has the form z(t; α) = 1 2πi Br e st s 2 + bs α + c ds, where Br denotes the Bromwich path, i.e., a straight line on the s plane from s = γ − i∞ to s = γ + i∞ such that the singularities of the integrand sit on its left side. For a further calculation by the residue theorem, we need to find out the zeros of the characteristic equation where, when α is not integers, the one-valued branch s α = e α(ln |s|+i arg s) is taken and the zeros are limited on the principal Riemann surface −π < arg s ≤ π. Notice that the characteristic Equation (4) can be directly written out from Equation (1) like an integerorder case. We note that the characteristic equation in this paper is introduced through the inverse Laplace transform and is different from the characteristic equation Dubovski and Slepoi [27,28] introduced to determine the parameter in a series solution. Moreover, we note that the homogeneous fractional vibration equation belongs to the types investigated in [27,28], where the method of series solutions was proposed.
In [21], Naber considered Equation (4) in the case of 0 < α < 1 and proved that for any specified b > 0, c > 0 and α ∈ (0, 1), Equation (4) has a pair of conjugated complex roots with a negative real part on the principal Riemann surface −π < arg s ≤ π. In [22], Wang and Du investigated the case 1 < α < 2 and proved that the conclusion is still true.
It is known that when α = 1, the roots of Equation (4) have three different forms clarified by the three cases of the discriminant b 2 − 4c: > 0, = 0 and < 0, which correspond to the over damping, critical damping and under damping for the classical vibration equation z (t) + bz (t) + cz(t) = F(t). How the three cases of the discriminant b 2 − 4c affect the couple conjugated complex roots in the circumstance of α = 1, and further on the solution of the fractional vibration equation, has not been reported to our knowledge.
In particular, for the case of b 2 − 4c > 0, when α approaches 1, how the couple conjugated complex roots approach the two different real roots of the case α = 1 is not an apparent status. This is the motivation for this work.
In this paper, we consider the roots of the characteristic Equation (4) and the impulsive response of the fractional vibration Equation (1) in the three cases of the discriminant b 2 − 4c, where the range of the order α covers the integers 0, 1 and 2, and also the noninteger cases 0 < α < 1 and 1 < α < 2, simultaneously. The text is organized as follows. In Section 2, we prove that for any b > 0, c > 0 and α ∈ (0, 1) ∪ (1, 2), the characteristic Equation (4) has exactly a pair of conjugated simple complex roots with a negative real part on the principal Riemann surface. For the three cases of the discriminant b 2 − 4c, the variations of the argument and modulus of the roots according to α are clarified. In Section 3, the trajectories of the roots s(α) on the upper-half complex plane are analyzed and simulated for the three cases of the discriminant b 2 − 4c. In Section 4, a particular discussion for the case of b 2 − 4c < 0 is conducted, where the trajectories of the roots s(α) are further clarified according to the change rates of the argument θ(α), real part λ(α) and imaginary part ω(α) of the root s(α) at α = 1. In Section 5, for the three cases of the discriminant b 2 − 4c, the residue contribution and Hankel integral contribution to the solution of the fractional vibration equation resulted from the inverse Laplace transform are studied. Section 6 summarizes our conclusions.

Argument and Modulus of Roots of the Characteristic Equation
The considered characteristic Equation (4) with three parameters, two coefficients b, c and one power exponent α is a transcendental equation when α is not integers. First, we list results for cases of α being integers. If α = 0 or α = 2, then Equation (4) has a pair of conjugated pure imaginary roots, or Here, the values of α are appended in parentheses to emphasize the dependence of the roots on α, and by putting a bar, we denote the conjugate. If α = 1, then the roots of Equation (4) have the following three forms clarified by the discriminant, The absolute values of the imaginary parts in Equations (5), (6) and (9) have the following relationship. Proposition 1. The imaginary part in Equation (9) satisfies the property: Proof. If −b < b 2 − 4c < 0, then we have 4c − b − b 2 < 0, and further (4c − b 2 )(1 + b) < 4c, from which it follows that 0 < 4c−b 2 4 < c 1+b . Thus, the inequations in (10) hold.
is obtained. In addition, it is apparent that The proof is completed.
We pay particular attention to the case of α = 1, and for this case, by the residue theorem, Equation (3) can be expressed as z(t; 1) = Sum of residues of e st s 2 + bs + c .
For the case of b 2 − 4c > 0, calculating the residues at the two simple poles s 1,2 (1) in Equation (7), we obtain For the case of b 2 − 4c = 0, calculating the residue at the second-order pole s(1) in Equation (8) leads to For the case of b 2 − 4c < 0, calculating the residues at the two conjugated simple poles s(1) ands(1) in Equation (9) yields For the noninteger case α ∈ (0, 1) ∪ (1, 2), s α is multivalued. Here, we take the one-valued branch s α = e α(ln |s|+i arg s) and look for roots of Equation (4) on the principal Riemann surface −π < arg s ≤ π. Because s satisfies Equation (4) if and only if its conjugatē s satisfies Equation (4), we only need to discuss the problem on the half complex plane 0 ≤ arg s ≤ π.
Let the exponential form of root be s = re iθ , where r > 0 and 0 ≤ θ ≤ π. Substituting it into Equation (4) leads to r 2 e i2θ + br α e iαθ + c = 0. Separating the real part and the imaginary part, we have the two equations From the two equations, we determine the argument θ and modulus r of root on the upper-half complex plane for any specified b > 0, c > 0 and α ∈ (0, 1) ∪ (1, 2). It follows from Equation (15) that θ = 0. From Equation (16), θ = π and θ cannot belong to the interval (0, π 2 ]. Thus, the argument of the root s is confined to the interval π 2 < θ < π, and the root, if any, has a negative real part. From Equation (16), we have where the denominator is positive. In order to guarantee sin (αθ) > 0, we must require Thus, we can express the modulus as By inserting Equation (18) into Equation (15) and applying trigonometric formulas, we derive the equation only involves θ free from r as Here, the requirement sin((2 − α)θ) > 0 is necessary, which implies the limitation for the argument The limitations for the argument in Equations (17) and (20) can be combined together as This range of limiting the argument of root is shown in Figure 1.  Rewrite Equation (19) into the following form and define the left-hand side as the function f (θ, α), Our aim is to show the argument θ of the root can be uniquely determined by Equation (22) for any specified b > 0, c > 0 and α ∈ (0, 1) ∪ (1, 2). In Figure 2, sur- For fixed α and b, the function f (θ, α) satisfies the following limits, as shown in Figure 2. Furthermore, the partial derivative with respect to θ is can be estimated as g(θ, α) > 4 sin 2 (αθ) + α 2 sin 2 (2θ) + 4α sin(αθ) sin(2θ) = (2 sin(αθ) + α sin(2θ)) 2 ≥ 0.
Because sin(2θ) < 0, so we have f θ (θ, α) < 0. Thus, by the theorem of implicit function, for any b > 0 and c > 0, Equation (22) uniquely determines an implicit function of the argument of the root vs. α, θ = θ(α), α ∈ (0, 1) ∪ (1, 2), which has the continuous derivative θ (α). With the argument, the modulus of the root is determined from (18) as We summarize the above deduction as follows, taking into account the conjugated part in the lower semi-complex plane.
The graphs of the implicit function θ = θ(α) determined by Equation (22) are shown in Figure 3 for b = 4 and for c = 1, 2, 4, 8 and 16. The outermost thin red dash curve is the boundary of θ, which corresponds to the limit c → 0 + . The six curves arrayed outside-in from the outermost in Figure   In Figure 3, dash line and dot-dash line are for the case of b 2 − 4c > 0, solid line is for the case of b 2 − 4c = 0 and dot line and dot-dot-dash line are for the case of b 2 − 4c < 0.
As c → 0 + , the curve of θ(α) approaches the boundary (24). We note that the left-and right-hand derivatives of the boundary (24) at α = 1 are The function r = r(α) in Equation (23) is composed through the implicit function θ = θ(α). We note that the argument and the modulus of the root in the upper-half plane, θ(α) and r(α), can be extended to define the closed interval 0 ≤ α ≤ 2 by using the results of the integer cases in (5)- (9). That is, for α = 0 and 2, we have It is worth noting that for the case α = 1, b 2 − 4c > 0, the moduli of roots are doublevalued with r 1 (1) > r 2 (1).
We depict the function r = r(α) on the discrete values of α with the step size 0.02: The data plots are shown in Figure 4 for the cases of b 2 − 4c > 0 and b 2 − 4c = 0 and in Figure 5 for the case of b 2 − 4c < 0, where data for the integer cases, α = 0, 1 and 2, are depicted by red '+'. We note that the circular dot line and square dot line in Figure 4 belong to the case b 2 − 4c > 0, and jump discontinuities are shown at α = 1.  and c = 4 (red rhombic dots). For the two implicit functions θ(α) and r(α), the derivatives can be given as follows. First, from (22) follows the partial derivative, By the formula of derivative of implicit function, we obtain and furthermore, from Equation (23), we derive Note that in Equations (30) and (31), θ denotes the implicit function θ(α). The two derivatives will be used in the next section. Proposition 3. Suppose 0 < α < 1, then s = re iθ is a root of equation s 2 + bs α + c = 0 if and only if s = 1 r e iθ is a root of equation Proof. If s = re iθ is a root of equation s 2 + bs α + c = 0, then r 2 e i2θ + br α e iαθ + c = 0, which is equivalent to its conjugate r 2 e −i2θ + br α e −iαθ + c = 0. Then, multiplied by 1 cr 2 e i2θ on the two sides, the equation becomes This means s = 1 r e iθ is a root of Equation (32). The reverse is also true. The proof is completed.
From Figure 3, the curve of θ(α) is not symmetrical about the straight line α = 1 in general. However, from Proposition 3, a sufficient condition for symmetry can be given as follows.

Root Trajectories in Three Cases
In this section, we consider the variation of the root s(α) = r(α)e iθ(α) with respect to α in three cases clarified by the discriminant b 2 − 4c. Our discussion concentrates upon the upper-half complex plane and is based on the fact that the argument of the root, θ(α), is continuous on the interval 0 ≤ α ≤ 2. Special attention is paid to the variation of the root s(α) near α = 1. Case 1. b 2 − 4c > 0 In this case, the limit of the argument of the root is lim α→1 θ(α) = θ(1) = π, which is the peak value of θ(α). For the modulus in Equation (23), by L'Hospital rule Similarly, Due to that θ(1) is the maximum value of θ(α), so θ − (1) ≥ 0 and θ + (1) ≤ 0. Thus, from Equations (33) and (34), we have lim α→1 − r(α) ≥ lim α→1 + r(α). The two limits should equal the absolute values of the two unequal real roots of the case α = 1 in Equation (7), i.e., the following equations hold, Thus, r(α) has a downward jump when α increasingly passes through 1, as shown in Figure 4. Substituting the two limits into Equations (33) and (34), we obtain The limit case of c → 0 + is consistent with the boundary in Equation (25). Therefore, for the root s(α) = r(α)e iθ(α) , we have the following result. Proposition 5. If b 2 − 4c > 0, then in the upper-half complex plane, the root s(α) as a function of α is discontinuous at α = 1, and Taking b = 4, c = 3 and α as the values in Equation (29), the trajectory of roots s(α) in the upper-half complex plane is shown in Figure 6, where s(0) = √ 7 i, s(2) = √ 0.6 i, s 1 (1) = −3 and s 2 (1) = −1, indicated in red crosses, are the roots for the integer cases α = 0, 2 and 1 in the upper-half complex plane.
Re(s) Im(s) i  The root s(α) is continuous at α = 1, and also differentiable, which will be further discussed in the next section.   From the trajectories of roots, s(α) has a faster change rate as α approaches 1. For the cases of b 2 − 4c ≥ 0, the imaginary part of the root s(α) can become arbitrarily small as α approaches 1. For the case of b 2 − 4c < 0, the smaller the value b 2 − 4c, the larger the imaginary part of the root s(1).

Further Discussion for the Case of
From the last section, we know that the root trajectory s(α) has better behavior at α = 1 in the case of b 2 − 4c < 0 than other cases. We confine to the case b 2 − 4c < 0 and further clarify the variation of roots s(α) according to the change rates of the argument θ(α), real part λ(α) and imaginary part ω(α) of roots s(α) at α = 1. The discussion is confined to the upper-half complex plane.
Substituting Equations (28), (39) and (40) into Equation (41), we acquire the derivative of the real part at α = 1 as a function of b and c, Further, the partial derivative is calculated from (43) as Thus, the equation p(b, c) = 0 uniquely determines a continuous implicit function Calculating the partial derivative with respect to b, we have It follows from the equation p(b, c) = 0 that Replacing ln c in Equation (44) by using Equation (45), we obtain Hence, the derivative of the implicit function c = ξ(b) is Therefore, the implicit function c = ξ(b) is monotonically increasing. Finally, by letting b → 0 + in Equation (45), we obtain ξ(0 + ) = 1. Thus, the curve L 1 of the implicit function c = ξ(b) divides the domain b 2 − 4c < 0 in the (b, c) plane into two regions. The signs of λ (1) are deduced from the equation λ (1) = p(b, c). The proof is completed.
Proof. In the domain b 2 − 4c < 0, if c ≥ 1, then ω (1) < 0 from Equation (46). Now, we look into the function q(b, c) in the region 0 < b < 2, b 2 /4 < c < 1. For any fixed b (0 < b < 2), the limits hold Calculating the partial derivatives, we have where the equation q(b, c) = 0 is used. Now, we estimate the second partial derivative as follows, From the above analyses, the curves L 1 , c = 1 and L 2 subdivide the domain b 2 − 4c < 0 in the (b, c) plane into four disjoint regions, as shown in Figure 10. The region IV is narrow and is magnified in Figure 11. The characteristics of the four regions in the domain b 2 − 4c < 0 are listed in Table 1.  Figure 11. Local zoom of the rectangular area 0 ≤ b ≤ 2, 0 ≤ c ≤ 1 in Figure 10. Table 1. Characteristics of four regions in the domain b 2 − 4c < 0.

Regions Constraint Conditions
Characteristics of Roots at α = 1 Samples (b, c) and Root Trajectories     In region I, we take a sample b = 3, c = 15 and depict the trajectory of roots s(α) in Figure 12. In Figure 8, b = 4, c = 4.4, and in Figure 9, b = 4, c = 8, which are all in region II. In region III, we take a sample b = 0.9, c = 0.6, and the trajectory of roots is shown in Figure 13. In region IV, we take b = 0.6, c = 0.1 and display the trajectory of roots in Figure 14.
Note that region IV is in the immediate vicinity of the critical case b 2 − 4c = 0. This means that only when s(1) is close enough to the real axis, its imaginary part may satisfy ω (1) > 0.

Conclusions
We consider the roots of the characteristic Equation (4) of the fractional vibration Equation (1) on the principal Riemann surface for the three cases of the discriminant b 2 − 4c: > 0, = 0 and < 0, and with the range of α covering the interval 0 ≤ α ≤ 2. Particular attention is paid to the varying of the roots s(α) as a function of the order α near α = 1 on the upper-half complex plane. We find that the root trajectories of the characteristic equation with α varying have different behaviors in the three cases of the discriminant. The residue contribution and Hankel integral contribution serve as two solution components of the impulsive response of the fractional vibration equation. It is found that the changing pattern of the solution components as α → 1 depends on the sign of the discriminant b 2 − 4c.
In Section 2, we prove that for any b > 0, c > 0 and α ∈ (0, 1) ∪ (1, 2), the characteristic Equation (4) has a pair of conjugated simple complex roots with a negative real part. For the three cases of the discriminant b 2 − 4c, the variations of the argument and modulus of the roots according to α are clarified. For the integer-order cases, α = 0, 1 and 2, the inverse Laplace transform gives the known analytic results. In Section 3, the trajectories of the roots s(α) on the upper-half complex plane are analyzed and simulated for the three cases of the discriminant b 2 − 4c. For the case of b 2 − 4c > 0, the root s(α) as a function of α is discontinuous at α = 1, and the left and right limits equal the two different real roots of the case α = 1 as in Equations (37) and (38). For the case of b 2 − 4c = 0, the root s(α) is continuous at α = 1 but non-differentiable at α = 1. For the case of b 2 − 4c < 0, the root s(α) is continuous and also differentiable at α = 1. In Section 4, the particular analyses for the case of b 2 − 4c < 0 are conducted. The trajectories of the roots s(α) are further clarified in the domain b 2 − 4c < 0 on the (b, c) plane according to the change rates of the argument θ(α), real part λ(α) and imaginary part ω(α) of the root s(α) at α = 1. For this purpose, the domain b 2 − 4c < 0 is subdivided into four disjoint regions to clarify the trajectories of the roots.
In Section 5, the residue contribution and Hankel integral contribution to the impulsive response of the fractional vibration equation are considered. For the cases of b 2 − 4c ≥ 0, the left and right limits of the residue contribution z R (t; α) at α = 1 produce a jump and neither equals z(t; 1). The Hankel integral contribution z H (t; α) does not vanish as α → 1 but fills up the jump such that z(t; α) = z R (t; α) + z H (t; α) → z(t; 1) as α → 1. Nevertheless, for the case of b 2 − 4c < 0, the variations of z R (t; α) and z H (t; α) are different from the cases b 2 − 4c ≥ 0, and we derive that the residue contribution z R (t; α) approaches z(t; 1) as α → 1, and so the Hankel integral contribution z H (t; α) approaches 0 as α → 1.
For the two solution components in Equations (50) and (51), the residue contribution represents a damping oscillation, while the Hankel integral contribution provides a monotonic recovering. Moreover, asymptotic behaviors can be obtained conveniently from the solution components [23]. In the residue contribution (50), the real part λ(α) and imaginary part ω(α) of the root s(α) describe the decaying rate of the amplitude and the oscillation frequency, respectively. These results could be helpful for revealing the relationship between the model parameters b, c and α and the oscillation properties. Finally, we indicate that there are few reports on the root trajectories of an equation with a fractional power.