Different Stochastic Resonances Induced by Multiplicative Polynomial Trichotomous Noise in a Fractional Order Oscillator with Time Delay and Fractional Gaussian Noise

A general investigation on the mechanism of stochastic resonance is reported in a timedelay fractional Langevin system, which endues a nonlinear form multiplicative colored noise and fractional Gaussian noise. In terms of theoretical analysis, both the expressions of output steady-state amplitude and that of the first moment of system response are obtained by utilizing stochastic averaging method, fractional Shapiro and Laplace methods. Due to the presence of trichotomous colored noise, the excitation frequency can induce fresh multimodal Bona fide stochastic resonance, exhibiting much more novel dynamical behaviors than the non-disturbance case. It is verified that multimodal pattern only appears with small noise switching rate and memory damping order. The explicit expressions of critical noise intensity corresponding to the generalized stochastic resonance are given for the first time, by which it is determined that nonlinear form colored noise induces much more of a comprehensive resonant phenomena than the linear form. In the case of slow transfer rate noise, a newfangled phenomenon of double hypersensitive response induced by a variation in noise intensity is discovered and verified for the first time, with the necessary range of parameters for this phenomenon given. In terms of numerical scheme, an efficient and feasible algorithm for generating trichotomous noise is proposed, by which an algorithm based on the Caputo fractional derivative are applied. The numerical results match well with the analytical ones.


Introduction
Based on the Markov process without after-effect, classical stochastic resonance theory grounded on ordinary Langevin equation and FPE tends to be perfect gradually. However, more and more systems under complex environments exhibit memory effects, such as longrange correlation and historical dependence in time. At this time, the traditional stochastic resonance theory with ODE is obviously no longer applicable as it often needs to be modeled by non-traditional mathematical models such as fractional order differential equation. Early scholars tried to use the research method corresponding to a normal diffusion system to approximate this type of issue, e.g., fractional relaxation equation [1], fractional diffusion equation [2] and fractional Fokker-Planck equations [3]. Due to the characteristics of the fractional calculus equation itself, the systematic description of stochastic resonance in the fractional system cannot be achieved by using analytical algorithm in most cases at present. γ(t − t )v(t )dt [25] rather than a constant coefficient, the damping force at any time t is acquired by integrating the velocity at all times before t according to the weight function γ. From the above consideration, it can be further turned into the generalized Langevin equation where γ(t) is the memory damping kernel, describing delayed effect of the friction. In many physical and biological environments, a power-law form memory kernel is often used to describe the dependence between viscosity force and particles' historical velocity, as follows The shape of damping kernel function with the friction coefficient δ = 1 is shown in Figure 1, it finds that γ(t) has a power-law attenuation trend with the speed of q, and it would degenerate into Delta function when q approach 1. Inserting (6) into (5) and using the definition of Caputo derivative, Equation (5) can be further transformed to the following FDE The well-known Brownian motion B(t) is generally used in the traditional stochastic dynamics, to model the normal diffusion with mean-square displacement (MSD) x 2 (t) = t induced by the random factor accompanied with all sorts of systems, it has further been confirmed to be quite serviceable in many situations. However, it is clunky when dealing with some other ubiquitous situations, i.e., the Brownian motion theory does not apply to the anomalous diffusion process [26], such as the long-range correlation behavior in historydependent viscoelastic medium under the background of anomalous diffusion, with a magnitude of MSD the same as t α , the diffusion exponents α lies between zero and one for a subdiffusion situation and beyond one for a superdiffusion case, respectively [27]. By this reason, a noise with a long time memory effect is frequently adopted to describe the non-Markov dynamic behavior of anomalous diffused particles, and the correlation function of W H (t) tends to be in a power-law form, as in Equation (7) [26]. According to the fluctuation dissipation theorem [25], the friction produced by the microscopic forces of the system depends entirely on the correlation of the random forces: W H (t)W H (s) = k B Tγ(|t − s|), with k B the Boltzmann constant and T the absolute temperature. Substituting (6) into the fluctuation dissipation equality gives the power-law form autocorrelation function of random force H D denotes the noise Here we consider W H (t) to be the fractional Gaussian noise (FGN), and similar to the definition pattern of GWN, FGN is given by the generalized differential of the fractional Brownian motion (FBM): W H (t) = √ 2D H dB H (t)/dt. D H denotes the noise strength, H denotes the Hurst index with 0 < H < 1. B H (t) is reduced to a standard Brownian motion for H = 0.5, which has independent increments. For H > 0.5 the increments are positively correlated, while they are negatively correlated for the opposite situation, and more details of FBM can be found in Ref. [28]. The mean value and autocorrelation function of FBM read as: where the correlation function of FGN can be obtained by the covariance of FBM When t = s, 0 < H < 1/2 corresponds to the subdiffusion movement with a negative correlation function, 1/2 < H < 1 corresponds the opposite situation, and it turns to the impulse function W H (t)W H (s) = 2D H δ(t − s) for H = 1/2, say, the normal diffusion case. Utilizing the fluctuation dissipation theorem, again one gets The Hurst parameter should be constrained in (0.5, 1) to guarantee the mean-square value W H 2 (t) to be positive, thus the memory kernel is positive as well. Combining (10) and (8) one has q = 2 − 2H, that is, the attenuation parameters of memory kernel function correspond exactly to the order of the fractional derivative. An intensity setting of FGN as D H = k B Tγ 0 /Γ(2H + 1) can provide a result consistent with (6), which reveals the practical significance of considering FGN in an anomalous diffusion system.
In this paper, we consider a kind of anomalous diffusion system with small time delay x(t − τ), disturbed by nonlinear colored noise and internal FGN, which is given by the following integro-differential equation: x(t)dt + ω 0 2 g[x(t), x(t − τ)] = A cos Ωt + W H (t) (11) ..
x(t) and .
x(t) denote the acceleration and velocity, respectively. A cos Ωt is external harmonic excitation, ω 0 2 the inherent frequency. Ψ[ζ(t), N] (N ≥ 2) represents the nonlinear trichotomous noise disturbance on damping intensity, which manifests itself as the polynomial form Ψ[ζ(t), N] = N ∑ k=1 γ k ζ k (t). ζ(t) is symmetric trichotomous colored noise with three states, W H (t) is zero-mean FGN with autocorrelation function given by (9) and memory damping kernel given by (6), respectively. Considering the potential field force where ρ represents delay intensity. ρ = 0 and ρ = 1 correspond, respectively, to non-delay and global item delay systems [17], 0 < ρ < 1 corresponds to the system containing delay feedback, wherein ρ is also referred as the feedback gain [29]. For computational convenience sake assume m = 1, then based on the generalized Langevin equation given by Equation (11) and the definition of Caputo fractional derivative [30], the system can be rewritten as the following fractional Langevin equation ..

Theoretical Analysis of Response
With the constructed equation ζ(t)(ζ(t) − ∆)(ζ(t) + ∆) = 0, and the arbitrary power of ζ(t) can be expressed by m is a positive integer, with (13) the polynomial noise in Equation (12) can be simplified to (14) where " x " denotes just the value x rounded down to the nearest integer. For the case of N = 1, the effect of the linear Markovian symmetric trichotomous noise on stochastic resonance has been investigated [6]. From the expansion expression (13) it is clear that a N-order (N > 2) polynomial trichotomous noise can be simplified to an equivalent quadratic polynomial one, i.e., the first two orders play a post role in the nonlinear trichotomous noise with arbitrary (adequate) high-level power, thus for convenience, this paper adopts the quadratic form given in (14), and all the other cases with high-power exponents can be expressed by the first two orders, thus can be discussed similarly. Moreover, for notational convenience, hereinafter the subscript "N" in α N and β N will be dropped, and the simplified versions α and β will be used to describe the linear coefficient and nonlinear coefficient of the polynomial function Ψ(ζ(t), N) of the trichotomous colored noise, respectively. As for the function g(x(t), x(t − τ)) of the small time delay, when ρ = 0 no delay appeared, and this mainstream situation has attracted widespread attention. Meanwhile the global delay case happens for ρ = 1 and the stiffness term disappears [15,17]. The scenario of 0 < ρ < 1 refers to the time-delayed feedback in the system [29], the small time-delay pre-assumption enables us to expand the small delay approximation method (SDAM) [31] to the current fractional system case, say, g(x(t), x(t − τ)) and can be expressed as a Taylor expansion series with respect to the time delay around x(t), with the convergence precision the same order of magnitude as O(τ 2 ). In fact, such a Taylor expansion approximation of a function with respect to x(t) has been proved workable when dealing with a small delay term, no matter whether in an integer-order differential system [32] or a fractional-order system [33]. Making use of this approximation method, one gets: where ω 0 2 = ω 0 2 1 − τρω 0 2 . Combining Equation (15) and, as discussed below, Equation (14), we obtain a nondelayed stochastic differential equation, which is equivalent to Equation (12): As mentioned above in Introduction section, there are many alternative assessment indexes to quantify the SR phenomenon in stochastic systems. Among them, the amplitude gain of the output signal is often considered as the one that could be utilized the most, from an analytical point of view. When describing the generalized SR and bona fide one, naturally many researchers have adopted this index in dealing with SR problems for its simplicity and effectiveness in practice. In this paper an analytical expression of the amplitude gain of the output signal in regard to the system given by Equation (12) will be discussed in detail, however, before doing this, we first calculate the exact solution of the first-order moment of the fractional system. This operation is feasible since the displacement process x(t) could be always stationary [33]. In fact, the fractional particle in Equation (12) will always be pulled back to the origin sooner or later due to the existence of the harmonic potential V(x) = ( ω 0 2 x 2 )/2. For this purpose, the Shapiro-Loginov formula [8] should be taken account for the exponential correlative trichotomous process: and the fractional Shapiro-Loginov formula: Here, x(t) is considered as a function of ζ(t). Taking average of Equation (16), one gets The prerequisite condition η H (t) = 0 has been used in Equation (18), considering the properties of the trichotomous noise and using the transition probabilities and stationary probabilities, the mean-square value and the correlation function of ζ 2 (t) are given by ζ 2 (t) = 2∆ 2 p and It is apparent from the above autocorrelation function of the mean value ζ 2 (t) that ζ 2 (t) does not satisfy the conditions of the Shapiro-Loginov formula, however, the situation becomes gratifying for the adjustive variable ζ 2 (t) − 2∆ 2 p, that Thus, the random variable ζ 2 (t) − 2∆ 2 p is applicable for various kinds of Shapiro-Loginov formulas (of course also including the fractional version (3.5)) (19) leads to the new correlation factor ζ 2 (t)D q x(t) appear in (18): By substituting Equations (17) and (20) into (18), one gets the equation To solve the new factor ζ(t)x(t) , multiplying Equation (16) by ζ(t) and averaging over the gained equation, one arrives at in Equation (22) the uncorrelated postulate between ζ(t) and the FGN η H (t) has already been adopted, i.e., ζ k (t)η H (s) = 0, k = 1, 2 . . . Reusing the Shapiro-Loginov formula, we have Substituting (17), (20) and (23) into Equation (22) yields In order to solve the term ζ 2 x , multiplying Equation (16) by ζ 2 (t) and the averaging the gained results, one obtains: the first averaged term in Equation (25) can also be derived by (19) and (23) From (26) the averaging value ζ 2 d 2 dt 2 x can be obtained as the expression Inserting (17), (20) and (27) into Equation (25), meanwhile, replacing ζ 3 (t), ζ 4 (t) for ∆ 2 ζ(t) and ∆ 2 ζ 2 (t), respectively, then (27) can be transformed to the harmonic terms on the right side of Equation (28) would not further be eliminated unless recurring to Equation (21), in fact, multiplying Equation (21) by 2p∆ 2 and subtracting the gained equation from Equation (28), one can obtain: The collection of (21), (24) and (29) truly consists of a close linear system of secondorder differential equations for three variables m 1 (t) = x(t) , m 2 (t) = ζ(t)x(t) , m 3 (t) = ζ 2 (t)x(t) . By virtue of the Laplace transform technique, this equation set can be equivalently changed to a corresponding linear algebraic equation set, given as follows Displacement's first-order moment can be formally expressed formally as bellow: (32) Applying the inverse Laplace transformation technique on Equation (31) and based on the corresponding uniqueness condition [34], the averaging mean value (the first-order moment) of the displacement of Equation (12) can be represented in the following form: " * " the convolution operation, with the inverse Laplace operator the primitive function y(t) can be calculated by the integral of its corresponding Laplace transformation function Y(s): Based on the stability theory [35], the stability of the solution expressed as (33) can be guaranteed under the circumstances that all the possible solutions of the equation I(s) = 0 have roots with a negative or zero real part. When considering the asymptotic behavior of the functions on the right hand of (33), i.e., at a long-time limit t → ∞ , by making use of the Tauberian theorem [36], one gets the estimation: 1+q) ), thus the memorizing effect of the initial conditions can be disregarded when the system has a long enough evolute, i.e., the stationary solution of the first-order moment of the displacement in a system described by Equation (12) can be obtained by (33).
the complex susceptibility of the stochastic system given by Equation (12) can be given based on the form of (34) [37]: ℵ(Ω) = ∞ 0 H 10 (t)e iΩt dt = H(−iΩ) = ℵ 1 (Ω) + ℵ 2 (Ω) · i, where i denotes the imaginary unit with i 2 = −1, ℵ 1 and ℵ 2 representing the real and imaginary parts of the complex susceptibility, respectively. Moreover, with the help of the linear response theory [38], (34) can be rewritten in the following harmonic pattern as well x(t) as = F cos(Ωt + φ). With the output amplitude F = A/|ℵ| and the phase shift given as, φ = arctan(−ℵ 2 /ℵ 1 ), it is worth nothing that H 10 (iΩ) = H 10 (−iΩ) = ℵ = ℵ 1 − ℵ 2 i, thus F and φ can also be calculated by F = A H 10 (iΩ) , φ = arg( H 10 (iΩ)). The response amplitude gain, say, the ratio of the amplitude of the output signal to the amplitude of input signal, then can be expressed as Through the expressions of the transfer function given in (31) and the coefficients below (29), after a series of simplification calculations, one can obtain R = e 1 2 + e 2 2 e 3 2 + e 4 2 , φ = arctan( where e k , k = 1, . . . , 4 are listed as follows: potential well of the system. The third one is output amplitude gain (OAG), which was proposed by Gitterman [41]. Traditional stochastic resonance refers to the optimal response of the system at a critical noise intensity when the system is subjected to both external and random forces. In recent years the concept has been expanded: SR in broad sense [41] proposed by Gitterman refers to the non-monotonic dependence of some functions of the output signal (such as SNR, autocorrelation function, amplitude gain, etc.) on noise or system parameters, showing that traditional SR is actually included in the concept of generalized stochastic resonance (GSR). Another extended concept is bona fide SR, which indicates the peak phenomenon of output signal amplitude induced by the fluctuation of signal frequency [42]. In this section, based on the theoretical output amplitude gain results in the previous section, different types of stochastic resonance cases will be analyzed for the fractional perturbation system with small hysteresis.
Before discussing the possible stochastic resonance phenomenon in the system (12), we first examine the case of no polynomial trichotomous noise disturbance, i.e., α = β = 0 or ∆ = 0, then the corresponding system equation reads ..
using the small time delay approximation method (15), Equation (36) is further transformed into a delay-free equation, taking the average of both sides of the obtained equation, one has d 2 after the Laplace transformation of Equation (37) with respect to the variable x(t) , the mean value of steady-state for t → ∞ can be calculated as A * and φ * denote the corresponding amplitude and phase shift, respectively.
for any fractional order satisfying 0 < q < 1, the steady-state output amplitude of the system (37) can be obtained as In fact, according to the linear response theory, if one substitutes x(t) = F * cos(Ωt + φ * ) into (37) directly, a same result as (38) can also be obtained by utilizing the undetermined coefficient method. In addition, the same situation happens when inserting α = β = 0 into the expression (35), showing that system (36) is actually a special case of system (12). Figure 2 shows the bona fide SR phenomenon in system (36) that was induced by external excitation frequency Ω under different flatness parameter q, with system parameters A = 1, ρ = 0.5, γ 0 = 1, ω 0 = 1, τ = 0.2. It finds that the resonance of the system output signal is striking when q is small, while for the opposite situation F * changes slowly despite the occurrence of peak Figure 2b describes the combined effect of external excitation frequency and flatness parameters on OAG of the system. With the increase in q, the single-resonance point Ω SR gradually increases and F * tends to be flat. Indeed, the critical resonant frequency Ω SR must satisfy flat. Indeed, the critical resonant frequency On the one hand, expression (38) when the index is close to one. According to (39), the critical SR frequency of the system c On the one hand, expression (38) implies F * ≈ A/ ω 0 2 − Ω 2 + γ 0 when the Hurst index is close to one. According to (39), the critical SR frequency of the system can be obtained as Ω SR ≈ ω 0 2 + γ 0 , According to the form of ω 0 2 , single-peak SR will certainly emerge in the system for small delay, and F * max → ∞ when q approaching to 0. Therefore, in Figure 2b F * exhibits a saltation within a certain frequency range for q → 0 . On the other hand, con-  Figure 3 shows the bona fide SR in the disturbed system (12) with a small flatness parameter q = 0.05 under different noise transfer rates, three-peak resonance, double-peak resonance and single-peak resonance appear at v = 0.01, v = 0.2 and v = 3, respectively. It is called stochastic multiresonance (SMR) if SR manifests as more than one peak. Thus in Figure 3 one can say that the external excitation frequency induces bona fide SMR. The noise transfer rate is sometimes called the correlation rate, which reflects the memory effect of noise. A small correlation rate means a long correlation time, and the switching between different states of noise will be extremely slow, leading to excessive correlation time results in the resonance peak at three different frequencies. In fact, in this case a state transition of the colored noise will be accompanied by a long-time stagnation, during which the damping coefficient of the system is actually fixed. Therefore, the random damping disturbance coefficient of the original system (12) can be divided into three cases in a long time.
Replacing (40) for γ 0 in the results without damping disturbance (38), and using resonance generation condition (39), it can be estimated that the positions of the three peak points that may appear in the slow-switching-rate system for Hurst index approaching to one ( Figure 3a) are Replacing (40) for 0 γ in the results without damping disturbance (38), and resonance generation condition (39), it can be estimated that the positions of the thre points that may appear in the slow-switching-rate system for Hurst index approach one ( Figure 3a) are  2 When H is close to 1 2 , the external excitation frequency corresponding three resonance extreme points can also be calculated in the same way, given as As v increases, the noise switching between the three states becomes mo more frequent, resulting in the three peak points gradually converging and even merging into one overlapping peak point. The influence of the transfer rate on th fide SR peak is given in Figure 3d, it can be seen that when v keep growing, the u resonance peak of the system will also increase gradually. The system presents at single-peak stochastic resonance phenomenon within the range of the combined p ter pairs ( ) , v Ω . It should be pointed out that all the system parameters considered section meet the stability condition ( ) The output amplitude of the system can also be approximately estimated for th switching noise scenarios, as When H is close to 1/2, the external excitation frequency corresponding to the three resonance extreme points can also be calculated in the same way, given as As v increases, the noise switching between the three states becomes more and more frequent, resulting in the three peak points gradually converging and eventually merging into one overlapping peak point. The influence of the transfer rate on the bona fide SR peak is given in Figure 3d, it can be seen that when v keep growing, the unique resonance peak of the system will also increase gradually. The system presents at least a single-peak stochastic resonance phenomenon within the range of the combined parameter pairs (Ω, v). It should be pointed out that all the system parameters considered in this section meet the stability condition I(0) > 0.
The output amplitude of the system can also be approximately estimated for the fast-switching noise scenarios, as Through calculation, it is not difficult to find that (41) exactly corresponds to the amplitude of the steady-state solution of the noiseless system (12) with the equivalent damping coefficient by solving dF(Ω)/dΩ = 0 in Equation (41), one has Ω SR = ω 0 2 + γ 0 + 2p∆ 2 β, the single-peak-resonance always exists in the system, which is consistent with the single peak point position Ω c = 1.76 in Figure 3c. For q → 1 , only single-resonance may occur at In the case of the slow-switching case v = 0.01, and the bona fide SR phenomenon under different flatness parameters is shown in Figure 4. For q = 0.01, the bona fide SR occurs at three critical frequencies; 1.38, 1.87 and 2.08, as q increases, triple resonant peaks gradually converge to double and single peaks, which is similar to the situation in Figure 3. The difference is that F decreases with the increase in q. When q is large enough, the resonance phenomenon disappears finally. Consider the other case q → 1 (we might as well set q = 0.99), α = 2, v = 3, Figure 5 shows SR phenomena under different polynomial noise coefficients. Since the condition (43) is not met for β = 2, the variation in Ω cannot induce bona fide SR. In fact, due to the effect of polynomial nonlinear multiplicative trichotomous noise, the resonance conditions are determined by many parameters. When nonlinear noise coefficient and damping coefficient were taken as control parameters and other parameter values fixed in (43), the parameter distribution of the occurrence of bona fide SR under the influence of combined parameters was given in Figure 5b. The critical value of γ 0 corresponding to the resonant peak is almost monotonically dependent on the critical value of β. When γ 0 > 1.82, any β will no longer meet the condition (43), and the external excitation frequency of the system will no longer induce the bona fide SR phenomenon.
In the case of the slow-switching case 0.01 v = , and the bona fide SR phen under different flatness parameters is shown in Figure 4. For 0.01 q = , the bona occurs at three critical frequencies; 1.38 , 1.87 and 2.08 , as q increases, triple r peaks gradually converge to double and single peaks, which is similar to the situ Figure 3. The difference is that F decreases with the increase in q . When q enough, the resonance phenomenon disappears finally. Consider the other case (we might as well set Figure 5 shows SR phenomen different polynomial noise coefficients. Since the condition (43) is not met for  variation in  cannot induce bona fide SR. In fact, due to the effect of poly nonlinear multiplicative trichotomous noise, the resonance conditions are determ many parameters. When nonlinear noise coefficient and damping coefficient we as control parameters and other parameter values fixed in (43) In the case of the slow-switching case 0.01 v = , and the bona fide SR phenomenon under different flatness parameters is shown in Figure 4. For 0.01 q = , the bona fide SR occurs at three critical frequencies; 1.38 , 1.87 and 2.08 , as q increases, triple resonant peaks gradually converge to double and single peaks, which is similar to the situation in Figure 3. The difference is that F decreases with the increase in q . When q is large enough, the resonance phenomenon disappears finally. Consider the other case 1 q → (we might as well set Figure 5 shows SR phenomena under different polynomial noise coefficients. Since the condition (43) is not met for 2  = , the variation in  cannot induce bona fide SR. In fact, due to the effect of polynomial nonlinear multiplicative trichotomous noise, the resonance conditions are determined by many parameters. When nonlinear noise coefficient and damping coefficient were taken as control parameters and other parameter values fixed in (43), the parameter distribution of the occurrence of bona fide SR under the influence of combined parameters was given in Figure 5b. The critical value of 0  corresponding to the resonant peak is almost monotonically dependent on the critical value of  . When 0 1.82

 
, any  will no longer meet the condition (43), and the external excitation frequency of the system will no longer induce the bona fide SR phenomenon.    Figure 6 shows the influence of noise intensity on bona fide SR of the system with respect to the parametric region q − v, under the polynomial trichotomous noise coefficient setting α = 0.4, β = 2. Considering the system parameters ρ = 0.5, γ 0 = 1, ω 0 = 1, p = 0.3, τ = 0.2, it can be seen from Equation (38) that the system always exhibits monostable bona fide SR phenomenon when the trichotomous noise is absent. By comparing Figure 6a-c, it is apparent that appropriate noise intensity is conducive to the occurrence of SMR. When the noise amplitude is small (Figure 6b), only single-resonance appears in the highest range of parameter values (light blue area in the figure), with non-resonance appearing around 0.88. When (q, v) are close to zero, a small range of two-peak resonance (cyan) and three-peak resonance (green) regions appear, indicating that thedditionn of trichotomous colored noise has enriched the stochastic resonance phenomenon of the system. When ∆ increases to one, the proportion of the double-resonance region and triple-resonance region increase significantly, which reflects the fact that the SMR phenomenon of the system can be enhanced by the appropriate noise amplitude. However, it is not the case that the stronger the noise is, the stronger the resonance of the system will be: excessive noise intensity (Figure 6c) results in the disappearance of SMR in the system. Moreover, excessive noise intensity will even lead to the non-stochastic resonance (white) with a wider range of parameters than the noise-free case. Figure 6 shows the influence of noise intensity on bona fide SR of the system with respect to the parametric region qv − , under the polynomial trichotomous noise coefficient setting 0.4  = , 2  = . Considering the system parameters , it can be seen from Equation (38) that the system always exhibits monostable bona fide SR phenomenon when the trichotomous noise is absent. By comparing Figure 6a-c, it is apparent that appropriate noise intensity is conducive to the occurrence of SMR. When the noise amplitude is small (Figure 6b), only single-resonance appears in the highest range of parameter values (light blue area in the figure), with nonresonance appearing around 0.88 . When

( )
, qv are close to zero, a small range of twopeak resonance (cyan) and three-peak resonance (green) regions appear, indicating that the addition of trichotomous colored noise has enriched the stochastic resonance phenomenon of the system. When  increases to one, the proportion of the doubleresonance region and triple-resonance region increase significantly, which reflects the fact that the SMR phenomenon of the system can be enhanced by the appropriate noise amplitude. However, it is not the case that the stronger the noise is, the stronger the resonance of the system will be: excessive noise intensity (Figure 6c) results in the disappearance of SMR in the system. Moreover, excessive noise intensity will even lead to the non-stochastic resonance (white) with a wider range of parameters than the noisefree case.

Parameter Induced GSR
Firstly, the influence of time delay  is investigated. In the case of a system with deterministic fractional damping (42), the peak point SR  of the system output amplitude corresponds to the root of the equation solving (44) and giving rise to the critical delay value ( ) Inserting (45) into Equation (38), the optimal output amplitude at the critical delay SR  is given by , which can also be referred to as the time-delay induced intrinsic resonance peak occasionally. The dependence of the output amplitude of the system without trichotomous noise disturbance on the time delay is depicted in Figure 7. The single-peak of GSR appears at 0.64  = with the maximum output amplitude max 12.8 , which is consistent with the critical delay analysis result in (45).
Examining the disturbed situation to be a comparison, based on the prediction result (35) of system (12), Figure 8 shows the relationship between system output amplitude and delay, and other system parameters

Parameter Induced GSR
Firstly, the influence of time delay τ is investigated. In the case of a system with deterministic fractional damping (42), the peak point τ SR of the system output amplitude corresponds to the root of the equation solving (44) and giving rise to the critical delay value Inserting (45) into Equation (38), the optimal output amplitude at the critical delay τ SR is given by F * = 1/γ 0 Ω q sin(qπ/2), which can also be referred to as the time-delay induced intrinsic resonance peak occasionally. The dependence of the output amplitude of the system without trichotomous noise disturbance on the time delay is depicted in Figure 7. The single-peak of GSR appears at τ = 0.64 with the maximum output amplitude F * max = 12.8, which is consistent with the critical delay analysis result in (45). Examining the disturbed situation α = 0.5, β = 1.2, p = 0.35, ∆ 2 = 0.6 to be a comparison, based on the prediction result (35) of system (12), Figure 8 shows the relationship between system output amplitude and delay, and other system parameters read A = 1, ρ = 0.8, γ 0 = 0.5, ω 0 = 1, q = 0.1, Ω = 0.99. For a different noise correlation rate and noise intensity, variation in delay always causes GSR in the system, and can also induce the double-and triple-SMR phenomenon. When v = 0.01, the damping coefficient of the system can be regarded as undisturbed over long periods, as the transfer rate of the trichotomous noise is very slow. By substituting the three cases of (40) into (45), one can obtain the possible critical time-delay of the resonance peak corresponding to ζ(t) = ∆, ζ(t) = 0 and ζ(t) = −∆.
with γ 1 = γ 0 − α∆ + β∆ 2 , γ 2 = γ 0 + α∆ + β∆ 2 . Substituting the parameters into (46), the result is consistent with the data annotation value in Figure 8. Triple-peak GSR is converted to a double-peak GSR with the increase in v. When v is large enough, the critical delay of the single-peak GSR under the fast-switching noise case can be calculated as follows: rate and noise intensity, variation in delay always causes GSR in the system, and can also induce the double-and triple-SMR phenomenon. When 0.01 v = , the damping coefficient of the system can be regarded as undisturbed over long periods, as the transfer rate of the trichotomous noise is very slow. By substituting the three cases of (40) into (45), one can obtain the possible critical time-delay of the resonance peak corresponding to  The fact that (47) is consistent with the peak position in Figure 8a. Figure 8b shows that despite the long correlation time of noise, the three equivalent damping coefficients corresponding to (45) will be very close due to the low noise intensity ( 2 0.01 = ) in the system, so that the inherent resonance peak could only exist at the critical delay given by  , the damping coefficient of the system can be regarded as undisturbed over long periods, as the transfer rate of the trichotomous noise is very slow. By substituting the three cases of (40) into (45), one can obtain the possible critical time-delay of the resonance peak corresponding to ( )   The fact that (47) is consistent with the peak position in Figure 8a. Figure 8b shows that despite the long correlation time of noise, the three equivalent damping coefficients corresponding to (45) will be very close due to the low noise intensity ( 2 0.01 = ) in the system, so that the inherent resonance peak could only exist at the critical delay given by  The fact that (47) is consistent with the peak position in Figure 8a. Figure 8b shows that despite the long correlation time of noise, the three equivalent damping coefficients corresponding to (45) will be very close due to the low noise intensity (∆ 2 = 0.01) in the system, so that the inherent resonance peak could only exist at the critical delay given by (45). When ∆ 2 = 0.3, γ 0 is very close to γ 1 ≈ 0.58, while there is a large gap between γ 0 and γ 2 ≈ 1.14. Therefore, the peak position of GSR should be τ SR1 and τ SR3 in (46), and it is not difficult to get that to τ SR1 = 0.642 and τ SR3 = 1.453, which is consistent with the resonance points of numerical results labeled in the figure. When ∆ 2 = 0.6, the intensity of the slow switching noise is large enough, and the three equivalent damping coefficients γ 0 , γ 1 and γ 2 , differ greatly from each other. At this time, three optimal peaks appear, respectively, at the three positions in (46), thus the system delay induces the three-peak SMR of the output amplitude. Based on the above analysis, it can be determined that when the system time delay is the only control parameter, a small noise transfer rate and an appropriate noise intensity would make the occurrence of multi-peak GSR phenomena more possible. In order to analyze the influence of parameters v and ∆ more carefully, different fractional orders are selected to give the joint dependence between GSR induced by delay and parameter v, ∆ 2 . It is found that SMR tends to appear for small q (cyan and green areas in Figure 9a), the single-peak GSR region becomes larger and larger as q increases, while the phenomenon of multi-peak GSR gradually disappears.
SMR of the output amplitude. Based on the above analysis, it can be determined tha the system time delay is the only control parameter, a small noise transfer rate appropriate noise intensity would make the occurrence of multi-peak GSR phen more possible. In order to analyze the influence of parameters v and  more ca different fractional orders are selected to give the joint dependence between GSR in by delay and parameter ( ) 2 , v  . It is found that SMR tends to appear for small q and green areas in Figure 9a), the single-peak GSR region becomes larger and large increases, while the phenomenon of multi-peak GSR gradually disappears. The aforementioned results confirmed that the noise intensity has a sign influence on both the bona fide SR induced by external excitation frequency (Figure GSR induced by time delay (Figure 8b). As such, we examined the change in the output amplitude when 2  varied with the other system parameters that were f As in Figure 10, variation in noise intensity leads to peak phenomenon of system amplitude under the slow noise switching rate condition of 0.01 v = , analogous traditional form of SR that occurs in the system under the measurement index SNR. 10a shows that the noise-induced stochastic resonance phenomenon presents d forms under the different nonlinear trichotomous noise factor  : One uniqu The aforementioned results confirmed that the noise intensity has a significant influence on both the bona fide SR induced by external excitation frequency ( Figure 6) and GSR induced by time delay (Figure 8b). As such, we examined the change in the system output amplitude when ∆ 2 varied with the other system parameters that were fixed as A = 1, γ 0 = 0.5, α = 1, ω 0 = 1, p = 0.35, q = 0.5, v = 0.01, Ω = 0.99, τ = 0.2. As shown in Figure 10, variation in noise intensity leads to peak phenomenon of system output amplitude under the slow noise switching rate condition of v = 0.01, analogous to the traditional form of SR that occurs in the system under the measurement index SNR. Figure 10a shows that the noise-induced stochastic resonance phenomenon presents different forms under the different nonlinear trichotomous noise factor β: One unique peak F max = 3.938 exists if the nonlinear part of colored noise is not considered (β = 0), and the single-peak SR only occurs when the trichotomous noise switches to ζ(t) = ∆. Indeed, by inserting β = 0 and (40) into Equation (38), one can obtain a real solution only if the equivalent damping coefficient is equal to γ 2 Two amplitude spikes appear in the system when β increases to 0.5. Similarly, by virtue of (38) and (40), it can be calculated that the system output amplitude may have a peak only when the trichotomous noise switches to −∆, and all possible extremum points are With κ = α 2 − 4βγ 0 Ω 2q + 4β Ω q Ω 2 − ω 0 2 cos(qπ/2), two resonance points corresponding to β = 0.5 were calculated as 0.279 and 2.136, respectively, by using Equation (49). In addition, there was a minimum point at ∆ = α/(2β), which was consistent with the numerical results of (12) (data marked in Figure 10a). This shows that, compared to the system with linear form trichotomous noise disturbed damping, appropriate nonlinear noise will induce new resonance phenomena. The resonant peak on the right side disappears when β = 1 as α Ω q < √ κ, while the left peak also disappears for β = 2, so no SR analogous to the traditional form appears in the system, which is due the fact that κ < 0 Figure 10b describes the resonance under different delay intensity: the only peak of output amplitude appears for the none-delay case (ρ = 0), which is due to κ < 0; when time delay exists, the phenomenon of double-peak stochastic resonance happens, with the positions of the two optimal points moving apart and decreasing gradually with the increase in delay intensity. Similar to the double-resonance in Figure 10a, the optimal points of SR are, respectively, given by ∆ SR1 and ∆ SR2 given in (49), and the minimum points of output amplitude all appear at ∆ = α/(2β) = 1. It means that on one hand the existence of time delay enriches the stochastic resonance induced by noise intensity; on the other hand, it also proved that compared with the full-delay case (ρ = 1), delay feedback system with both x(t) and x(t − τ) would exhibit more evident double-resonance. , which was consistent with the numerical results of (12) (data marked in Figure 10a). This shows that, compared to the system with linear form trichotomous noise disturbed damping, appropriate nonlinear noise will induce new resonance phenomena. The resonant peak on the right side disappears when 1  = as q   , while the left peak also disappears for 2  = , so no SR analogous to the traditional form appears in the system, which is due the fact that 0   Figure 10b describes the resonance under different delay intensity: the only peak of output amplitude appears for the none-delay case ( 0  = ), which is due to 0   ; when time delay exists, the phenomenon of double-peak stochastic resonance happens, with the positions of the two optimal points moving apart and decreasing gradually with the increase in delay intensity. Similar to the double-resonance in Figure 10a, the optimal points of SR are, respectively, given by  Under low noise switching rate condition, examining the effect of damping order q and steady-state probability p of noise on the relationship curve of F − ∆ 2 , we find the hypersensitive response induced by noise intensity, i.e., the system output amplitude shows a sensitive dependence on the noise intensity. Analytical results are shown in Figure 11a, with system parameters set at A = 1, ρ = 0.7, γ 0 = 0.5, α = 1, β = 0.5, ω 0 = 1, v = 0.01, Ω = 1, τ = 0.2. As the noise intensity increases from zero, the output amplitude reaches the resonant peak at ∆ 2 = 0.1923, followed by a sharp decrease in F with a slight increase in noise intensity, the first hypersensitive response occurs. The second hypersensitive response appears around ∆ 2 = 2.174 after reaching a minimum at ∆ 2 = 0.2765 and then slowly passes another peak point, as the output amplitude increases rapidly to the last peak point. The results obtained by substituting the parameters into (49) are consistent with the positions of the three resonance peak points in the figure, and the twice hypersensitive response phenomena appear after ∆ SR1 and before ∆ SR2 , respectively. Further comparison of Figure 11b,c shows that a two-time hypersensitive response only appears when the damping order and steady-state probability of noise are both small. When one of the two, or both of them, is larger, the amplitude curves before and after resonance positions ∆ SR1 and ∆ SR2 exhibit a more gentle property than those in Figure 11a, and the peak point ∆ SR3 = α/(2β) gradually disappears and eventually becomes the minimum point. At this juncture, noise-induced triple-SR is transformed to a double-SR, thus hypersensitive response disappears. In order to show the joint influence of parameter q and p concretely, different resonance regions that may appear in F − ∆ 2 relation curve under all parametric conditions of 0 < q < 1 and 0 < p < 1/2 are given in Figure 11d. The results suggest that two-time hypersensitivity response may occur only when parameters satisfy p < 0.219 and q < 0.276 simultaneously.
gentle property than those in Figure 11a, and the peak point ( ) 3 2 SR  = disappears and eventually becomes the minimum point. At this juncture, noi triple-SR is transformed to a double-SR, thus hypersensitive response disappea to show the joint influence of parameter q and p concretely, different regions that may appear in 2 F − relation curve under all parametric con 01 q  and 0 1 2 p  are given in Figure 11d. The results suggest tha hypersensitivity response may occur only when parameters satisfy 0 p  0.276 q  simultaneously. Considering the influence of fractional order q on the system output amplitude, it is difficult to give the explicit expression of the critical value q SR corresponding to the zero of dF(q)/dq since the output amplitude F is a complex function of q. Even so, dependence of the amplitude of the system output signal on the flatness parameter can be obtained by calculating the expression numerically, given in Figure 12 with parameters A = 1, ρ = 0.8, α = 1, ω 0 = 1, p = 0.35, v = 0.01, Ω = 0.99, τ = 0.2. The system response manifests as a monotone increasing function of q under the linear trichotomous noise disturbance. When the colored noise is in nonlinear form (β = 0.5), the output amplitude of the system peaks at q SR = 0.42, that is, the fractional damping order does make the system produce the phenomenon of GSR. When different damping intensity values are considered (Figure 12b), we determine that γ 0 has a great influence on the GSR phenomenon induced by damping order. Moreover, the output amplitude decreases with the increase in damping order for γ 0 = 0 and γ 0 = 0.75, while stochastic resonance occurs at q = 0.445 for γ 0 = 0.5, and a minimal resonance point F min = 1.844 emerges when γ 0 = 0.25. considered (Figure 12b), we determine that 0  has a great influence on the GSR phenomenon induced by damping order. Moreover, the output amplitude decreases with the increase in damping order for

Numerical Verification Scheme
With these aforesaid preliminaries mentioned in Section 2, we are now in a position to generate realizations of a stochastic process for the trichotomous noise, which is done in the following way: 1. Without loss of generality, assume the noise is located initially at

Numerical Verification Scheme
With these aforesaid preliminaries mentioned in Section 2, we are now in a position to generate realizations of a stochastic process for the trichotomous noise, which is done in the following way: 1. Without loss of generality, assume the noise is located initially at ∆ i at time t n , i.e., ζ n = ∆ i , i = 1, 2, 3. At (n + 1)th time point, to determine whether the process moves to the other two sites or merely remains at the same site, we consider the conditional probability matrix as given by (3). Then, it generates a uniformly distributed random number r n that is falling in [0, 1], which is then compared against the three values in the ith line of the conditional probability matrix. If r n < P i1 (h), then we accept the value of the noise ∆ 1 , i.e., ζ n+1 = ∆, and if r n > 1 − P i3 (h) we accept the value ∆ 3 or ζ n+1 = −∆, otherwise, it will be located at ∆ 2 , i.e., ζ n+1 = 0; 2. Taking the site at time t n+1 as the starting location again, marked as ζ n+1 = ∆ k , another new uniformly distributed random number r n+1 between [0, 1] is generated, and we conduct a comparison operation on r n+1 with the elements in the kth line of (3). If r n+1 < P k1 (h) or r n+1 > 1 − P k3 (h) then the value of noise at time t n+2 is confirmed as ζ n+2 = −∆ or ζ n+2 = ∆, otherwise we accept ζ n+2 = 0. By repeating this two-step procedure one can generate a sequence of random numbers ζ(t) switching between three values ∆, 0 and −∆.
It should be noted that the time interval between the two steps (h) is always fixed and be equal to ∆t, which should be much smaller than the correlation time of the noise (h << τ cor ). Let h = 10 −4 , ∆ = 5, and the total time T = 10, Figure 13 shows four realizations of the trichotomous noise under different steady-state probability parameters and different correlation times. Closer observation would reveal that a larger τ cor is often accompanied by a longer residence time of a particular state on average, with a larger p value that would lead to a shorter residence time at the origin during the noise switching between ∆ and −∆. Comparing the left and right figures, one finds that smaller relaxation time corresponds to a faster average switching rate between the three states. In order to examine the accuracy of the above numerical algorithm, we first calculate the arithmetic average value of 2000 sample results; the first moment of generating time series is determined to satisfy the zero-mean characteristic of trichotomous colored noise ( ζ(t) ∼ 10 −3 ). In addition, to check whether the results satisfy the theoretical autocorrelation function, we give the theoretical and numerical results of the normalized auto-correlation function ζ(t)ζ(t + s) / ζ 2 (t) under different characteristic relaxation time, as shown in Figure 14, wherein red markers represent the simulation results for 2000 samples. τ cor ( f it) is the exponential order of the fitting curves of discrete results. It can be seen that the relaxation time is very consistent with the theoretical setting, which in turn proves the feasibility (of this algorithm). ( ) cor fit  is the exponential order of the fitting curv discrete results. It can be seen that the relaxation time is very consistent with theoretical setting, which in turn proves the feasibility (of this algorithm).  Introducing the velocity variable ( ) yt, the original second-order FDE (12) can be reduced to first-order equation set:

auto-correlation function
According to the Caputo fractional derivative approximation method used in the previous work [43], (50) can be discretized as follows: is the trichotomous noise state at n t .

( )
is the state value of FBM at n t , which is generated by the Wood-Chan algorithm [44]. Based on the discrete Fourier method, this method provides a high computational speed, even in the Introducing the velocity variable y(t), the original second-order FDE (12) can be reduced to first-order equation set: .
According to the Caputo fractional derivative approximation method used in the previous work [43], (50) can be discretized as follows: n = 0, 1, 2 . . . N, N = T/h, h denotes the time step, T the total time, M = [τ/h], t n = nh, x n = x(t n ), y n = y(t n ), the initial condition in [0, τ] is set as x(t) = y(t) = 0. A n = A cos Ωt n , ζ n = ζ(t n ) is the trichotomous noise state at t n . B Hn = B H (t n ) is the state value of FBM at t n , which is generated by the Wood-Chan algorithm [44]. Based on the discrete Fourier method, this method provides a high computational speed, even in the case of large time step. Figure 15 shows the simulated time series of the process of FBM and their corresponding FGN under different Hurst indices. It can be seen that compared with the classical White Gaussian noise (H = 0.5), the larger H is, the smoother the sample path is. The situation of H < 0.5 is opposite to that of H > 0.5, moreover, the range of FGN decreases with the increase in Hurst index. Select h = 10 −3 and T = 1000, the simulation of discrete algorithm (51) for the output signal of the system under different intensity of nonlinear trichotomous noise are shown in Figure 16. Time domain and frequency domain simulation results is displayed for the first 1000 s, with parameters A = 1, ρ = 0.8, γ 0 = 0.5, α = 1, β = 1, ω 0 = 1, p = 0.35, v = 0.1, Ω = 0.99, τ = 0.2. By comparing (a) and (c), it can be seen that larger noise intensity does not drive the system to output signals with larger amplitude, and the results of ∆ 2 = 0.1 show stronger periodic characteristics than those of ∆ 2 = 1. (b) and (d) show the amplitude of the system output signal at each frequency, and it is found that it peaks at the external excitation frequency Ω = 0.99 in both cases. In the case of ∆ 2 = 1 and ∆ 2 = 0.1, the numerical results of the amplitude at the external excitation frequency are F num = 1.587 and F num = 2.545, respectively, and the theoretical results obtained by the analytical expression (35) are F = 1.592 and F num = 2.545, respectively, fully showing that the numerical results of a single sample are quite close to the theoretical prediction results.
In order to verify the feasibility of the numerical simulation algorithm (51), independent Monte Carlo repeated experiments should be conducted to obtain I samples, x i (t n ), i = 1, 2, . . . I, n = 1, 2, . . . N, N represents the total steps. The first-order moment of the system response at time t n can be given by the following formula: The Fourier coefficient of the response is taken as the measurement of the output amplitude, and the time series (52) is used to give the simulation results of the output amplitude as follows x i (t n ) sin Ωt n /N I where FFT denotes the fast Fourier transform.
noise intensity does not drive the system to output signals with larger amplitude results of 2 0.1 = show stronger periodic characteristics than those of 2 1 = (d) show the amplitude of the system output signal at each frequency, and it is fo it peaks at the external excitation frequency   Figure 17a. The time consuming T cost increases exponentially with the decrease in the time step h, and the approximate fitting relationship reads T cost ≈ 164e −1253h , the calculation time of a single sample increases greatly especially for h ≤ 10 −2 , which is bound to lead to a sharp increase in the experimental time of a large number of independent samples. In order to verify the convergence of numerical simulation algorithms (51) and (53), taking the output amplitude as a function of h, the simulation results based on (53) and theoretical results based on (35) are compared in Figure 17b. When h decreases, the numerical simulation amplitude gradually approaches the theoretical result F = 1.592, and tends to 1.603 without significant fluctuation when h ≤ 10 −2 , the inherent error is caused by the numerical discretization of the fractional system (12). Numerical results are no longer significantly affected by step size when h ≤ 10 −2 , so time step h = 0.01 is selected in the following part.
Under the same parameter setting as Figure 16, I = 1000, h = 10 −3 , the relationship between F num and the total simulation time T for ∆ 2 = 0.1 and ∆ 2 = 0.1 was examined by (53), depicted by Figure 18a,b, respectively. When T < 100, the simulation results differ greatly from the theoretical results due to the transient influence of the initial condition. When T is large enough, the accuracy of the numerical results is no longer significantly affected by the initial condition, yet fluctuates in a small range near the theoretical ones rather than gradually approaching them. In fact, due to the existence of nonlinear trichotomous noise, the output amplitude is actually a random variable. The statistical result (53) with a fixed sample number is bound to bring certain error. To measure this error, we still 0.1 = was examined by (53), depicted by Figure 18a,b, respectively. When   condition. When T is large enough, the accuracy of the numerical results is no significantly affected by the initial condition, yet fluctuates in a small range n theoretical ones rather than gradually approaching them. In fact, due to the exis nonlinear trichotomous noise, the output amplitude is actually a random varia statistical result (53) with a fixed sample number is bound to bring certain e measure this error, we still use the Monte Carlo method to calculate the mean error (MSE) between num F and the theoretical results given by (35) ( ) In virtue of (54), the mean square error of simulation results corresponding to ∆ 2 = 0.1 and ∆ 2 = 1 are given in Figure 18b,d, respectively. It is found that the mean square error tends to be stable with the increase in simulation time, σ → 7.3 × 10 −2 for ∆ 2 = 0.1 and σ → 4.2 × 10 −2 for ∆ 2 = 1, which reflects the reliability of the algorithm (51). Under the same parameter conditions, the influence of the number I of the Monte Carlo experiment on the steady-state output amplitude F and mean square error σ was investigated, and the simulation time was fixed T = 500. Comparison between ∆ 2 = 0.1 (Figure 19a,b) and ∆ 2 = 1 (Figure 19c,d) shows that, when the number of experimental samples is large enough, the simulated output amplitude oscillates within a small range of the theoretical one, and the mean square error is around 6 × 10 −2 , which indicates that the algorithm in this paper (51) is basically consistent with the theoretical analysis result given by (35).
Fractal Fract. 2022, 6, x FOR PEER REVIEW was investigated, and the simulation time was fixed 500 T = . Comparison 2 0.1 = (Figure 19a,b) and 2 1 = (Figure 19c,d) shows that, when the nu experimental samples is large enough, the simulated output amplitude oscillates small range of the theoretical one, and the mean square error is around 6 10 −  indicates that the algorithm in this paper (51) is basically consistent with the th analysis result given by (35).

Conclusions
In this paper, a general investigation into stochastic resonance in the f Langevin system with internal FGN is reported theoretically and numerically, the fractional Brown particle is driven by an external periodic excitat multiplicative polynomial colored trichotomous noise. The theoretical analys steady-state amplitude of the system output signal and the first moment of th response, relying on a stochastic vibrational mechanism, are obtained by exten stochastic average method and utilizing a fractional Shapiro-Loginov form fractional Laplace transformation law. Bona fide stochastic resonance induced frequency of periodic signal and generalized stochastic resonance induced by d noise intensity 2  and damping order q are analyzed in detail. An effic feasible algorithm for generating trichotomous noise is proposed, by which a n discretization algorithm for the fractional Langevin system is presented, numerical results agree well with the theoretical prediction.
The theoretical treatment provides a convenient and shortcut approach to r mechanism of two kinds of stochastic resonance. For bona fide SR, new SMR phe are found in the system with trichotomous noise (especially, the one with a slow s rate), which does not appear in the system without colored noise. Under appropri intensity, bimodal and trimodal bona fide SR are more likely to occur with switching colored noise and small fractional damping order ( Figure 6). In oth

Conclusions
In this paper, a general investigation into stochastic resonance in the fractional Langevin system with internal FGN is reported theoretically and numerically, wherein the fractional Brown particle is driven by an external periodic excitation and multiplicative polynomial colored trichotomous noise. The theoretical analysis of the steady-state amplitude of the system output signal and the first moment of the system response, relying on a stochastic vibrational mechanism, are obtained by extending the stochastic average method and utilizing a fractional Shapiro-Loginov formula and fractional Laplace transformation law. Bona fide stochastic resonance induced by the frequency of periodic signal and generalized stochastic resonance induced by delay τ, noise intensity ∆ 2 and damping order q are analyzed in detail. An efficient and feasible algorithm for generating trichotomous noise is proposed, by which a numerical discretization algorithm for the fractional Langevin system is presented, and the numerical results agree well with the theoretical prediction.
The theoretical treatment provides a convenient and shortcut approach to reveal the mechanism of two kinds of stochastic resonance. For bona fide SR, new SMR phenomena are found in the system with trichotomous noise (especially, the one with a slow switching rate), which does not appear in the system without colored noise. Under appropriate noise intensity, bimodal and trimodal bona fide SR are more likely to occur with a slow switching colored noise and small fractional damping order ( Figure 6). In other ways, dependence of GSR on joint parameters v, ∆ 2 ( Figure 9) indicates that multimodal GSR induced