Time-Delay Effects on the Collective Resonant Behavior in Two Coupled Fractional Oscillators with Frequency Fluctuations

Abstract

In this study, we propose coupled time-delayed fractional oscillators with dichotomous fluctuating frequencies and investigate the collective resonant behavior. Firstly, we obtain the condition of complete synchronization between the average behavior of the two oscillators. Subsequently, we derive the precise analytical expression of the output amplitude gain. Based on the analytical results, we observe the collective resonant behavior of the coupled time-delayed system and further study its dependence on various system parameters. The observed results underscore that the coupling strength, fractional order, and time delay play significant roles in controlling the collective resonant behavior by facilitating the occurrence and optimizing the intensity. Finally, numerical simulations are also conducted and verify the accuracy of the analytical results.

1. Introduction

Ever since Benzi et al. [1] first introduced the concept of “stochastic resonance” (SR), it has garnered extensive attention across various disciplines such as physics, biology, and chemistry [2,3,4,5,6,7,8,9]. Early investigations indicated that the essential prerequisites for inducing SR involve a nonlinear system, a weak periodic signal, and the presence of additive noise. Consequently, earlier research predominantly concentrated on exploring SR behavior in nonlinear systems driven by additive noise [10,11,12,13]. However, subsequent studies have revealed that SR phenomena can occur in linear systems driven by multiplicative noise [14,15,16,17]. Here, the preference for linear systems in subsequent research mainly stems from their ease of solution and widespread applicability [15].

A coupled system refers to a system consisting of finite (or infinite) local (or global) coupled elements [18], such as complex networks, where the particles are interconnected, making the coupled model crucial for accurately describing real-world problems. Consequently, there is a growing interest among scholars in exploring SR within coupled systems [19,20,21,22,23,24,25]. For instance, Vishwamittar et al. [19] delved into the study of collective resonant behavior in two coupled fractional oscillators subjected to quadratic asymmetric dichotomous noise. Yu et al. [20] proposed a model featuring two coupled fractional harmonic oscillators with fluctuating masses and investigated the SR phenomena through a combination of numerical and theoretical solutions. Lin et al. [21] explored the collective resonant behavior of coupled oscillators with fluctuating masses and incorporated a tempered Mittag–Leffler memory kernel. The coupling mechanism employed in this study involves linear interaction between the two oscillators, akin to the relationship between any two particles in a globally coupled system [22,23] or the connection between the two nearest particles in a nearest-neighbor coupling system [24,25].

In the context of evolving complex networks, the role of time delays becomes pivotal in the exploration of dynamic phenomena. Time delays are omnipresent in nature due to limitations in signal propagation, response times, and switching speeds [26]. Notably, Guillouzi et al. [27,28] introduced the small delay approximation, which is a novel method utilized for the first time to address stochastic differential equations featuring delay terms. Zhong et al. [29] extended the exploration by proposing a time-delayed fractional oscillator with fluctuation damping and investigating the generalized SR phenomena. In a related vein, Li et al. [30] observed that SR can be either enhanced or suppressed by manipulating the delay time and feedback strength. Concurrently, Jin et al. [31] investigated coherence resonance and SR in a time-delayed bistable system driven by both additive and multiplicative white noises. In addition, fractional models [32,33,34] are more suitable for describing memory effects and non-local behavior, and some numerical methods [35,36,37] have been proposed for solving fractional integro–differential equations. Therefore, generalized SR in time-delay fractional oscillators with random masses [38,39], random damping [40,41], and random frequencies [42] has been recently investigated through theoretical analysis and numerical calculations.

Furthermore, it is worth noting that all the aforementioned studies focused on the time-delayed stochastic differential model of individual particles, neglecting the consideration of interactions or coupling between particles. In reality, both time-delay and inter-particle coupling are prevalent in natural systems, leading to the proposition and exploration of numerous coupled models incorporating time delays. For instance, Kim et al. [43] delved into the dynamic behavior of coupled oscillators with time-delayed interaction: specifically under the influence of a pinning force. In a related study, Banerjee et al. [44] introduced the synchronization phenomenon in coupled hyperchaotic electronic oscillators featuring time delays. Additionally, Zhang et al. [45] investigated the SR phenomena exhibited in coupled systems with time delays while considering the impacts of mass and frequency fluctuations and further applied these findings to the diagnosis of bearing faults. Meanwhile, the influence of time delays on SR and discharge frequency oscillations within cortical neural networks is thoroughly analyzed in [46].

Nevertheless, limited attention has been given in the existing literature to the combined impact of time delays and frequency fluctuations in coupled fractional systems. Therefore, in this paper, we propose a novel system model featuring two coupled fractional oscillators with both time delay and frequency fluctuations and study the collective resonant behavior through a comprehensive exploration employing both theoretical and numerical methods. The primary outcomes of this investigation encompass the statistical synchronization observed between the two particles and the influences of coupling, time delay, and fractional order on the collective resonant behavior.

The subsequent sections are organized as follows: In Section 2, a detailed presentation of the system model outlining the characteristics of two coupled time-delayed fractional oscillators with frequency fluctuations is provided. Then, the analytical expression of the output amplitude gain is derived. Section 3 employs numerical results to scrutinize the collective resonant behavior in the proposed system. In Section 4, a numerical simulation is conducted to corroborate the analytical results. Finally, the succinct conclusions are summarized in Section 5.

2. System Model

We consider two coupled time-delayed fractional oscillators with frequency fluctuations that are described by the following fractional Langevin equations:

$$m{\ddot{x}}_1\left(t\right)+{\gamma }_0^C{D}_t^{\alpha }{x}_1\left(t\right)+\left[{\omega }^2+{\xi }_1\left(t\right)\right]x_1(t-\tau )=\epsilon \left[x_2\left(t\right)-x_1\left(t\right)\right]+Rcos(\Omega t)+{\zeta }_1\left(t\right),$$

(1)

$$m{\ddot{x}}_2\left(t\right)+{\gamma }_0^C{D}_t^{\alpha }{x}_2\left(t\right)+\left[{\omega }^2+{\xi }_2\left(t\right)\right]x_2(t-\tau )=\epsilon \left[x_1\left(t\right)-x_2\left(t\right)\right]+Rcos(\Omega t)+{\zeta }_2\left(t\right),$$

(2)

where $x_i\left(t\right)$ is the displacement of the ith particle at time t, $i=1,2$, and $\gamma >0$ is the damping coefficient, $\omega$ represents the system’s intrinsic frequency, $\tau >0$ is the time delay, and R and $\Omega$ represent the amplitude and frequency, respectively, of the periodic driving force. Moreover, $\epsilon$ is the coupling strength, and the two particles interact through the linear coupling term $\pm \epsilon \left(x_2-x_1\right)$. The viscous damping force related to the memory effect is defined by fractional derivatives with $\alpha$-order Caputo form:

$${\gamma }_0^C{D}_t^{\alpha }{x}_i\left(t\right)=\gamma {\int }_0^t\frac{1}{\Gamma (1-\alpha )}{|t-s|}^{-\alpha }{\dot{x}}_i\left(s\right)\mathrm{d}s,0<\alpha \le 1.$$

(3)

In the Equations (1) and (2), the internal noise terms ${\zeta }_1\left(t\right)$ and ${\zeta }_2\left(t\right)$ are fractional Gaussian noise and are supposed to satisfy the generalized second fluctuation–dissipation theorem [47]:

$$〈{\zeta }_i\left(t\right)〉=0,〈{\zeta }_i\left(t\right){\zeta }_i\left(s\right)〉=\frac{\gamma {\kappa }_{\mathrm{B}}T}{\Gamma (1-\alpha )}{|t-s|}^{-\alpha },i=1,2.$$

(4)

Here, the frequency fluctuations ${\xi }_1\left(t\right)$ and ${\xi }_2\left(t\right)$ are modeled as uncorrelated symmetric dichotomous noises taking two values ${\xi }_i\left(t\right)\in \{-\sigma ,\sigma \}$, $i=1,2,$, for which the statistical properties are represented as follows:

$$〈{\xi }_i\left(t\right)〉=0,〈{\xi }_i\left(t\right){\xi }_j\left(s\right)〉={\delta }_{ij}{\sigma }^2{\mathrm{e}}^{-\lambda |t-s|},i=1,2,$$

(5)

where $\sigma$ and $\lambda$ are the noise intensity and correlation rate, respectively. Further, ${\xi }_i\left(t\right)$ and ${\zeta }_i\left(t\right)$, $i=1,2$, are reasonably supposed to be uncorrelated due to their different origins; that is,

$$〈{\xi }_i\left(t\right){\zeta }_j\left(s\right)〉=0,i,j=1,2.$$

(6)

As mentioned in [20], the internal noises ${\zeta }_i\left(t\right)$ lose their effect on the first moments $〈x_i\left(t\right)〉$, $i=1,2$, based on the general assumptions (Equations (4) and (6)). In this study, we only focus on the first-order moment of the system’s steady response; therefore, ${\zeta }_1\left(t\right)$ and ${\zeta }_2\left(t\right)$ can be reasonably ignored for simplicity in form. In addition, we first perform the $O\left({\tau }^2\right)$ Taylor expansion around $\tau =0$ on the term $x(t-\tau )$ and obtain the following equivalent fractional Langevin equations without time delays:

$$m{\ddot{x}}_1\left(t\right)+{\gamma }_0^C{D}_t^{\alpha }{x}_1\left(t\right)+\left[{\omega }^2+{\xi }_1\left(t\right)\right]\left[x_1\left(t\right)-\tau {\dot{x}}_1\left(t\right)\right]=\epsilon \left[x_2\left(t\right)-x_1\left(t\right)\right]+Rcos(\Omega t),$$

(7)

$$m{\ddot{x}}_2\left(t\right)+{\gamma }_0^C{D}_t^{\alpha }{x}_2\left(t\right)+\left[{\omega }^2+{\xi }_2\left(t\right)\right]\left[x_2\left(t\right)-\tau {\dot{x}}_2\left(t\right)\right]=\epsilon \left[x_1\left(t\right)-x_2\left(t\right)\right]+Rcos(\Omega t).$$

(8)

It should be noted that the Equations (7) and (8) is an approximation of the Equations (1) and (2) only in the situation of a small time delay.

2.1. Complete Synchronization between the Average Behavior of the Two Oscillators

To analyze the collective resonant behavior of the coupled system (Equations (7) and (8)), it is imperative to assess whether the average behavior of the two particles is synchronized, which leads us to need to calculate $〈x_1\left(t\right)-x_2\left(t\right)〉$. In order to analyze the stationary state response of the Equations (7) and (8), we introduce the widely recognized Shapiro–Loginov formulas [48] and their generalized form [18] as preparation:

$$〈{\xi }_i\left(t\right)\frac{{\mathrm{d}}^n{x}_j\left(t\right)}{\mathrm{d}{t}^n}〉={\left(\frac{\mathrm{d}}{\mathrm{d}t}+\lambda \right)}^n〈{\xi }_i\left(t\right)x_j\left(t\right)〉,i,j=1,2,$$

(9)

$$〈{\xi }_1\left(t\right){\xi }_2\left(t\right)\frac{{\mathrm{d}}^n{x}_j\left(t\right)}{\mathrm{d}{t}^n}〉={\left(\frac{\mathrm{d}}{\mathrm{d}t}+2\lambda \right)}^n〈{\xi }_1\left(t\right){\xi }_2\left(t\right)x_j\left(t\right)〉,j=1,2,$$

(10)

and

$$〈{\xi }_i{\left(t\right)}_0^C{D}_t^{\alpha }{x}_j\left(t\right)〉={\mathrm{e}}^{-\lambda t}{ }_0^C{D}_t^{\alpha }\left(〈{\xi }_i\left(t\right)x_j\left(t\right)〉{\mathrm{e}}^{\lambda t}\right),i,j=1,2,$$

(11)

$$〈{\xi }_1\left(t\right){\xi }_2{\left(t\right)}_0^C{D}_t^{\alpha }{x}_j\left(t\right)〉={\mathrm{e}}^{-2\lambda t}{ }_0^C{D}_t^{\alpha }\left(〈{\xi }_1\left(t\right){\xi }_2\left(t\right)x_j\left(t\right)〉{\mathrm{e}}^{2\lambda t}\right),j=1,2,$$

(12)

which are important in the calculation process.

Firstly, by subtracting Equation (8) from Equation (7) we obtain

$$\begin{array}{cc}\hfill & \left(m\frac{{\mathrm{d}}^2}{\mathrm{d}{t}^2}+{\gamma }_0^C{D}_t^{\alpha }+{\omega }^2-\tau {\omega }^2\frac{\mathrm{d}}{\mathrm{d}t}+2\epsilon \right)\left[x_1\left(t\right)-x_2\left(t\right)\right]+\left[{\xi }_1\left(t\right)x_1\left(t\right)-{\xi }_2\left(t\right)x_2\left(t\right)\right]\hfill \\ \hfill & -\tau \left[{\xi }_1\left(t\right)\frac{\mathrm{d}}{\mathrm{d}t}{x}_1\left(t\right)-{\xi }_2\left(t\right)\frac{\mathrm{d}}{\mathrm{d}t}{x}_2\left(t\right)\right]=0.\hfill \end{array}$$

(13)

Then, we perform three operations on Equation (13): (I) averaging it with respect to the noise, (II) multiplying it by $\left({\xi }_1\left(t\right)+{\xi }_2\left(t\right)\right)$ and then averaging, and (III) multiplying it by $\left({\xi }_1\left(t\right){\xi }_2\left(t\right)\right)$ and then averaging. By using the Shapiro–Loginov formulas (9)–(12), we obtain

$$\left(m\frac{{\mathrm{d}}^2}{\mathrm{d}{t}^2}+{\gamma }_0^C{D}_t^{\alpha }+{\omega }^2-\tau {\omega }^2\frac{\mathrm{d}}{\mathrm{d}t}+2\epsilon \right)〈x_1\left(t\right)-x_2\left(t\right)〉+\left[1-\tau \left(\frac{\mathrm{d}}{\mathrm{d}t}+\lambda \right)\right]〈{\xi }_1\left(t\right)x_1\left(t\right)-{\xi }_2\left(t\right)x_2\left(t\right)〉=0,$$

(14)

$$\begin{array}{cc}\hfill & \left[m{\left(\frac{\mathrm{d}}{\mathrm{d}t}+\lambda \right)}^2+{\omega }^2-\tau {\omega }^2\left(\frac{\mathrm{d}}{\mathrm{d}t}+\lambda \right)+2\epsilon \right]\left[〈{\xi }_1\left(t\right)x_1\left(t\right)-{\xi }_2\left(t\right)x_2\left(t\right)〉+〈{\xi }_2\left(t\right)x_1\left(t\right)-{\xi }_1\left(t\right)x_2\left(t\right)〉\right]\hfill \\ \hfill & +\gamma {\mathrm{e}}^{-\lambda t}{ }_0^C{D}_t^{\alpha }\left[\left(〈{\xi }_1\left(t\right)x_1\left(t\right)-{\xi }_2\left(t\right)x_2\left(t\right)〉+〈{\xi }_2\left(t\right)x_1\left(t\right)-{\xi }_1\left(t\right)x_2\left(t\right)〉\right){\mathrm{e}}^{\lambda t}\right]\hfill \\ \hfill & +{\sigma }^2\left(1-\tau \frac{\mathrm{d}}{\mathrm{d}t}\right)〈x_1\left(t\right)-x_2\left(t\right)〉+\left[1-\tau \left(\frac{\mathrm{d}}{\mathrm{d}t}+2\lambda \right)\right]〈{\xi }_1\left(t\right){\xi }_2\left(t\right)\left(x_1\left(t\right)-x_2\left(t\right)\right)〉=0,\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill & \left[m{\left(\frac{\mathrm{d}}{\mathrm{d}t}+2\lambda \right)}^2+{\omega }^2-\tau {\omega }^2\left(\frac{\mathrm{d}}{\mathrm{d}t}+2\lambda \right)+2\epsilon \right]〈{\xi }_1\left(t\right){\xi }_2\left(t\right)\left(x_1\left(t\right)-x_2\left(t\right)\right)〉\hfill \\ \hfill & +\gamma {\mathrm{e}}^{-2\lambda t}{ }_0^C{D}_t^{\alpha }\left[〈{\xi }_1\left(t\right){\xi }_2\left(t\right)\left(x_1\left(t\right)-x_2\left(t\right)\right)〉{\mathrm{e}}^{2\lambda t}\right]+{\sigma }^2\left[1-\tau \left(\frac{\mathrm{d}}{\mathrm{d}t}+\lambda \right)\right]〈{\xi }_2\left(t\right)x_1\left(t\right)-{\xi }_1\left(t\right)x_2\left(t\right)〉=0.\hfill \end{array}$$

(16)

Next, we multiply the Equations (7) and (8) by ${\xi }_1\left(t\right)$ and ${\xi }_2\left(t\right)$, respectively, and then we subtract the two results. We further average the obtained equation, utilize the Shapiro–Loginov formulas (9)–(12) again, and have

$$\begin{array}{cc}\hfill & \left[m{\left(\frac{\mathrm{d}}{\mathrm{d}t}+\lambda \right)}^2+{\omega }^2-\tau {\omega }^2\left(\frac{\mathrm{d}}{\mathrm{d}t}+\lambda \right)+\epsilon \right]〈{\xi }_1\left(t\right)x_1\left(t\right)-{\xi }_2\left(t\right)x_2\left(t\right)〉+\epsilon 〈{\xi }_2\left(t\right)x_1\left(t\right)-{\xi }_1\left(t\right)x_2\left(t\right)〉\hfill \\ \hfill & +\gamma {\mathrm{e}}^{-\lambda t}{ }_0^C{D}_t^{\alpha }\left(〈{\xi }_1\left(t\right)x_1\left(t\right)-{\xi }_2\left(t\right)x_2\left(t\right)〉{\mathrm{e}}^{\lambda t}\right)+{\sigma }^2\left(1-\tau \frac{\mathrm{d}}{\mathrm{d}t}\right)〈x_1\left(t\right)-x_2\left(t\right)〉=0.\hfill \end{array}$$

(17)

Hence, we get closed Equations (14)–(17) with four newly defined variables: $y_1\triangleq 〈x_1\left(t\right)-x_2\left(t\right)〉$, $y_2\triangleq 〈{\xi }_1\left(t\right)x_1\left(t\right)-{\xi }_2\left(t\right)x_2\left(t\right)〉$, $y_3\triangleq 〈{\xi }_2\left(t\right)x_1\left(t\right)-{\xi }_1\left(t\right)x_2\left(t\right)〉$, $y_4\triangleq 〈{\xi }_1\left(t\right){\xi }_2\left(t\right)\left(x_1\left(t\right)-x_2\left(t\right)\right)〉$. In the long-time limit of $t\to \infty$, the influence of the initial conditions on the system’s stationary state response gradually disappears. Subsequently, we solve closed Equations (14)–(17) by using Laplace and inverse Laplace transforms, and we obtain the following important results:

$$\underset{t\to \infty }{lim}{y}_1=\underset{t\to \infty }{lim}{y}_2=\underset{t\to \infty }{lim}{y}_3=\underset{t\to \infty }{lim}{y}_4=0,$$

(18)

which can be rewritten as

$$\begin{array}{cc}\hfill & {\langle x_1\left(t\right)\rangle }_{\mathrm{as}}={\langle x_2\left(t\right)\rangle }_{\mathrm{as}}=\underset{t\to \infty }{lim}\langle x_1\left(t\right)\rangle ,\hfill \\ \hfill & {\langle {\xi }_1\left(t\right)x_1\left(t\right)\rangle }_{\mathrm{as}}={\langle {\xi }_2\left(t\right)x_2\left(t\right)\rangle }_{\mathrm{as}}=\underset{t\to \infty }{lim}\langle {\xi }_1\left(t\right)x_1\left(t\right)\rangle ,\hfill \\ \hfill & {\langle {\xi }_2\left(t\right)x_1\left(t\right)\rangle }_{\mathrm{as}}={\langle {\xi }_1\left(t\right)x_2\left(t\right)\rangle }_{\mathrm{as}}=\underset{t\to \infty }{lim}\langle {\xi }_2\left(t\right)x_1\left(t\right)\rangle ,\hfill \\ \hfill & {\langle {\xi }_1\left(t\right){\xi }_2\left(t\right)x_1\left(t\right)\rangle }_{\mathrm{as}}={\langle {\xi }_1\left(t\right){\xi }_2\left(t\right)x_2\left(t\right)\rangle }_{\mathrm{as}}=\underset{t\to \infty }{lim}\langle {\xi }_1\left(t\right){\xi }_2\left(t\right)x_1\left(t\right)\rangle .\hfill \end{array}$$

(19)

It can be concluded by Equation (19) that the average of the mean field is equal to the average of any single particle’s displacement in the long-time limit, i.e., ${\langle (x_1\left(t\right)+x_2\left(t\right))/2\rangle }_{\mathrm{as}}={\langle x_1\left(t\right)\rangle }_{\mathrm{as}}={\langle x_2\left(t\right)\rangle }_{\mathrm{as}}$, which shows that the average behavior of the two particles is completely synchronized. In addition, the results ensure the rationality and validity to study the mean field through $\langle x_1\left(t\right)\rangle$. Therefore, we only analyze one particle’s stationary state response $\langle x_1\left(t\right)\rangle$ in the following subsection.

2.2. Output Amplitude Gain Of System

We perform four operations on Equation (7): (I) averaging it with respect to the noise, (II) multiplying it by ${\xi }_1\left(t\right)$ and then averaging, (III) multiplying it by ${\xi }_2\left(t\right)$ and then averaging, and (IV) multiplying it by ${\xi }_1\left(t\right){\xi }_2\left(t\right)$ and then averaging. Then, we utilize the Shapiro–Loginov formulas (9)–(12) and the results in Equations (19) again, and thus we obtain the closed equations with four new variables as follows:

$$\left(m\frac{{\mathrm{d}}^2}{\mathrm{d}{t}^2}+{\gamma }_0^C{D}_t^{\alpha }-\tau \omega \frac{\mathrm{d}}{\mathrm{d}t}+{\omega }^2\right)z_1\left(t\right)+\left[1-\tau \left(\frac{\mathrm{d}}{\mathrm{d}t}+\lambda \right)\right]z_2\left(t\right)=Rcos(\Omega t),$$

(20)

$$\begin{array}{cc}\hfill & \left[m{\left(\frac{\mathrm{d}}{\mathrm{d}t}+\lambda \right)}^2-\tau {\omega }^2\left(\frac{\mathrm{d}}{\mathrm{d}t}+\lambda \right)+{\omega }^2+\epsilon \right]z_2\left(t\right)+\gamma {\mathrm{e}}^{-\lambda t}{ }_0^C{D}_t^{\alpha }\left(z_2\left(t\right){\mathrm{e}}^{\lambda t}\right)\hfill \\ \hfill & +{\sigma }^2\left(1-\tau \frac{\mathrm{d}}{\mathrm{d}t}\right)z_1\left(t\right)-\epsilon z_3\left(t\right)=0,\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill & \left[m{\left(\frac{\mathrm{d}}{\mathrm{d}t}+\lambda \right)}^2-\tau {\omega }^2\left(\frac{\mathrm{d}}{\mathrm{d}t}+\lambda \right)+{\omega }^2+\epsilon \right]z_3\left(t\right)+\gamma {\mathrm{e}}^{-\lambda t}{ }_0^C{D}_t^{\alpha }\left(z_3\left(t\right){\mathrm{e}}^{\lambda t}\right)\hfill \\ \hfill & -\epsilon z_2\left(t\right)+\left[1-\tau \left(\frac{\mathrm{d}}{\mathrm{d}t}+2\lambda \right)\right]z_4\left(t\right)=0,\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill & \left[m{\left(\frac{\mathrm{d}}{\mathrm{d}t}+2\lambda \right)}^2-\tau {\omega }^2\left(\frac{\mathrm{d}}{\mathrm{d}t}+2\lambda \right)+{\omega }^2\right]z_4\left(t\right)+\sigma \left[1-\tau \left(\frac{\mathrm{d}}{\mathrm{d}t}+\lambda \right)\right]z_3\left(t\right)\hfill \\ \hfill & +\gamma {\mathrm{e}}^{-2\lambda t}{ }_0^C{D}_t^{\alpha }\left(z_4\left(t\right){\mathrm{e}}^{2\lambda t}\right)=0,\hfill \end{array}$$

(23)

where $z_1\left(t\right)\triangleq 〈x_1\left(t\right)〉$, $z_2\left(t\right)\triangleq 〈{\xi }_1\left(t\right)x_1\left(t\right)〉$, $z_3\left(t\right)\triangleq 〈{\xi }_2\left(t\right)x_1\left(t\right)〉$, $z_4\left(t\right)\triangleq 〈{\xi }_1\left(t\right){\xi }_2\left(t\right)x_1\left(t\right)〉$.

Solving closed Equations (20)–(23) by the Laplace transform, we obtain

$$\left(\begin{array}{cccc}{a}_{11}& a_{12}& 0& 0\\ a_{21}& a_{22}& a_{23}& 0\\ 0& a_{32}& a_{33}& a_{34}\\ 0& 0& a_{43}& a_{44}\end{array}\right)\left(\begin{array}{c}{Z}_1\\ Z_2\\ Z_3\\ Z_4\end{array}\right)=\left(\begin{array}{c}\frac{Rs}{s^2+{\Omega }^2}\\ 0\\ 0\\ 0\end{array}\right),$$

(24)

where the elements $a_{ij}$ are given by

$$\begin{array}{cccc}\hfill & a_{11}=m{s}^2+\gamma s^{\alpha }-\tau {\omega }^{2}s+{\omega }^2,\hfill & \hfill a_{22}& =m{\left(s+\lambda \right)}^2-\tau {\omega }^2\left(s+\lambda \right)+\gamma {\left(s+\lambda \right)}^{\alpha }+{\omega }^2+\epsilon ,\hfill \\ \hfill & a_{12}=1-\tau \left(s+\lambda \right),\hfill & \hfill a_{33}& =m{\left(s+\lambda \right)}^2-\tau {\omega }^2\left(s+\lambda \right)+\gamma {\left(s+\lambda \right)}^{\alpha }+{\omega }^2+\epsilon ,\hfill \\ \hfill & a_{21}={\sigma }^2-\tau {\sigma }^{2}s,\hfill & \hfill a_{32}& =-\epsilon ,\hfill \\ \hfill & a_{23}=-\epsilon ,\hfill & \hfill a_{34}& =1-\tau \left(s+2\lambda \right),\hfill \\ \hfill & a_{43}={\sigma }^2\left[1-\tau \left(s+\lambda \right)\right],\hfill & \hfill a_{44}& =m{\left(s+2\lambda \right)}^2-\tau {\omega }^2\left(s+2\lambda \right)+\gamma {\left(s+2\lambda \right)}^{\alpha }+{\omega }^2,\hfill \end{array}$$

and $Z_i\left(s\right)=\mathcal{L}\left\{z_i\left(t\right)\right\}={\int }_0^{+\infty }{z}_i\left(t\right){\mathrm{e}}^{-st}\mathrm{d}t$ is the Laplace transform of $z_i\left(t\right)$, $i=1,2,3,4$.

To derive the steady-state response of the system, we focus only on the solution for $Z_1\left(s\right)$, which can be express by

$$Z_1\left(s\right)=H_1\left(s\right)\frac{Rs}{s^2+{\Omega }^2},$$

(25)

where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & H_1\left(s\right)=\frac{g\left(s\right)}{f\left(s\right)},\hfill \\ \hfill & g\left(s\right)=a_{22}{a}_{33}{a}_{44}-a_{22}{a}_{34}{a}_{43}-a_{23}{a}_{32}{a}_{44},\hfill \\ \hfill & f\left(s\right)=a_{11}{a}_{22}{a}_{33}{a}_{44}-a_{11}{a}_{22}{a}_{34}{a}_{43}-a_{11}{a}_{23}{a}_{32}{a}_{44}-a_{12}{a}_{21}{a}_{33}{a}_{44}+a_{12}{a}_{21}{a}_{34}{a}_{43}.\hfill \end{array}\end{array}$$

Applying the inverse Laplace transform to Equation (25), we get the analytical expression of $〈x_1\left(t\right)〉$ as follows:

$$〈x_1\left(t\right)〉=R{\int }_0^t{h}_1(t-t^{\prime })cos(\Omega t^{\prime })\mathrm{d}{t}^{\prime },$$

(26)

where $h_1\left(t\right)$ is the inverse Laplace transform of $H_1\left(s\right)$. Similarly, in a long-time limit $t\to \infty$, the effect of the initial conditions on the system’s stationary state response ${〈x\left(t\right)〉}_{\mathrm{as}}$ vanishes, and thus, ${〈x\left(t\right)〉}_{\mathrm{as}}$ can be asymptotically expressed as

$${〈x\left(t\right)〉}_{\mathrm{as}}={〈x_1\left(t\right)〉}_{\mathrm{as}}=A_{\mathrm{as}}cos(\Omega t+{\phi }_{\mathrm{as}}),$$

(27)

where $A_{\mathrm{as}}$ and ${\phi }_{\mathrm{as}}$ are the amplitude and phase shift, respectively, of ${〈x\left(t\right)〉}_{\mathrm{as}}$. More specifically, $A_{\mathrm{as}}$ and ${\phi }_{\mathrm{as}}$ can be expressed as

$$A_{\mathrm{as}}=R\left|H_1(j\Omega )\right|=R\sqrt{\frac{u_1^2+u_2^2}{u_3^2+u_4^2}},$$

(28)

and

$${\phi }_{\mathrm{as}}=arg\left(H_1(j\Omega )\right)=arctan\left(\frac{u_2{u}_3-u_1{u}_4}{u_1{u}_3+u_2{u}_4}\right),$$

(29)

where j is the imaginary unit, and the related coefficients are arranged in Appendix A.

Finally, the output amplitude gain is given by

$$G=\frac{A_{\mathrm{as}}}{R}=\sqrt{\frac{u_1^2+u_2^2}{u_3^2+u_4^2}}.$$

(30)

3. Collective Resonant Behavior

This section is dedicated to analyzing the collective resonant behavior of the two coupled time-delayed fractional oscillators with fluctuating frequencies in terms of the analytical expression of output amplitude gain G (30). Specifically, we mainly focus on investigating the synergistic effects of the coupling strength $\epsilon$, time delay $\tau$, fractional order $\alpha$, and noise parameters (${\sigma }^2$, $\lambda$) on the collective resonant behavior, including collective parameter-induced SR and collective conventional SR behaviors. Here, the parameter-induced SR phenomenon is presented by changing system parameters with fixed noise [49], and the conventional SR phenomenon is presented by adjusting the noise intensity ${\sigma }^2$ and correlation rate $\lambda$ with fixed system parameters.

3.1. Collective Parameter-Induced SR based on System Parameters

Firstly, we present the phase diagram in the $\sigma$–$\alpha$ plane in Figure 1 to analyze the emergence of collective parameter-induced SR in $G\left(\epsilon \right)$ both with and without a time delay (i.e., $\tau =0,0.05$). The gray region indicates that single-peak collective parameter-induced SR takes place, while the white region signifies there is no collective resonant behavior. Figure 1a shows that collective parameter-induced SR only occurs in the strong-memory region ($\alpha <0.45$, except for $\alpha =0$). A smaller $\alpha$ with stronger memory can induce the collective parameter-induced SR phenomena. Figure 1b indicates that when a time delay exists in the system ($\tau =0.05$), the resonant region evidently gets larger, suggesting that the time delay also can induce collective parameter-induced SR. Actually, both small $\alpha$ values and large $\tau$ values reflect strong memory effects to control the emergence of collective parameter-induced SR. Without loss of generality, we take $\tau =0.05$ as an example and further analyze the effects of $\alpha$ and $\sigma$ on $G\left(\epsilon \right)$. As shown in Figure 1c, an increase in $\alpha$ results in a decrease in the peak value of $G\left(\epsilon \right)$, accompanied by an obvious leftward shift in the peak position. Conversely, as shown in Figure 1d, the resonance peak of $G\left(\epsilon \right)$ is observed to rise with an increase in noise intensity $\sigma$, accompanied by a slightly leftward shift in the peak position.

To further explore the effect of the small time delay $\tau$ on $G\left(\epsilon \right)$, we plot the nephogram of $G(\epsilon ,\tau )$ and the curves depicting $G\left(\epsilon \right)$ for various values of $\tau$. As depicted in Figure 2, G varies non-monotonously with increasing $\epsilon$, indicating the occurrence of the collective parameter-induced SR phenomenon with one peak. It is evident in Figure 2b that the peak value experiences a sharp increase with increasing $\tau$, and the peak position shifts to the right. It follows then that increasing $\tau$ and decreasing $\alpha$, reflecting strong system memory, can increase the maxima of $G\left(\epsilon \right)$ (see Figure 1c and Figure 2b). Furthermore, another observation from the figures reveals that in the weakly coupled region, G exhibits a non-monotonic variation with an increase in $\epsilon$, while in the strongly coupled region, the curves converge to a limit value. This phenomenon can be interpreted as follows: when $\epsilon$ reaches a sufficiently large value ${\epsilon }^{*}$, the coupling forces $\pm \epsilon (x_1-x_2)$ become relatively large, compelling the two particles to move cohesively. Consequently, any further increase in $\epsilon$ does not impact the motion of the two particles, resulting in a constant value of G rather than a change, indicating that it stabilizes at a fixed value.

In order to have a deeper insight on the influence of $\tau$ on G, we plot Figure 3, in which we show the curves of $G\left(\tau \right)$ for various $\epsilon$ and $\alpha$ values. Notably, all the curves exhibit a discernible peak, signifying that the collective parameter-induced SR appears. As shown in Figure 3a, there is complex change in $G\left(\tau \right)$. As $\epsilon$ increase, $G\left(\tau \right)$ rises first and falls rapidly later. After arriving at the bottom, $G\left(\tau \right)$ creeps up again and tends to a fixed limit (also see Figure 2b). Moreover, the peak position of $G\left(\tau \right)$ moves to the larger $\tau$, meaning stronger system memory. As shown in Figure 3b, as $\alpha$ increase, the peak value of $G\left(\tau \right)$ decreases, with the peak position moving to the larger value of $\tau$. We can conclude that there is a smallest system memory with a combination of a fractional order $\alpha$ and time delay $\tau$ that can induce optimal matching between the system and noise when other system parameters are fixed. Thus, the peak position ${\tau }_{cr}$ moves to the right as $\alpha$ increase. That is, the weakening of memory brought by increasing $\alpha$ offsets the enhancement of memory caused by increasing $\tau$.

Similarly, the impact of $\alpha$ on G is investigated in Figure 4, in which we present the curves of $G\left(\alpha \right)$ for various $\epsilon$ and $\tau$ values. The non-monotonous curves of $G\left(\alpha \right)$ mean that collective parameter-induced SR occurs. Comparing with Figure 3a, $G\left(\alpha \right)$ exhibits a similar variation trend as $\epsilon$ increases in Figure 4a, but the peak position of $G\left(\alpha \right)$ moves slightly to the smaller $\alpha$, which also means stronger system memory. As shown in Figure 4b, the peak value of $G\left(\alpha \right)$ first increases and then decreases with increasing $\tau$, and the peak position moves to the larger value of $\alpha$ significantly. Here, the enhancement of memory brought by increasing $\tau$ also offsets the decay to memory caused by increasing $\alpha$.

3.2. Collective Conventional SR based on Noise Parameters

Finally, we delve into the dependence of the output amplitude gain G on the noise intensity $\sigma$ and correlation rate $\lambda$. Figure 5 and Figure 6, respectively, illustrates the curves of $G\left(\sigma \right)$ and $G\left(\lambda \right)$ for varying values $\epsilon$, $\tau$, and $\alpha$; the figures show non-monotonous behavior signifying the occurrence of collective conventional SR phenomena.

As shown in Figure 5a, the peak value of $G\left(\sigma \right)$ varies non-monotonically as $\epsilon$ increases: specifically, increasing first and then decreasing until it remains unchanged. There is also an interesting phenomenon in the peak position in that it moves to the right initially and then moves to the left. We can see from Figure 5b that with an increase in $\tau$, the peak value of $G\left(\sigma \right)$ increases, and the peak position shifts to the larger value of $\sigma$. Furthermore, the similar influence of $\alpha$ on $G\left(\sigma \right)$ is shown in Figure 5c. As system memory is enhanced, as reflected by decreasing $\alpha$, the peak value of $G\left(\sigma \right)$ increases quickly, and the peak position shifts to the larger value of $\sigma$ slightly. Obviously, both increasing $\tau$ and decreasing $\alpha$ can enhance the intensity of collective conventional SR. The phenomena can be attributed to the fact that the increase in system memory is beneficial for transferring noise energy to periodic signals.

In Figure 6a, it is evident that the peak value of $G\left(\lambda \right)$ rises with an increase in $\epsilon$, and the peak position shifts to the left. This result signifies that an appropriately chosen coupling strength $\epsilon$ can enhance the collective resonant intensity. As in Figure 6b, the presence of a time delay leads to the occurrence of collective conventional SR with one peak and one valley. As $\tau$ increases, $G\left(\lambda \right)$ increases and reaches a maximum at a larger $\lambda$. In Figure 6c, there is no resonance behavior observed for $\alpha \ge 0.6$. However, for $\alpha <0.6$, a decrease in $\alpha$ leads to a increase in the peak value accompanied by a rightward shift in the peak position. The phenomena observed in Figure 6b,c can be elucidated as follows. Time delay $\tau$ and fractional order $\alpha$ reflect the memory characteristics of the system, and their effects are opposite. Meanwhile, the noise correlation rate $\lambda$ reflects the memory characteristics of the noise. Specifically, a larger $\lambda$ represents weaker memory of the noise. The system’s memory effect resulting from an increase in $\tau$ or a decrease in $\alpha$ necessitates a larger noise correlation rate to counteract it. In addition, it can be concluded from Figure 6 that $\epsilon$, $\tau$, and $\alpha$ play a significant role in controlling the collective conventional SR intensity of $G\left(\lambda \right)$.

4. Numerical Simulation

In this section, we present numerical simulations to validate the credibility and precision of the analytical results. Utilizing the Grünwald–Letnikov definition of a fractional derivative [50], the numerical simulations are performed for the two coupled fractional oscillators described by the fractional Langevin Equations (1) and (2).

Firstly, the Grünwald–Letnikov definition of a fractional derivative on integral interval $[a,t]$ is expressed as follows:

$${ }_a^G{D}_t^{\alpha }x\left(t\right)=\underset{h\to 0^{+}}{lim}\frac{1}{h^{\alpha }}\sum _{r=0}^{\left[\frac{t-a}{h}\right]}{\left(-1\right)}^r\left(\genfrac{}{}{0pt}{}{\alpha }{r}\right)x\left(t-rh\right),$$

(31)

where $\alpha \in (0,1]$ is the fractional order. Secondly, the fractional differential equation ${ }_a^G{D}_t^{\alpha }x\left(t\right)=f\left(x\left(t\right),t\right)$ can be discretized as follows:

$${ }_{t_1}^G{D}_{t_k}^{\alpha }x\left(t\right)=f\left(x\left(t_k\right),t_k\right),k=1,2,\cdots ,K.$$

(32)

where K is the total number of sampling points, $h=\frac{t-a}{K}$ is the time step of the calculation, and $t_k=kh$, $k=1,2,\cdots ,K$. Then, Equation (32) can be rewritten as

$$\frac{1}{h^{\alpha }}\sum _{r=0}^{k-1}{\left(-1\right)}^r\left(\genfrac{}{}{0pt}{}{\alpha }{r}\right)x\left(t_{k-r}\right)=f\left(x\left(t_k\right),t_k\right),k=1,2,\cdots ,K,$$

(33)

and we obtain

$$x\left(t_k\right)=h^{\alpha }f\left(x\left(t_k\right),t_k\right)-\sum _{r=1}^{k-1}{C}_r^{\left(\alpha \right)}x\left(t_{k-r}\right),k=1,2,\cdots ,K.$$

(34)

Substituting $f\left(x\left(t_k\right),t_k\right)$ with $f\left(x\left(t_{k-1}\right),t_{k-1}\right)$, we obtain the recurrence formula for $x\left(t_k\right)$:

$$x\left(t_k\right)=h^{\alpha }f\left(x\left(t_{k-1}\right),t_{k-1}\right)-\sum _{r=1}^{k-1}{C}_r^{\left(\alpha \right)}x\left(t_{k-r}\right),k=1,2,\cdots ,K,$$

(35)

where the binomial coefficients $C_r^{\left(\alpha \right)}$ are expressed as follows:

$$C_0^{\left(\alpha \right)}=1,C_r^{\left(\alpha \right)}=\left(1-\frac{1+\alpha }{r}\right)C_{r-1}^{\left(\alpha \right)}.$$

(36)

Next, we reformulate the fractional Langevin Equations (1) and (2) to obtain the numerical solution of the two coupled fractional oscillators:

$$\left\{\begin{array}{c}\frac{{\mathrm{d}}^{\alpha }{x}_1\left(t\right)}{\mathrm{d}{t}^{\alpha }}=v_1\left(t\right),\hfill \\ \frac{{\mathrm{d}}^{\alpha }{x}_2\left(t\right)}{\mathrm{d}{t}^{\alpha }}=v_2\left(t\right),\hfill \\ m\frac{{\mathrm{d}}^{2-\alpha }{x}_1\left(t\right)}{\mathrm{d}{t}^{2-\alpha }}=-\gamma v_1\left(t\right)-\left[{\omega }^2+{\xi }_1\left(t\right)\right]x_1(t-\tau )+\epsilon \left[x_2\left(t\right)-x_1\left(t\right)\right]+Rcos(\Omega t),\hfill \\ m\frac{{\mathrm{d}}^{2-\alpha }{x}_2\left(t\right)}{\mathrm{d}{t}^{2-\alpha }}=-\gamma v_2\left(t\right)-\left[{\omega }^2+{\xi }_2\left(t\right)\right]x_2(t-\tau )+\epsilon \left[x_1\left(t\right)-x_2\left(t\right)\right]+Rcos(\Omega t).\hfill \end{array}\right.$$

(37)

Hence, the general discrete forms of the equations can be expressed as

$$\left\{\begin{array}{c}{x}_1\left(t_k\right)=h^{\alpha }{v}_1\left(t_{k-1}\right)-\sum _{r=1}^{k-1}{C}_r^{\left(\alpha \right)}{x}_1\left(t_{k-r}\right),\hfill \\ x_2\left(t_k\right)=h^{\alpha }{v}_2\left(t_{k-1}\right)-\sum _{r=1}^{k-1}{C}_r^{\left(\alpha \right)}{x}_2\left(t_{k-r}\right),\hfill \\ v_1\left(t_k\right)=h^{2-\alpha }\left[-\gamma v_1\left(t_{k-1}\right)-\left({\omega }^2+{\xi }_1\left(t_{k-1}\right)\right)x_1(t_{k-1}-\tau )+\epsilon \left[x_2\left(t_{k-1}\right)-x_1\left(t_{k-1}\right)\right]\right]/m\hfill \\ +h^{2-\alpha }\left[Rcos(\Omega t_{k-1})\right]/m-\sum _{r=1}^{m-1}{C}_r^{\left(2-\alpha \right)}{v}_1\left(t_{k-r}\right),\hfill \\ v_2\left(t_k\right)=h^{2-\alpha }\left[-\gamma v_2\left(t_{k-1}\right)-\left({\omega }^2+{\xi }_2\left(t_{k-1}\right)\right)x_2(t_{k-1}-\tau )+\epsilon \left[x_1\left(t_{k-1}\right)-x_2\left(t_{k-1}\right)\right]\right]/m\hfill \\ +h^{2-\alpha }\left[Rcos(\Omega t_{k-1})\right]/m-\sum _{r=1}^{k-1}{C}_r^{\left(2-\alpha \right)}{v}_2\left(t_{k-r}\right).\hfill \end{array}\right.$$

(38)

Furthermore, we are concerned with the average trajectory of particle i, which can be approximated by the Monte Carlo method as

$$〈x_i\left(t\right)〉=\frac{1}{N}\sum _{n=1}^N{x}_i^{\left(n\right)}\left(t\right),i=1,2,$$

(39)

where N is the number of realizations, and $x_i^{\left(n\right)}\left(t\right)$ signifies the nth implementation of particle i. For this purpose, we set the Monte Carlo experiment time $N=100$ and choose the simulation time $T=100$ s with a time step $h=0.01$ s. By setting the system parameters with $\omega =1$, $\Omega =0.6\pi$, $\gamma =0.5$, $\alpha =0.8$, $\tau =0.1$, $\epsilon =1$, $\lambda =0.5$, and $\sigma =0.5$, the absolute errors at these selected points with different mass values m are shown in

Table 1. The absolute error between analytical and numerical results of $G\left(m\right)$.

m	0.5	0.7	0.9	1.0	1.1	1.3	1.5
$\left|G_N-G_A\right|$	0.0062	0.0030	0.0006	0.0017	0.0025	0.0000	0.0004

, which shows that the numerical and analytical results maintain consistency well. In the following numerical simulations, we set $m=1$ for simplicity.

Figure 7a demonstrates the time domain waveform of symmetric dichotomous noises ${\xi }_1\left(t\right)$ and ${\xi }_2\left(t\right)$. Figure 7b shows the motion trajectories of the two particles over an extended period, showcasing perfect synchronization of their long-term average displacement. Figure 7c presents a comparison between the analytical and numerical results of $〈x\left(t\right)〉$, confirming that the simulation result is in good agreement with the analytical result. To further verify the reliability and accuracy of the analytical result in Equation (30), we depict the curves of G versus $\epsilon$, $\alpha$, and $\tau$ in Figure 8, juxtaposing numerical results with the theoretical values of G. Consequently, we can safely conclude from Figure 8 that within an acceptable margin of error, the numerical results align well with the analytical results for small time delays. Notably, Figure 8c reveals a high degree of consistency between theoretical and numerical results for $\tau \le 0.3$, which indicates again that the analytical solution is valid for the two coupled fractional oscillators for a small time delay.

Finally, we provide a simple case to reveal the positive effect of the proposed system on weak signal detection. It is assumed that two coupled particles in the system are respectively driven by weak periodic signal $f\left(t\right)=Rcos\left(2ft\right)$ (with $R=0.1$, $f=0.2$) and strong noise ${\zeta }_1\left(t\right)$ and ${\zeta }_2\left(t\right)$ (with $D=4$). Figure 9a shows the time-domain waveform of weak signal $f\left(t\right)$, and the corresponding spectrum is depicted in Figure 9b, where a weak frequency feature can be clearly identified at the driving frequency $f=0.2$ Hz. However, in the strong noise background, the weak signal is completely submerged in the time domain, as shown in Figure 9c,e. Moreover, the weak frequency feature cannot be identified in the spectra, as shown in Figure 9d,f. By regulating system parameters $m=1$, $\omega =0.8$, $\gamma =1$, $\alpha =0.4$, $\tau =0.1$, $\epsilon =1$, $\lambda =0.5$, and $\sigma =0.5$, we solve for system responses and obtain the mean field. It is observed from Figure 9g,h that the interferences are effectively suppressed in the time domain, and the component at $f=0.2$ Hz can also be clearly identified. It reflects the positive effect of the proposed system on amplifying weak signals. Actually, this phenomenon has received widespread attention in mechanical fault diagnosis.

5. Conclusions

In this study, we propose two coupled time-delayed fractional oscillators with fluctuating frequencies to explore collective resonant behavior, with a specific focus on the influence of the coupling strength $\epsilon$, time delay $\tau$, and fractional order $\alpha$. Employing the stochastic average method, we obtain complete synchronization between the average behaviors of the two particles and derive the analytical expression for the output amplitude gain G.

Based on the analytical results, we observe the collective parameter-induced SR phenomena in $G\left(\epsilon \right)$, $G\left(\tau \right)$, and $G\left(\alpha \right)$. More concretely, time delay $\tau$ is conducive to the resonant region expansion of $G\left(\epsilon \right)$. Obviously, increasing $\tau$ and decreasing $\alpha$ can enhance the intensity of the collective parameter-induced SR in $G\left(\epsilon \right)$. In addition, $\epsilon$ has a similar complex influence on the collective parameter-induced SR in $G\left(\tau \right)$ and $G\left(\alpha \right)$. We also observe that collective conventional SR phenomena take place in the coupled system. It is worth emphasizing that $\epsilon$, $\alpha$, and $\tau$ play a significant role in controlling the collective conventional SR behavior of $G\left(\sigma \right)$ and $G\left(\lambda \right)$ by facilitating the emergence and optimizing the intensity of collective resonant behavior. Subsequently, we corroborate the accuracy and validity of the analytical results through numerical simulation. Finally, we anticipate that the results of this study can furnish theoretical support for subsequent investigations for some potential applications.