Revealing fractionality in the Rössler system by recurrence quantification analysis

This study discusses the results of a Recurrence quantification analysis (RQA) of the Rössler system with a fractional order (q1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_1$$\end{document}) of the derivative in the first equation. The fractional order q1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_1$$\end{document} changes slightly in the range q1∈⟨0.9,1.0⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_1 \in \langle 0.9,1.0\rangle$$\end{document}. Even with such relatively small changes in the q1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_1$$\end{document} derivative, significant changes in the dynamics of the system are observed between the bifurcation diagrams determined for the bifurcation parameter a. Nevertheless, as q1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_1$$\end{document} decreases one can notice the preservation of some structures of the bifurcation diagram, in particular the main periodic windows of the integer-order Rössler system. The RQA shows clear differences between various regular windows of the integer system and only slight changes in these windows are caused by an increase in the system’s fractionality. Nonetheless, by selecting appropriate recurrence variables it is possible to expose the changes occurring in the regular windows under the influence of the fractionality of the system. This approach allows for the detection of the fractional character of the system through a recurrence analysis of the time series taken from periodic regions.


Introduction
Systems of differential equations with a fractional derivative order (fractional systems) are increasingly used to model real phenomena and technical processes [1][2][3][4][5][6][7][8][9][10][11][12]. They describe the evolution of dynamical systems for which the current state depends not only on the current state but also on previous states (i.e. on the history of the system), which is clearly visible in their mathematical definitions [13]. Since the end of the 20th century fractional order derivatives have also been used more and more in modelling highly complex phenomena [5]. The mathematical apparatus of fractional derivatives plays an important role in solving many engineering problems [1][2][3][4]6]. Quite often it is used to model and control magnetic systems [7][8][9], and it has also been successfully implemented in completely different disciplines, such as biology [3], medicine [11], and physiology [10]. The fractional models are applied to describe systems and objects of various scales, from the nanoscale [12], up to cosmic scales [14].
Methods for detecting nonlinearities and chaos in integral and fractional systems are constantly being developed and tested [15,16]. In the case of exploration and/or control of such systems, the ability to a e-mail: a.rysak@pollub.pl (corresponding author) detect the fractional nature of the system, identify fractional terms, and evaluate the order of fractional derivatives may be very important. System identification is crucial for steering problems. The introduction of various definitions of fractional derivatives and methods of solving them was followed by the development in the identification of fractional systems [17]. These early works assumed a given order of the fractional derivative and were focused on finding the system parameters using the standard least square method [18,19]. Further works extended the scope of identification of the order of the derivative of the fractional system [20,21]. Some identification methods cope very well even with noisy time series, where the proposed scheme offers a high degree of accuracy for erroneous data [20]. Similarly like in systems with integer order derivatives, the identification of the fractional system is performed using methods working in the time domain [22,23], and in the frequency domain [2,[24][25][26]. Recently further refinements to methods of identification of fractional order systems have been developed [27][28][29][30]. In the presented work we show a numerical method of detection of fractionality in the Rössler system, which is based on the recurrence analysis of a periodically oscillating system. The proposed approach does not call for complex mathematical algorithms. It is based on the statistics of the recurrence plot dots, depicting the return of the system to the vicinity of its earlier states. The advantage of this approach is its computational simplicity, which 84 Eur. Phys. J. Spec. Top. (2023) 232:83-98 can be performed on relatively short time series [31]. This method also gives good results in the analysis of complex time series [32], signals with noise [33], nonstationary [34] and transients [35]. Recurrence quantification analysis (RQA) is a method which has already been used effectively to study changes in the dynamics of various systems, using a variety of measures, procedures, and algorithms [36][37][38][39]. Despite this, studies of recurrences in fractional systems are hard to find. Of the three found [40][41][42], two were written by the authors of the current study Considering such a wide and effective application of RQA in the study of various dynamical systems, this method was chosen as a potential tool to evaluate the fractionality of dynamical systems.
The current numerical tests check whether the results of recurrence analysis can show changes depending on the system fractionality large enough to enable the evaluation of the fractional character of a dynamical system. The Rössler system [43,44] is chosen as the research subject, with weak nonlinearity appearing only in the third equation. In this study, an additional slight fractionality is adopted in the first equation of the system. Based on the changes in the bifurcation diagram, resulting from the change in the system fractionality, two windows of the bifurcation parameter a were selected for further recurrence analysis. These windows correspond to intervals of a where the behaviour is regular. The RQA results change slightly due to slight changes in the weak fractionality of the system under study. In chaotic and quasi-periodic regions, recurrence variables are difficult to separate due to relatively high level of noise. But in the periodic modes, the noise level decreases. The obtained results showed that in some periodic intervals of the regular windows it is possible to find a RQA subspace in which the changes caused by the (relatively weak) fractionality of the system are clearly visible. This confirms that the recurrence quantification analysis can be considered as a potential tool for the detection and evaluation of fractionalities in dynamical systems.
The paper is organized as follows. Section 2 gives the equations defining the fractional Rössler system. Section 3 describes the numerical methods and techniques used for solving the fractional differential equations (Sect. 3.1), computing the recurrence variables (Sect. 3.2), and determining the bifurcation diagrams (Sect. 3.3). The results of the analysis of the bifurcation diagrams, determined for different system fractionalities, are presented in Sect. 4.1. Recurrence results are discussed in Sects. 4.2 and 4.3. The part of the study devoted to the discussion of numerical results ends with Sect. 4.3, in which the presented approach was applied to time series with different noise levels. The last section concludes and summarises.

Rössler system
The Rössler system was developed in 1976 by Otto Rössler [43], and is the simplest three-dimensional chaotic system, with nonlinearity only in the third equation. In the generalization with fractional order derivatives, the system is defined as follows [44,45]: where we adopt b = 2, c = 4, and the initial values at t = 0 equal to x(0) = 0, y(0) = 1, z(0) = 9. The operator D α t means the fractional derivative of order α. Its meaning and the method of numerical calculation are given in Sect. 3.1. q n , (n = 1, 2, 3) are in general the fractional orders of derivatives with values q 2 = q 3 = 1 and q 1 ∈< 0.9, 1.0 > adopted in the current study. The bifurcation parameter a is swept in a wide interval with bounds depending on the q 1 value.
Numerical solutions of the fractional system (1) for a given a value were obtained by implementation of the Grünwald-Letnikov method (Sect. 3.1), solutions of the integer system (q 1 = 1) were obtained using a simple Euler algorithm with a small time step. The following parameters affect the results of all simulations: the integration step h, the settling time ST and the next period of time AT from which the data were collected for further numerical analysis. For the integer system h = 2 · 10 −4 , ST = 5000, and AT = 3000. In the simulations of the fractional system, it was assumed: h = 10 −3 , ST = 10, 000, and AT = 300.

Fractional derivative
In the numerical calculations of this work, the Grünwald-Letnikov method [46,47] was adopted to solve the fractional differential equation. Grünwald-Letnikov's approach is a direct extension of the derivative of integer order n to that of a fractional order q and is given by the formula Discretization of the independent variable t in time limits 0, t n (t n = nh) gives the approximation The binomial coefficients c which simplifies further numerical calculations. Applying the 'short memory' principle [49], transforms Eq.
(3) into where L = mh is a short time interval from which the values of previous states are taken to compute the solution at time t n , p = n for t n < L and p = m for t n ≥ L.
Assuming the general form of the fractional equation the solution at the n-th time step is given by In calculating the solutions of fractional systems, the system history length L = 17, 500h was adopted.

Recurrence analysis
Recurrence plot analysis was initiated in 1987 by Eckmann et al. [50], with the idea of representing the dynamics of a time-varying system using a graphical chart called a Recurrence Plot (RP ). A recurrence plot is a square matrix of dots with two axes representing time. The black dots (t i , t j ) of the RP correspond to the moments of time in which detects the proximity of system states V(t k ) in the phase space. ε is the threshold parameter and the i, j, k = 1, 2, . . . , N indices number the successive values of the time t for which the values of the vector V of system states were recorded. The recurrence quantification analysis is an extension of RP analysis by developing statistical variables measuring distribution of dots in the RP . This idea was developed by Zbilut and Webber [51,52] and was further developed by Casdagli [53], Marwan [32,54] and others. In this study, we examine the three-dimensional Rössler system, described by the full set of three differential equations, so the time delay embedding [55] is not required here. Before the construction of the RP , all system variables are normalised according to the following scheme: x N = (x − x)/σ x , where x is a chosen system variable and σ x is the standard deviation of the variable x . In the presented studies of the system the following recurrence variables were implemented to test the dynamics of the Rössler system: where are the probability of occurrence of a diagonal line of length l and a vertical line of length v , respectively. H D and H V are histograms of diagonal (D) and vertical (V ) line length distribution. In the above formulas, the main diagonal of RP is consistently omitted, because it represents the closeness of the system states with themselves and therefore does not add any valuable recurrence information. RR (recurrence rate) is the dots density of the RP ; D is the average length of the diagonal lines and TT (trapping time) is the average length of the vertical lines. LAM (laminarity) is defined as the ratio of the number of dots on vertical lines to all recurrence dots. LEnt and VEnt are the Shannon entropy of the diagonal and vertical line length distributions, respectively. Definitions and interpretations of recurrence measures (9)-(14) are given in detail in [56]. As shown in [57] for the Duffing system and in [58] for the Rössler system, recurrence variables also strongly depend on the density of time series from which the RP is built. In this study, the time series taken for the RQA analysis were AT = 300 long and were diluted to a density of 20 points per unit of t, giving 6000 points used to build recurrence plots.

Bifurcation diagrams
The bifurcation diagram gives the general characteristics of a dynamical system in which one of the parameters, called the bifurcation parameter, changes. In our study, the a parameter was selected for this role. To build a bifurcation diagram, one must determine the puncture points of the Poincaré plane (x = 0) through the time series which is the solution of the system (1) for a given a, in the time interval of NS seconds after the settlement time (ST ). The bifurcation diagram presents the projections of Poincaré points on the Y axis for all values of the bifurcation parameter a. Each bifurcation diagram was computed for 1920 values of the parameter a. For the entire range of the bifurcation parameter, its values changed with the h a step, the value of which was selected taking into account the range of changes of the bifurcation parameter and the system characteristics.

Bifurcation diagrams of fractional system
The main inspiration for this research were the modifications of the bifurcation diagram caused by the increasing fractionality of the system (due to the decreasing of q 1 ). The selected results obtained for decreasing value of the q 1 parameter are presented in Fig. 1.
The step Δa for the bifurcation diagrams q 1 = 1.00, q 1 = 0.97 and q 1 = 0.94 was 1.15 · 10 −04 , and for the bifurcation diagram q 1 = 0.92 it was reduced to 0.9 · 10 −04 (to compensate for the diagram compression effect due to the increase in the fractionality). Increasing fractionality in the first equation shifts the bifurcation diagram towards higher values of the bifurcation parameter a. It is interesting that as the fractionality of the system increases, the windows R2 and R3 (marked in Fig. 1) maintain approximately the same position in relation to the diagram, while many other regular windows appear and disappear during this process. The scale of the a axis in the subsequent diagrams in Fig. 1 is changed in such a way that the regular windows R2 and R3 are aligned vertically, thus indicating stability To determine the changes in the position and size of the bifurcation diagram following the decrease of q 1 value, two reference points were selected: BP-the point of the first bifurcation marked by a circle in the diagram obtained for q 1 = 1 (Fig. 1) and EoR3-the end of interval of the regular window R3.
As can be seen in Fig. 2a, decreasing the value of q 1 shifts both reference points towards higher values of the bifurcation parameter. At the same time, the distance between these characteristic points is decreasing (Fig. 2b). This shows that the bifurcation diagram of the system, when decreasing the order of the derivative of the first equation, shifts towards the higher values of the bifurcation parameter and becomes compressed.
Despite the deformation of the bifurcation diagram with respect to the bifurcation parameter, the result in Fig. 1 shows that wide regular windows do not change their bifurcation characteristics and relative position in the diagram. Taking into account previous studies [58] which showed that recurrences of regular windows are very sensitive to changes in the dynamics of the system, two intervals corresponding to two regular windows: R2 and R3, marked in Fig. 1, were selected for further recurrence tests. Further details of bifurcation diagrams are given in Appendix A.

Regular windows of integer-order system
When a system parameter is changed, its dynamics change, which is reflected in the results of the recurrence analysis. It is shown in Fig. 3a for the R3 window.
In this case, the parameter a is scanned in the range a ∈< 0.5000, 0.5121 >, which shifts the dynamics of the system through the R3 window. The colors of the markers reflect the values of the a parameter. The system at   the beginning is chaotic (dark blue markers in the upper right corner), next it jumps sharply to the periodic mode (lower left corner), and returns to chaos (yellow markers) through the period doubling cascade and the quasi-periodic region. Figure 3b compares trajectories of R2 and R3 regular windows in the {RR, D ,LEnt} RQA subspace. The results obtained for both regular windows are clearly separated in this space. Thus the RQA can potentially be used to identify different regular windows of the Rössler system.

RQA of the fractional system in regular windows
In order to check whether RQA can detect differences in system fractionality, the calculations of the RQA variables were performed for both selected regular windows (R2,R3), according to the formulas given in Sect.
3.2. Figure 4 compares the results obtained in the {RR, D ,LEnt} RQA subspace for integer and fractional systems in R2 (a) and R3 (b) regular windows for low fractionality defined by q 1 = 0.99. Under these conditions, a slight difference can be noticed in the results obtained for the regular window R2 (Fig. 4a), while the results of the integer and fractional systems obtained for the R3 window almost coincide (Fig. 4b). This corresponds to the results of bifurcation, which in the range of higher values of the parameter a are more resistant to the fractional nature of the system.
By increasing the fractionality of the system, a slight shift and modification in the individual recurrence results are observed. This is shown in Fig. 5  As can be noted, the largest differences in the results caused by the change in the fractionality of the system are visible for the variable V Ent (Fig. 5). Similar analyses performed for the R2 window gave qualitatively similar results-the biggest changes in the recurrence results due to the change in the system fractionality occur for the variable V Ent (Fig. 6). The areas marked on Figs. 5c and 6 correspond to periodic modes in which the RQA variables stabilize, thus they become better separated and better identify solutions that differ in fractionality.
The differences that are clearly visible in the results of V Ent obtained for different values of q 1 occur in a limited range of changes in the bifurcation parameter a. In search of location of these intervals, the V Ent results obtained in the periodic area of the R3 window for individual q 1 values were plotted as a function of the parameter a and compared with the bifurcation diagram (Fig. 7).
A detailed analysis has shown that the largest differences in V Ent occur in a narrow window of the bifurcation parameter (Δa), which is located just after the first bifurcation point of the period doubling cascade. This area is marked with vertical lines on Fig. 7a for q 1 = 0, 99 and for q 1 = 0, 96 on Fig. 7b. The regular window R3 was selected for further analysis, because its bifurcation characteristics are stable in the whole range of changes in the fractionality of the system which are examined in the current study. In order to find the range of the biggest sensitivity of the variable V Ent to changes in the fractional structure of the system, additional calculations were performed in the periodic range of the R3 window (Sect. 4.4). The question why the V Ent measure works effectively in this case is discussed in Appendix B.

Revealing fractionality by the R3 periodic window
The preliminary results showed that the method of determining the Δa q1 intervals (for individual q 1 values) on the basis of the results presented in Fig. 7 is not accurate enough. To obtain a satisfactory separation of results in the subspace {RR, V Ent}, the intervals Δa q1 should be determined based on the area of the highest differences in {RR, V Ent}. Such areas are boxed in Figs. 5c and 6. This approach is demonstrated on the example of the series of results obtained for the extended set of fractional order values in the R3 window, for q 1 ∈ {.90, 0.92, 0.94, 0.96, 0.97, 0.98, 0.99, 1} and for ε = 0.06 (Fig. 8a). The value of ε was selected as a non-optimal value to demonstrate that sharp optimisation of the threshold parameter is not required, see Appendix C for more details.
As in the case of Figs. 5c and 6, the separation area is marked with a black frame. For each result, the frame limits the ΔRR range in which it is best separated  from the rest of the results in the {RR, V Ent} subspace. As can be noticed in Fig. 8a, only a small fraction of the points are in the marked area. Here too, we intentionally want to show the effectiveness of the procedure when using non-optimal initial data. In this case it means that we can start the analysis from the scan of the bifurcation parameter through the whole regular window. Based on the separation area, ΔRR q1 intervals are determined, corresponding to the results obtained for the individual q 1 values. The ΔRR q1 intervals are then converted to Δa q1 intervals using the RR(a) relationships. Figure 8b illustrates this transformation for q 1 = 0.98. In the next step, the calculations are repeated for all q 1 values in their respective Δa q1 intervals. For each q 1 , 128 time series were calculated for 128 a values, uniformly covering the respective Δa q1 interval. Then, the full RQA was performed for these points, which gave the results presented in Fig. 9a, in the {RR, V Ent} RQA subspace on which the current study is focused.
The frame drawn in Fig. 9a has the same coordinates as the frame in Fig. 8a. The mismatch of some results to the framed area is the result of a (intentionally) nonoptimal selection of the initial data (Fig. 8a). Figure 9b, presents points obtained by averaging of the RR and V Ent variables for each result indexed by q 1 . As one can see, positions of these points in the {RR, V Ent} RQA subspace correlate with the fractionality of the system. The dependence of the averaged V Ent values on q 1 (Fig. 10) confirms the effectiveness of the proposed approach in detecting the weak fractionality of the Rössler system without noise. The effects of noise are detailed in Appendix C.  Fig. 9 a The results of the V Ent variable obtained for individual values of q1 in narrow Δa intervals, which were determined using the procedure illustrated in Fig. 8a,b.

Summary
In this study, bifurcation diagrams and recurrence analyses were used to study the fractional Rössler system. It was assumed that the system fractionality is relatively weak and is caused by the fractional order q 1 of the derivative occurring only in the first equation of the system and varying in the range q 1 ∈ 0.9, 1.0 . The bifurcation diagram of the integer system determined for the bifurcation parameter a changes as the system fractionality changes. As q 1 decreases, the bifurcation diagram shifts towards higher values of a, undergoes compression and local modifications. However, when comparing the diagrams with respect to selected specific reference points, it can be seen that some large regular windows are not affected by weak fractionality of the system. The recurrence results calculated for these regular areas show clear differences between the various regular windows and change slightly in response to slight changes in the system fractionality.
The use of the proposed approach as a tool for detecting and assessing the fractionality of the Rössler system requires the predefinition of appropriate values for the main parameters of the analysis, i.e. for the density of time series points and the threshold parameter ε. The ε value should be chosen depending on the noise level of the analysed signal, based on appropriate calibration. As the detailed analysis of the R3 regular window has shown, an important element of the procedure is to determine the appropriate range of the bifurcation parameter value (Δa), in which the recurrence results show the largest changes under the influence of the change in the Rössler system fractionality.
In conclusion, as a result of numerical tests of both selected regular windows, recurrence variables were found which in the dynamic region correspond to the period doubling cascade and which show clear differences between systems that slightly differ in fractionality. In this way, we showed the potential possibilities of using recurrence analysis for detection and evaluation of the weak fractionalities in the Rössler dynamical system.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and 92 Eur. Phys. J. Spec. Top. (2023) 232:83-98 your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Data availibility statement Data availability statement: no data associated in the manuscript.

A Details of the bifurcation diagram modified by weak fractionality
The change in the Rossler system fractionality implies many interesting changes in its bifurcation diagram, determined for the bifurcation parameter a. Some of them, of particular interest for the analysis under consideration, are discussed below.
For example, the R2 window is poorly detectable for q 1 = 0.92, but it is an example of an interesting phenomenon in the bifurcation representation of changes in system dynamics. Figure 11 presents bifurcation diagrams made for the R2 window range, for q 1 = 1 (Fig. 11a), q 1 = 0.96 (Fig. 11b) and q 1 = 0.92 (Fig. 11c).
In cases (a) and (b), the window R2 is immersed in a chaotic surroundings. However, for q 1 = 0.92, the R2 window gives the impression of being superimposed on changes in dynamics taking place in its surroundings. It looks like the 'fixed' R2 window is superimposed on the background dynamics of the system. The 'background' dynamics change with increasing parameter a from chaotic, through quasi-periodic (invisible), period halving cascade ('thin branches'-partially visible), to periodic. Finally, there is a sudden return of system dynamics to a chaotic mode. It seems that the original regular windows result from the intrinsic properties of the system and it is difficult to eliminate them by changing the dynamics of the system by modifying its parameters.
In Fig. 1 there is another wide window, located for a ≈ 4.1 (for q 1 = 1). Let's call this window R1. It can be seen from Fig. 1 that the window R1 disappears for q 1 > 0.94, whereas the R2 and R3 windows still exist at this level of the system fractionality. However, with a further increase in the fractionality of the system (by reducing q 1 ), the arrangement of the two fixed windows is disturbed. The R3 regular window is intact in the results obtained for q 1 = 0.9 (Fig. 12), while the R2 window is clearly modified in this case by strong changes in bifurcations.
These results (especially Fig. 12) show that the increase in the fractional structure of the system causes strong changes invading the bifurcation diagram, especially from the side of low values of parameter a.

B Suitability of V Ent for analysis of the periodic mode
It is obvious that for the periodic mode and a small ε value, the statistics of vertical lines are very simple. To answer the question stated above, this section is dedicated to discuss the specific properties of the V Ent a) measure using the results obtained in Sect. 4.4 (Fig. 9a).
Since we compare here the RQA results computed from single time series obtained for different q 1 , for each q 1 one result lying in the middle of its corresponding Δa q1 interval is chosen. When an RP is built for a periodic system, its dots form a pattern of thin diagonal lines. When the tested time series is periodic and complex, gaps and local thickenings may appear on some diagonal lines of the RP . In Fig. 13, the exemplary RP is presented, calculated for the time series obtained from the fractional system (q 1 = 0.96) in the time interval Δt = 25. When a small value of the threshold parameter is used, the local thickenings of the diagonal lines are only a few dots wide. As a result, the distribution of the length of the vertical lines of such a layout is very simple and the statistics of these lines do not contain much information. Therefore, the special role played by V Ent in the present study may be surprising, although the results of studies of various dynamical systems were previously reported in which the interesting characteristics of this variable were also explored [58]. In order to check the reasons for the effectiveness of the V Ent measure in the assessment of some phenomena, the results of the a) RQA were analysed, calculated for the data selected from the sets discussed in the previous section (shown in Fig. 9a). The next two graphs (Fig. 14) show the length histograms of the diagonal (a) and vertical (b) lines, calculated for two cases: q1 = 1 and q1 = 0.98. The calculations were made for Δt = 300. Diagonal line histograms (Fig. 14a) are very complex and show complex differences between the two results being compared. Such characteristics of the diagonal lines histogram mean that when changing the q 1 parameter, the LEnt values have a large dispersion and a very weak correlation with the q 1 values. Unlike LEnt, the histograms of both V Ent scores are very simple (Fig. 14b). Extending this result to all tested fractionalities of the system gives the set of histograms presented in Fig. 15.
It shows that the histograms of vertical lines, although they are very simple, change systematically with the change of the value of q 1 . This result is transferred to the statistical measure of V Ent, the values of which also change systematically with the change of q 1 .
The statistical results confirm the effective operation of the V Ent measure in distinguishing Rössler systems that differ slightly in fractionality. To interpret the efficiency of this variable in the context of the phase space and basics of RQA, the phase portraits of three systems differing in the q 1 value are compared in Fig. 16.
Interpretation can be derived from two features of phase portraits. First, all results presented in Fig. 16 are obtained for the value of the bifurcation parameter slightly higher than the value corresponding to the first bifurcation point (Fig. 7). Which means that the trajectories in this figure have a double loop structure, and local distance between these loops varies in both direction and value. Second, phase portrait obtained for a given q 1 shows turning points near which the density of trajectory points decreases. As a result, when the RP is constructed, the number of points crossing the ε ball changes markedly at the turning points, which is further enhanced by the slight separation of points lying on different trajectory loops. Changes in q 1 cause changes in the trajectory of the system, which entails a change in both discussed features of the phase portrait. This implies changes in the pulsation of the thickness of the diagonal lines, resulting in a different value of the V Ent measure.
This interpretation gives rise to the most important conditions that must be taken into account in the presented analysis. First of all, the density of the points in the time domain must be appropriate to the complexity of the periodic trajectory. Time series should be obtained from a narrow periodic area, located slightly after the first bifurcation point. Finally, the ε parameter should be selected so as to obtain the greatest possible difference in the statistics of vertical lines for changes in trajectories caused by weak changes in the fractionality of the system. Figure 17 shows the V Ent(RR) results obtained for the periodic range of the R3 window, for different system fractionalities and for three values of the threshold parameter.  Fig. 17 The results of V Ent obtained for a sweep of a through the R3 windows for different q1 and ε values: a ε = 0.05, b ε = 0.06, and c) ε = 0.07

C Importance of the ε parameter
In the Fig. 17a-c, the V Ent axis has the same range, which makes it possible to notice the decreasing separation of the results with increasing ε. To show that the proposed approach gives valuable results for a range of ε (i.e. it does not require sharp optimisation of the threshold parameter), a non-optimal value ε = 0.06 was adopted in further calculations.

D Influence of noise
The method presented above was also applied to the Rössler system, where noise was introduced to each variable v ∈ {x, y, z}. Before the recurrence analysis each of the variables was normalised according to the formula v N = v−v σv . Then Gaussian noise of order 5%, 10% and 15% was added, which resulted in signals with an SNR of 26, 20, and 16.5 dB, respectively.
Due to the assumptions limiting the analysis to periodic modes, its application in the real world is limited to stationary signals. Possible value shifts and scale differences of individual variables are corrected by normalizing their values. The addition of noise brings the analysed series very close to real periodic signals, but many features of real (measured) time series were not included in this study. The potential adaptation of the proposed approach to the analysis of real signals requires further research and preferably confronting the presented numerical method with real oscillators described by both integer and fractional equations.
The RQA results of the signals with noise are presented in Fig. 18 in a manner analogous to the results obtained for the time series without noise (Fig. 9).
On the left hand side the results for the subspace {V Ent, RR} are shown for (a) SN R = 26dB, (b) SN R = 20dB, and (c) SN R = 16.5dB. On the right hand side the results for the corresponding average values of V Ent and RR can be seen: (d), (e), and (f), respectively. The introduction of noise into the investigated time series necessitates higher values of the threshold parameter when calculating the entries of an RP matrix. In the analysis of the signal with noise the values of were obtained in simple comparative tests. No systematic investigation towards optimisation of ε for different levels of noise and fractionality were conducted. Based on the comparative tests, an appropriate value of the threshold parameter ε was established for each of the three different noise levels. The layout of the points in subspace {RR, V Ent} obtained by averaging ( Fig. 18d-f) depends on the adopted ε value. Hence, the criterion used in the comparative tests was based on finding the ε value for which the points are arranged in an orderly and monotonically variable manner, best showing the differences between the results obtained for different q 1 values. In this way, the following threshold parameter values were determined for the individual noise levels: (a) ε = 0.10 for SN R = 26dB, (b) ε = 0.20 for SN R = 20dB, and (c) ε = 0.27 for SN R = 16.5dB. The results for the average values show good data point separation in the recurrence subspace, and that they are correlated with corresponding values of q 1 is clearly visible. However, when a noise level is added to the system vibrating in periodic mode, the diagonal lines of its RP become disturbed, which change the statistics of the diagonal and vertical lines and results in modification of the RQA variables. This entails a change in the {RR, V Ent} results, which are used in the current study to assess the fractionality of the system. As can be seen in Fig. 18d-f, an increase in the noise level reduces the differences between the V Ent coordinates of the averaged points, and at the same time increases the differences between their RR coordinates. Thus, as the noise level increases, the result of the separation depends less and less on the V Ent variable and more and more on the RR variable. Unfortunately, the obtained results are not sufficient to make a thesis that the increase in the noise level will always be reflected in a systematic shift of the average value of the RR measure. This issue requires further detailed analysis.