# Scaling and fine structure of superstable periodic orbits in the logistic map

Frédéric Perrier, Frédéric Girault

► To cite this version:

Frédéric Perrier, Frédéric Girault. Scaling and fine structure of superstable periodic orbits in the logistic map. Chaos, Solitons and Fractals: the interdisciplinary journal of Nonlinear Science, and Nonequilibrium and Complex Phenomena, 2022, 165, 10.1016/j.chaos.2022.112767. insu-04462292

HAL Id: insu-04462292

https://insu.hal.science/insu-04462292v1

Submitted on 25 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1

2

3

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

# Scaling and fine structure of superstable

## periodic orbits in the logistic map

Frédéric Perrier* and Frédéric Girault

Université Paris Cité, Institut de Physique du Globe de Paris, CNRS, F-75005 Paris, France

Revised Version 3.0, September 2022, submitted to *Chaos, Solitons and Fractals*

## ABSTRACT

The logistic map, defined by the recurrence $x_{n+1}=\alpha x_n(1-x_n)$, has been a major tool in the emergence of chaos theory. When the parameter $\alpha$ is varied from 0 to 4, various bifurcations occur, with cycles of period 2, 4, 8 and so on, appearing successively, till cycles of period $2^k$ are all unstable and chaos emerges. Even after chaos has emerged, stable periodic cycles of period $p$ can appear in some parameter $\alpha$ windows, which include one superstable cycle, characterized by a zero first derivative. In this paper, all superstable values of $\alpha$ are obtained numerically for period $p$ up to 21, and the two smallest superstable values are established for $p$ values up to 1000, plus $p=1024$ and $p=2048$. These superstable parameters exhibit discrete fine structures, often different for even and odd $p$, combined with global and local scaling properties. The number of superstable orbits increases like $2^{p-1}/p$ for large values of $p$, with a distribution of values approximately following a gamma distribution with mean $\alpha_M=4-0.0753+0.0026p$. Our solution for the largest superstable value of $\alpha$ at a period $p$, written $\alpha=4-\epsilon_p$, follows Lorenz scaling prediction $\epsilon_p=\pi^2/4p$, to an accuracy better than $10^{-3}\%$. The

<footnote>* Corresponding author</footnote>

<footnote>Email addresses: perrier@ipgp.fr (F. Perrier), girault@ipgp.fr (F. Girault)</footnote>
22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

next superstable values $\alpha=4-\beta_{p,j}$ follow the $(1+2j)^2$ multiplier to $\epsilon_p$ predicted by Simó and Tatjer up to values of $j$ increasing with $p$. The width of the superstable orbits, the $x_n$ domain around initial point $1/2$ with periodicity $p$, determined numerically, obeys a scaling law $\sqrt{\beta}/2^p$, for $\beta$ sufficiently small. For $p$ larger than 15, the periodicity of the calculated time-series sometimes displays numerical instabilities. Thus, the numerical logistic map differs from the mathematical ideal, which invites for caution, but it also provides a generic imperfect model closer to reality and experimental data when investigating the role of higher order periodic orbits as underlying structures of chaotic behaviours.

Word count: 300 (maximum 300 words).

Keywords: logistic map, chaos, periodic orbits, superstability, scale invariance

**Highlights:**

* Periodic orbits are key steps in the onset and structure of chaotic time-series

* Superstable parameters calculated until period 21 and first two until period 1000

* Superstable parameters have scaling properties intertwined with fine structure

* Ultimate superstable orbits follow Lorenz and Simó-Tatjer scalings versus period

* Numerical instabilities of superstable time-series appear for period larger than 15
## 1. Introduction

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

The logistic function $f(x)=\alpha x(1-x)$ is popular since Verhulst showed in 1838 that it describes non-exponential population growth data [1], and this fact was illustrated again recently when, together with similar functions, it was used to describe COVID-19 epidemic data [e.g., 2,3 and references therein]. Apart from being useful in practical situations, the logistic equation emerged as an actor of one major scientific revolution as the discrete recurrence $x_{n+1}=\alpha x_n(1-x_n)$, where $\alpha$ is a real parameter. While non-linear recurrence equations were already considered in the second part of the 19th century by Schröder [4], the novel mathematics that was to impact the entire 20th century, including topological anomalies, fractal sets and functional equations, was discovered in the iteration of complex numbers by Julia in 1918 [5] and by Fatou in 1919 [6]. However, Lorenz, as he discovered chaos in low-dimensional differential systems, was also the pioneer who established the logistic map as the simplest non-linear map able to describe the emergence of irregular time-series from periodic cycles of increasing period [7], with a structure of bifurcations and properties of attractors that was revealed for the first time by May [8]. Deterministic chaos, discovered in the context of meteorology and fluids dynamics by Lorenz [7] or of population dynamics by May [8], created a paradigm shift similar to relativity and quantum mechanics, that invaded all fields of knowledge [9], in particular in the description of natural phenomena at human scale, for example in geosciences [10], hierarchical fracturation and earthquakes [11], geophysical and environmental processes [12] or geomorphology [13].

When the concepts of chaos and non-linear dynamics are applied to artificial and natural processes, for example earthquakes [14], one open question remains the relationship between scaling properties of parameters (i.e., Gutenberg-Richter magnitude law) and the predictability of temporal behaviours (i.e., earthquake cycles). Another most intriguing aspect, important to describe natural phenomena, is how chaotic behaviour emerges from the loss of stability of periodic orbits, and the logistic map, despite its apparent simplicity, is able to capture numerous aspects of this
complexity [7-10]. For $\alpha$ positive smaller than 1, the iteration converges to the origin. At $\alpha=1$, the
origin becomes unstable and a period 2 solution emerges, the 2 points of the cycle being solutions
of the equation $f(f(x))=f^{(2)}(x)=x$. At $\alpha=3$, the period 2 orbit becomes unstable and a cycle of period
4 emerges, then, when it becomes unstable (for $\alpha=1+\sqrt{6}$), a stable cycle of period 8 emerges and
so on. When the period increases, the stability window of the periodic cycles becomes smaller.
Then, when $\alpha$ is larger than about 3.569946, chaotic behaviour emerges, as revealed, for example,
by a property of Sensitivity to Initial Conditions (SIC) [8]. The way chaos emerges from period
doubling appears universal, namely independent of the particular function $f(x)$ considered, provided
some basic conditions are met (unimodal function with leading exponent larger than 1) [15].

The next intriguing property revealed by the logistic map is the fact that, even after the onset
of chaos, in some windows of the parameter $\alpha$, stable periodic orbits are still possible for periods
$p=3, 5, 6$ and larger. These periodic orbits of period $p$ are solutions of $f^{(p)}(x)=x$, and they are stable
when $|df^{(p)}(x)/dx|<1$. The existence of period 3 was actually shown by Li and Yorke to be related
to the emergence of chaos [16]. The periodic orbits appear in a definite order as a function of $\alpha$ that,
like the Feigenbaum scaling [15], seems to be universal [17,18]. In the windows of stable periodic
cycles, some particular values of $\alpha$ lead to superstable orbits, originally called “central” by Lorenz
[7], for which $|df^{(p)}(x)/dx|=0$. For $\alpha=4$, which is the maximum value that can be considered
keeping the iteration bounded, periodic solutions exist, in increasingly large numbers depending on
the integer $p$, but all these periodic cycles are unstable. For negative values of $\alpha$, no new dynamics
emerges; indeed, it was shown that a simple transformation leads to another logistic map with a
positive parameter [19]. Here, we therefore only consider $0 \le x_n \le 1$ and $1 \le \alpha \le 4$.

The topological properties of the stable cycles of quadratic recurrence relations, which include
the logistic map, and of their attraction basins, were established in details by mathematicians
between 1970 and 1990 [20-22], expanding the original result of Julia [5] stating that for any
parameter $\alpha$, at most one stable periodic cycle is possible. When no periodic cycle is stable, then
the behaviour is chaotic, and the sets of stable periodic and chaotic parameters have Lebesgue measures whose sum in absolute value is 4 [23]. While these mathematical insights are essential to understand the nature of chaotic phenomena, explicit solutions for the windows of periodic solutions of the logistic map between the onset of chaos and $\alpha=4$, important when considering practical applications, were only studied recently. Period 3, first, was extensively studied [24-28]: period 3 cycles appear for $\alpha$ between $1 + 2\sqrt{2} \approx 3.828427$ [24-26] and $\alpha \approx 3.841499$ [28], with a superstable orbit for $\alpha \approx 3.831874$, associated with a super-unstable orbit [28]. Exact analytical expressions are now known for these characteristic values [25,28]. While the parameter window for the period 3 stable cycle is unique, as for period 2, all higher order periods have multiple windows immersed in the approximate chaotic parameter interval $3.569946 < \alpha < 4$. While analytical solutions were given for period 5 and 7 [29-31], the periodic domains were studied numerically for $p$ larger or equal to 8 [32-35]. While some studies focused on the problem of the determination of rigorous bounds on the Lebesgue dimensions [32-33], complete numerical solutions for the periodic windows were obtained by Galias up to period 32 [34], then up to period 36 [35] and some first scaling relations versus period were also reported.

In this paper, as scale invariance plays a crucial role in the geophysical experimental data we are analysing, we focus on various scaling properties of the superstable parameters. Following the method proposed by Tucker and Wilczak [32], we use the fact that the high-order periodic cycles in the logistic map can be explored by focusing on the superstable periodic orbits, for which reliable numerical solutions are calculated from the roots of one polynomial. First, we recall the basic properties of the superstable orbits, and we present our numerical method. We explore the first two periodic orbits for $p \le 1000$, $p=1024$ and $p=2048$ and provide a complete set of period orbits for $p \le 21$. We use these results to study how various types of scaling behaviours are observed at global and local parameter scales, and we confirm that the superstable orbits with the largest values of $\alpha$ follow the scaling laws versus period $p$ proposed by Lorenz in 1964 [7] for the largest parameter and subsequently generalized for other largest values by Simó and Tatjer in 1991 [22]. In addition,
novel scaling laws are proposed for the width of the superstable orbits, also determined numerically.

For periods larger than $p=15$, we also observe that the numerical time-series for some superstable
parameters are unstable. The consequences of this limitation due to computer numerical precision
are discussed in the conclusion, together with potential perspectives of the scaling properties of
superstable orbits for the understanding of the various intermittent and chaotic regimes.

124

## 2. Basic properties of the logistic map

126

### 2.1 Stable periodic windows

128

The main properties of the logistic map were already established in the seventies [7,8]. In the following, we concentrate on values of $\alpha$ near the maximum 4, and therefore we define the variable $\beta=4-\alpha$. The distribution of iteration points, when shown with a logarithmic scale as a function of $\beta$ (Fig. 1), shows better the emergence of stable periodic windows, and the increasing amount of non-periodic domains for small values of $\beta$. The famous mechanism of period doubling [7-9,15] is seen from $\beta=1$, as the fixed point $x_1=1-1/\alpha=(3-\beta)/(4-\beta)$ becomes unstable for $\alpha \ge 3$. When chaotic regime starts, which occurs for $\beta$ smaller than about 0.430054, the iterated points $x_n$ are not able to converge to a discrete set of points, but are distributed between the minimum $x_{\min}$ and the maximum $x_{\max}$ possible values (Fig. 1), given by:

138

$$ \begin{cases} x_{\max} = \frac{\alpha}{4} = 1 - \frac{\beta}{4} \\ x_{\min} = f(x_{\max}) = \frac{\alpha^2(1-\alpha)}{4} = \frac{\beta}{4}(4-\beta)\left(1-\frac{\beta}{4}\right). \end{cases} \quad (1) $$

The distribution of points between $x_{\min}$ and $x_{\max}$ is non-uniform, a property which has been getting increasing attention lately [36-37]. After the initiation of the chaotic regime, which can be defined precisely by SIC [8], stable periodic orbits appear in definite intervals of $\beta$. The most conspicuous such interval is the period 3 window between $\beta \approx 0.1716$ and $\beta \approx 0.1585$ (Fig. 1), but

140

141

142
windows with $p=4, 7, 5$ and $6$ are clearly visible. The points of these low-order periodic orbits can be found numerically and shown as a function of $\beta$ (Supp. Mat. Fig. S1). When stable periodic orbits appear, superstable orbits are necessarily present [20,32].

146

147

## 2.2 Superstable periodic orbits

148

In the case of the logistic map, the condition for a superstable cycle with a set of points $x_i$, with $0 \le i \le p-1$, can be written [25,27]:

151

$$ \frac{df^{(p)}}{dx}(x_i) = \alpha^p \prod_{j=0}^{p-1} (1 - 2x_j) = 0. \qquad (2) $$

This leads to the simple consequence that one of the points in the set $x_i$ is equal to $1/2$, say $x_0$.

Consequently, the superstable orbits, in Fig. 1 (and in Supp. Mat. Fig. S1), appear when, within the window of stability, one of the branches crosses the value $x=1/2$. Another consequence is that the superstable orbits can be found directly without the prior knowledge of the periodic branches.

Indeed, it is sufficient to write that the set of points, starting with $x_0=1/2$, is closed [25,28,32]. This is the method we are using in the following.

158

Let us define, for period $p$, the superstability polynomial $P_p(\beta)$, defined by:

159

$$ P_p(\beta) = Q_p(\beta) - \frac{1}{2}, \qquad (3) $$

where $Q_k(\beta)$, with $1 \le k \le p$, is the series of polynomials in $\beta$ defined by:

161

$$ \left\{ \begin{array}{l} Q_1(\beta) = x_{\max} = 1 - \frac{\beta}{4} \\ Q_{k+1}(\beta) = (4 - \beta)Q_k(\beta)[1 - Q_k(\beta)] \end{array} \right. \qquad (4) $$

Superstability at period $p$ appears for the zeroes of the superstability polynomial $P_p(\beta)$, which has degree $D_p=2^p-1$. This degree increases fast: $D_1=1, D_2=3, D_3=7, \dots, D_{10}=1023$ and so on. Before discussing the general numerical solution for extracting the roots of the superstability conditions, we first discuss further general properties.

166
## 2.3 Properties of the superstability polynomials

The superstability condition can be displayed graphically (Fig. 2). Superstability for $p=1$ corresponds to $P_1(\beta)=Q_1(\beta)-1/2=(2-\beta)/4=0$, thus $\beta=2$, which is a trivial solution outside the range of Fig. 2. However, since 1 divides of all integer values of $p$, 2 is always a root of $P_p(\beta)$. More generally, if $\beta_p$ is a root of $P_p(\beta)$, then it will automatically be a root of all $P_q(\beta)$ where $q$ is a multiple of $p$.

The range of $P_p(\beta)$ is constrained (Fig. 2) between a minimum $P_{min}(\beta)$ and a maximum $P_{max}(\beta)$, following Eq. (1), given by:

$$ \left\{ \begin{array}{l} P_{\max}(\beta) = \frac{\alpha}{4} - \frac{1}{2} = \frac{1}{2} - \frac{\beta}{4} \\ P_{\min}(\beta) = \frac{\alpha^2(1-\alpha)}{4} - \frac{1}{2} = \frac{\beta}{4}(4-\beta)\left(1-\frac{\beta}{4}\right) - \frac{1}{2} \end{array} \right. \qquad (5) $$

When going from periodicity $p-1$ to $p$, there is a simple property that is worth pointing out. Let $\beta_{p-1}$ be a root of $P_{p-1}$, which means that $Q_{p-1}(\beta_{p-1})=1/2$. Then, using Eqs. (3) and (4), the value of $P_p$ at $\beta_{p-1}$ is given by:

$$ P_p(\beta_{p-1}) = (4 - \beta_{p-1})Q_{p-1}(\beta_{p-1})[1 - Q_{p-1}(\beta_{p-1})] - \frac{1}{2} = P_{\max}(\beta_{p-1}). \qquad (6) $$

Furthermore, the first derivative at that point is constrained:

$$ \begin{aligned} \frac{dP_p}{d\beta}(\beta_{p-1}) &= -Q_{p-1}(\beta_{p-1})[1 - Q_{p-1}(\beta_{p-1})] + (4 - \beta_{p-1})\frac{dQ_{p-1}}{d\beta}(\beta_{p-1})[1 - 2Q_{p-1}(\beta_{p-1})] \\ &= -\frac{1}{4} = \frac{dP_{\max}}{d\beta}(\beta_{p-1}). \end{aligned} \qquad (7) $$

Thus, for $\beta=\beta_{p-1}$, $P_p(\beta)$ is equal to the maximum value $P_{max}(\beta)$ and the two curves are tangent at that location. This property is illustrated in Fig. 2 for $p=8$ and $p=9$.

# 3. Numerical solution for the superstable orbits

## 3.1 Numerical method
190

191

192

193

194

195

196

The problem of finding the zeroes of a polynomial is hardly new. However, the particular nature of the superstability polynomial, namely the large degree and the large number of zeroes over a narrow parameter range, complicates the matter. The simple Newton method performed poorly in this case, with convergence problems. Indeed, not even half of the zeroes for small values of $\beta$ were found for $p$ larger than 3. Therefore, while a better Newton method can be used [32], a simpler but efficient algorithm was designed.

197

198

199

200

201

202

203

204

205

206

The algorithm, using real numbers in digital double precision [38], is an iterative zero crossing search by identifying intervals where the function gives opposite signs for the endpoints, and it proceeds as follows. A certain interval in $\ln\beta$ was divided into $N$ equal segments, and the segments $\beta_k$ to $\beta_{k+1}$ with $P_p(\beta_k)$ and $P_p(\beta_{k+1})$ of opposite sign are identified. Then, the process is repeated, dividing the identified segments with opposite signs into $M$ equal segments, then scanning these segments for opposite sign of the polynomial; this process is again repeated with $M$ segments. This way, the points of zero crossing of the function $P_p(\beta)$ are rapidly found with high precision. We considered values of $N$ and $M$ between $2^{20}$ and $2^{24}$, with $\alpha$ between 3.1 and 4 ($0 \le \beta \le 0.9$). The number of zeroes found in this range is noted $Z_p$. A number of tests were performed to check the reliability and the accuracy of this method for the considered situation.

207

208

## 3.2 Systematic tests

209

210

211

212

213

214

For $p=3$, an analytical solution for the unique superstable point is known, giving $\alpha_3=3.8318740552833155$ [28]. The solution given by our numerical algorithm differs from this value only at the level of the last (16th) significant given digit (Table 1), suggesting that the numerical precision for obtaining the zeroes of the superstability polynomial is at the level of $10^{-15}$, at least for small values of $p$.
215

216

217

218

219

220

221

222

For moderate values of *p*, the results of our algorithm were compared with the results given by the general Jenkins and Traub method (JT), as available in standard mathematical code libraries [39]. The results again coincide, but only for $p \le 5$. When the degree of the polynomial becomes larger than 30, the JT algorithm gives complex expressions for all $D_p$ zeroes of $P_p$. However, the precision given by the JT method for the real zeroes in our domain of concern, near $\alpha=4$, becomes poor, and the zeroes identified by our algorithm are not found purely real by the JT algorithm. Consequently, this comparison with the JT method does not provide a sufficiently convincing test for large values of *p*.

223

224

225

226

227

228

229

230

231

232

233

To test our algorithm for large values of *p*, first, we changed the values of the parameters *N* and *M*. We found that when *N* is larger than $2^{23}$ and *M* is larger than $2^{20}$, the number of identified zeroes in the domain of concern becomes stable for $p \le 21$. However, this was not the case for larger values of *p*, for which the number of zeroes and some numerical values were not stable when other values of *N* and *M* were chosen. This suggests that for large values of *p*, the zeroes could be closer to each other than the available numerical precision. We therefore considered that our method is not able to give all solutions for $p \ge 22$, and we stopped our systematic numerical exploration of superstability at $p=21$. However, there is no problem in the determination of the first zeroes of the superstability polynomial, which are sufficiently separated from each other, and, while complete solutions cannot be guaranteed for values of $p \ge 22$, robust values are certainly determined for the first zeroes of the superstability polynomial for *p* as large as 2048.

234

235

236

237

238

239

For $p \le 21$, however, we get further confidence based on another non-trivial test. Indeed, all zeroes found for $P_p(\beta)$ must be found for $P_q(\beta)$ where *q* is an integer multiple of *p*. For example, all zeroes found in a given range for $P_3(\beta)$ must be found for $P_6(\beta)$, $P_9(\beta)$, $P_{12}(\beta)$, $P_{15}(\beta)$, $P_{18}(\beta)$ and $P_{21}(\beta)$, within the above suggested precision of $10^{-15}$. For $N=2^{23}$ and $M=2^{20}$, this test was satisfied, with a few exceptions for $p=7$ to $p=21$ where the difference between the expected identical zeroes could increase to $10^{-14}$. Given the large number of combinations studied, the agreement can be
240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

considered as remarkable, especially given the fact that the $P_{10}(\beta)$ and $P_{20}(\beta)$ polynomials, for instance, are extremely different, and the numerical search for their roots largely independent.

# 4. Results

## 4.1 Overview and presentation of results

The parameters and numerical results for the superstable orbits for $p \le 21$ are summarized in Table 1, while the corresponding complete set of superstable parameters are given as Supp. Mat. Table S1. Another a-posteriori evidence for the satisfactory results given by our algorithm comes from the fact that the number $X_p$ of superstable orbits at period $p$ was calculated exactly by Guckenheimer in 1979 [40]. When $p$ is a prime number, we have $X_p=(2^{p-1}-1)/p$, giving $X_{11}=93$, $X_{13}=315$, $X_{17}=3855$ and $X_{19}=13797$, which coincide exactly with the values $Z_{11}$, $Z_{13}$, $Z_{17}$, and $Z_{19}$ given by our numerical algorithm (Table 1). In contrast, the number of periodic windows found by Tucker and Wilczak [32] are smaller than the theoretical prediction for $p=13$, 17 and 19. When $p$ is not a prime number, a simple formula for $X_p$ is not available, but the value of $B_p$ can be determined iteratively by a simple algorithm [40]. The obtained values, which exactly match the numbers given by Galias [34], are given in Table 1 and differ only slightly from the obtained value of $Z_p$, suggesting that some superstable values differing by the numerical resolution are either not distinguished or are artificially duplicated.

The smallest superstable $\alpha$ value at period $p$ is noted $\gamma_p$, considering the first non-trivial value, namely the first value which is not already a superstable value for a smaller period $q$ divisor of $p$. The values of $\gamma_p$ are given in Table 1 for $p \le 21$. The first two values are given in Supp. Mat. Table S2 for $p \le 1000$, with the addition of $p=1024$ and $p=2048$. Using these results, we can calculate the ratio $(\gamma_{1024}-\gamma_{12})/(\gamma_{2048}-\gamma_{1024}) \approx 4.6692012$, which is in good agreement with the known value 4.6692016091029 of the Feigenbaum constant $\delta[15,21]$. This fact gives some confidence that our
266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

numerical results for the first two superstable parameters (Table S2) are indeed reliable till $p=1000$ at least. Additional calculations were performed in selected $\alpha$ intervals, when studying particular properties (see below), then up to period $p=32$.

## 4.2 First superstable orbits and fine scaling structure

The $\gamma_p$ values are remarkably organized on the real axis (Fig. 3). The $\gamma_p$ values for $p=2^k$ appear in increasing order, following the well-known period-doubling scheme [8,15], with an approximative scaling law $\gamma_p \approx 3.57083-1.785/p^{2.30}$. The $\gamma_p$ values for even values of $p$ which are not powers of 2 then appear in two several bands, in decreasing order, as predicted by Sharkovskii ordering [17]. The $\gamma_p$ values for odd values of $p$, including all prime numbers except 2 (Fig. 3), also appear in decreasing order, in a single group, with an approximative scaling law $\gamma_p \approx 3.67808+1.8416/p^{2.28}$. Period 3 appears last, according to the well-known fact mentioned above [8,16,28]. An interesting domain emerges in Fig. 3 for $3.57083 \le \alpha \le 3.6786$, a domain where no odd period cycle seems to appear.

When the period $p$ increases, the number of superstable points (or the number of stable periodic windows) increases fast like $2^{p-1}/p$. Nevertheless, a detailed organisation of the stable points is preserved, with unexpected scaling laws, for example, when the difference $\Delta\gamma_p$ between the first and the second superstable parameter ($p \ge 4$) is considered (Fig. 4). When the period increases, $\Delta\gamma_p$ decreases fast for odd periods, reaching the computational resolution for $p$ smaller than 77. In contrast, even periods populate the domain between this minimum value and a maximum boundary value, obtained when the period is a power of 2, which decreases slowly with increasing period. This maximum difference follows a scaling law of the form $p^{-2.245}$, for which, nevertheless, no simple explanation could be found. The particular role of the periods of power 2 displays here another surprising form.
Fine scaling structures are found in other aspects of the superstable orbit parameters. For example, it is interesting to consider the interval $3.5972382<\alpha<3.6785736$ where no odd periodic orbit can be found (Fig. 5). In this interval, the number of superstable orbits with odd period increases exponentially, approximately following $10^{0.13p}$. In the interval $3.6785736<\alpha<3.69$, in contrast, both even and odd period superstable orbits are found. Similar scaling laws for the number of superstable orbits versus period are observed, but with slightly different slopes for even and odd periods (Fig. 5). In both intervals, the observed scaling relations for the smaller values of periods allows to make predictions for any value of $p$. Thus, while a complex fine structure can be present when focusing on particular properties, it is basically established for low periods, and is then scaled upwards without novel features.

301

302

## 4.3 Global scaling properties

303

304

305

306

307

308

309

310

311

Global scaling properties are observed when considering the overall set of superstable orbits. Indeed, while the distributions of the superstability points (zeroes of $P_p$), shown in Fig. 6 for $15 \le p \le 21$, bear an irregular structure, with some deficit zones, in particular for the largest values of $\beta$, for small values of $\beta$ ($10^{-9} \le \beta \le 10^{-3}$), the distributions display a scale-free structure, with a slope varying from 0.4866 for $p=15$ to 0.4964 for $p=21$, thus close to 1/2. The average value $B_p$ and the standard deviation $S_p$ of the distribution vary slowly with $p$ (Table 1). The average value $B_p$, which is here numerically similar to the median, varies from 0.0367 for $p=15$ to 0.0218 for $p=21$, and the following trend is found: $B_p = 0.0753252-0.002586p$.

312

313

These distributions of the bulk of superstability points are well represented by a gamma distribution $G_k(\beta)$ (Fig. 6):

314

$$G_k(\beta) = \frac{1}{\mu\Gamma(k+1)} \left(\frac{\beta}{\mu}\right)^k e^{-\frac{\beta}{\mu}}, \qquad (8)$$
where $k=0.5$, $\Gamma$ is the gamma function and $\mu$ the scale parameter. The parameter $\mu$ is estimated from the standard deviation $S_p$ from $\mu = S_p/\sqrt{1+k}$.

The largest interval between consecutive superstable orbits, however, does not display any scaling behaviour (Fig. 7). This quantity, indeed, is largely constrained by the largest period 3 to 5 windows, and rather it shows a semi-stable structure, with values of the order of 0.2 for even numbers and of 0.03 for odd numbers. The smallest interval, by contrast, it not dominated by the effect of prior periodic intervals; it decreases exponentially with increasing period, with slopes slightly different for even and odd periods (Supp. Mat. Fig. S2).

323

324

## 4.4 Scaling of the ultimate superstable orbits

325

On the left-hand side of the distribution (Fig. 6), the smallest values of $\beta$ for the identified superstability points seem to lie outside the scale-free tail of the gamma distribution. Among these outliers, the smallest value, noted $\epsilon_p$ (Table 1), shows a remarkable scaling behaviour versus periodicity $p$ (Fig. 8a), with a slope converging to -4. This scaling behaviour of $\epsilon_p$ for sufficiently large values of $p$, was predicted by Lorenz [7]:

331

$$ \epsilon_p = \frac{\pi^2}{4p} \qquad (9) $$

Eq. (9) gives $0.359054 \times 10^{-10}$, $0.089764 \times 10^{-10}$ and $0.022441 \times 10^{-10}$, for $p=19, 20$ and $21$, respectively. These values are almost identical to the numerical results of Table 1; this fact, incidentally, provides an additional check of the accuracy of our numerical method.

The proof given by Lorenz deserves to be reproduced here [7, page 8]. It relies on the fact that, for $\alpha=4$, the variable $x_n$ can be written $x_n=\sin^2(\pi y_n)$, with the $x_n$ recurrence transferring to dyadic transformation of $y_n$. Taking the set of superstable periodic points $x_i$, with $x_0=1/2$, we have, on the one hand, $x_1=\alpha/4=1-\epsilon_p/4$ and $x_2=\alpha^2/4(1-\alpha/4)=\epsilon_p$ keeping only first order terms in $\epsilon_p$. On the other hand, since $x_p=x_0=1/2=\sin^2(\pi/4)$, we can work out backwards the dyadic equivalence, which leads
to $x_{p-1}=\sin^2(\pi/8)$, $x_{p-2}=\sin^2(\pi/16)$ and so on, until $x_2=\sin^2(\pi/4/2^{p-2})=\sin^2(\pi/2^p)$. Identifying the two expressions of $x_2$ gives the limit of Eq. (9), rapidly converging when $p$ is increasing.

The theoretical prediction Eq. (9) is compared with the numerical calculation in Fig. 8a. In Fig 8b, we show the difference $\Delta_p$, versus $p$, of the numerical estimate minus the Lorenz prediction of Eq. (9). The numerical values found for $\epsilon_p$ converge fast to the theoretical Lorenz scaling expression. Without any prior knowledge, the scaling factor $4^p$ of $\epsilon_p$ and even an accurate expression for $\epsilon_p$ can be derived using an elementary method (Supp. Mat. Appendix S1).

Not only the largest superstable parameter can be calculated. Indeed, expanding the original result of Lorenz, Simó and Tatjer proved in 1991 that, for any quadratic function, when the period $p$ is large enough, all the largest superstable values near the maximum parameter value can be calculated exactly from a functional equation. In the case of the logistic map, when the maximum values are written $4-\beta_{p,j}$, from their Corollary 1 (page 339 of [22]), the following expression is obtained:

$$ \beta_{p,j} = \frac{\pi^2}{4p}(1+2j)^2 + o(\frac{1}{4^2p}), \qquad (10) $$
ignoring here the higher order terms scaling like $4^{-2p}$, whose explicit expression is also known [22].

To compare our numerical results with the Simó and Tatjer scaling Eq. (10), we calculate the expression $\Omega_p(j)$ defined by:

$$ \Omega_p(j) = \left( \frac{1}{1+2j} \sqrt{\frac{\beta_{p,j}}{\epsilon_p}} - 1 \right) \times 100, \qquad (11) $$

which quantifies, in %, the difference between the numerical result and the multiplier $(1+2j)^2$ predicted by Eq. (10). The obtained values of $\Omega_p(j)$ (Fig. 8c) for $j=1$ to 4 are all smaller than $10^{-2}$ % for period $p$ larger than 10, and, in general, Eq. (10) is compatible with the numerical results to the limit of the numerical precision. For $p=21$, for example, we have obtained $\beta_{21,1}=2.01967 \times 10^{-11}$, $\beta_{21,2}=5.61013 \times 10^{-11}$, $\beta_{21,3}=1.099605 \times 10^{-10}$ and $\beta_{21,4}=1.817719 \times 10^{-10}$ while Eq. (10) gives $2.01968 \times 10^{-11}$, $5.61022 \times 10^{-11}$, $1.099603 \times 10^{-10}$ and $1.817711 \times 10^{-10}$, respectively. In fact, a more
364 detailed comparison (Supp. Mat. Fig. S3-S5) shows that the ultimate numerical values of the
365 superstable parameters are in remarkable agreement with the expanded Simó and Tatjer scaling
366 prediction Eq. (10) [22] for larger values of j, even for moderately large values of the period p (p≈10
367 to 12), and for increasingly larger values of j for large periods (p≥16). For example,
368 β<sub>21,100</sub>=9.06634354×10<sup>-8</sup> while Eq. (10) gives 9.06634085×10<sup>-8</sup>.

369

370

4.5 Superstable periodic time-series

371

372 Once the superstable values are found, it is possible to display the periodic sets xᵢ, with
373 x₀=1/2 (Fig. 9), and systematic patterns are observed. For the largest superstable α=4-εₚ, the time-
374 series always show a slowly increasing phase, with a rapid growth to almost 1 at the end of the
375 cycle (red curves). This is the typical pattern of the unstable dyadic periodic solutions for α=4,
376 reminiscent of the temporal structure of a relaxation oscillator [41], with a slow growth phase
377 followed by a rapid return to the base level.

378 By contrast, the time-series for the first superstable periodic cycle (α = γₚ) is an oscillation
379 with a growing amplitude around the unstable fixed point x₁=1-1/α (blue curves in Fig. 9), and the
380 oscillation stops when the negative amplitude exceeds the available positive range. The growth
381 phase of the amplitude of the oscillation is exponential, with a slow initiation which is longer for
382 larger periods (case p=21 shown in Fig. 9e).

383 For superstable cycles with α near the maximum of the distribution (Fig. 6), the temporal
384 structure is intermediate, with internal subcycles, as shown in the case selected in Fig. 9c for p=13.
385 These various types of time-series may deserve to be looked into in more detail when confronted
386 with particular experimental time-series in mechanical or biological systems.

387 For large values of p, some time-series do not show stable periodicity, even for the
388 superstable value of the parameter. This behaviour is first observed for the two ultimate superstable
389 parameters for p=15. For example, the time-series obtained for ε₁₅ becomes unstable after
completing only five cycles (Fig. 10a). The amplitude spectrum corresponding to this time-series (Fig. 10b) barely shows the fundamental $1/p$ frequency and its harmonic peaks, almost invisible in blank noise. This is due to the presence of cumulated numerical instabilities, illustrating the fact that the numerical orbits do not necessarily coincide with the true orbits [42]. In contrast, the time-series corresponding to the third ultimate superstable orbit $\beta_{15,3}$ is stable, with an amplitude spectrum displaying only the fundamental and harmonic peaks of a rigorously periodic function (Supp. Mat. Fig. S6). When larger values of the period are considered, the instability of the time-series corresponding to superstable parameters is observed more and more frequently, and the instability appears more rapidly (after completing less cycles) for larger values of $\alpha$ (Supp. Mat. Fig. S7).

The instability of the superstable time-series for large periods illustrates the limitations due to digital precision [38,42], and raises the need for caution when drawing conclusions on the periodic or chaotic nature of particular states on the basis of time-series only. However, the fact that the numerical logistic map cannot be identified with the theoretical ideal logistic map can also be considered as an asset when trying to evaluate the capacity of chaotic models to approach reality. Indeed, while the theoretical ideal properties are interesting to study in details, for example to give explanations for observed scaling properties [7,22], only the properties maintained in the numerical time-series may be relevant to approach experimental situations or natural processes. The question of numerical stability also raises the question of the stability domain of the superstable orbits.

409

410

## 4.6 Stability domain of the superstable orbits

411

412

The condition of superstability only guarantees that the time-series is exactly periodic when the initial point $x_0$ is $1/2$, or, more generally, when one of the iterated $x_n$ is $1/2$ (superstable time-series). When starting from other initial values, the time-series may converge or not to the stationary periodic time-series, depending on the structure of the attraction basin of the considered cycle

413

414

415
[20,21]. While the detailed topological structure of the attraction basin is outside the scope of the present study, the behaviour near point can be established empirically relatively easily.

Writing this initial point $x_0=1/2+u$, after $p$ iterations, we arrive at the point $x_p=f^{(p)}(1/2+u)$. Periodicity of period $p$ and the superstability condition impose $f^{(p)}(1/2)=1/2$ and the fact that the first derivative of $f^{(p)}(x)$ at point $1/2$ is zero. Thus, the function $f^{(p)}(1/2+u)$ of $u$ must increase (in absolute value) from 0, with a zero slope at the origin. Consequently, it must cross the line $1/2+u$ at some value $u_s$, giving:

$$f^{(p)}\left(\frac{1}{2} + u_s\right) = \frac{1}{2} + u_s. \qquad (12)$$

Examples of $f^{(p)}(1/2+u)$ are displayed in Supp. Mat. Fig. S8. For $u$ smaller than $u_s$, the distance to $1/2$ decreases after $p$ iterations, thus the time-series converges to the superstable time-series, while the distance to $1/2$ increases when $u$ is larger than $u_s$. Thus, $u_s$ defines the half-width of the superstable orbit. When $p$ is increased, and, for a given $p$, when increasing values of the superstable $\alpha$ are selected, the function shows a faster increase and the intersection is reached for smaller values of $u_s$.

The value of the half-width $u_s$ can be first estimated numerically for each superstable value of $\alpha$ (Fig. 11). For large values of $p$ and small values of $\beta$, $u_s$ becomes smaller than $10^{-9}$, and numerical instabilities appear. Thus, a better estimate of the intersection is given by extrapolating the function $f^{(p)}(1/2+u)-1/2-u$ to zero from values of $u$ larger than $10^{-8}$. The results (Fig. 11) are remarkably organized as a function of $p$ and $\beta$, with conspicuous power laws for the smaller values of $\beta$, in agreement with theoretical predictions [22].

The half-width $u_s$ can also be approximated in the following simple manner. Consider the iterates of initial point $x_0=1/2+u$. First iterate is $f(1/2+u)=\alpha(1/4-u^2)$. Then, consider the next iterates with only terms up to $u^2$, thus writing:

$$f^{(k)}\left(\frac{1}{2} + u\right) = x_k^p + \alpha^k b_k u^2, \qquad (13)$$
where the $x_i^p$ are the periodic points of the considered superstable orbit. From the first iterate above, we have $b_1=-1$. The next iteration gives:

$$f^{(k+1)}(\frac{1}{2} + u) = x_{k+1}^p + \alpha^{k+1} b_{k+1} u^2 = \alpha(x_k^p + \alpha^k b_k u^2)(1 - x_k^p - \alpha^k b_k u^2). \quad (14)$$

We have $x_{k+1}^p = \alpha x_k^p (1 - x_k^p)$ and, identifying the $u^2$ terms, we get:

$$b_{k+1} = b_k (1 - 2x_k^p), \quad (15)$$

hence, applying recurrence Eq. (15) and initial value $b_1=-1$, we have:

$$b_k = - \prod_{i=1}^{k-1} (1 - 2x_i^p). \quad (16)$$

For $k=p$, we have $x_k^p = 1/2$, thus we obtain the half-width as:

$$u_s = \frac{\alpha^{-p}}{|b_p|} = \frac{\alpha^{-p}}{|\prod_{k=1}^{p-1}(1 - 2x_k^p)|}. \quad (17)$$

The predictions from Eq. (17) agree remarkably well (Fig. 10, full lines) with the previous numerical estimates. For small values of $\beta$, the theoretical prediction of Eq. (17) displays a simple scaling behaviour:

$$u_p(\beta) = \frac{\sqrt{\beta}}{2^{p+1}}. \quad (18)$$

For the ultimate superstable orbit ($\beta=\epsilon_p$), using Lorenz scaling Eq. (9), the half-width becomes:

$$u_p = \frac{\pi}{2^{2p+1}}. \quad (19)$$

This prediction agrees well with the numerical estimates and with the more general scaling given by Simó and Tatjer [22]. For large periods, this half-width becomes extremely small ($3.6 \times 10^{-13}$ for $p=21$), comparable or smaller than the numerical accuracy. Thus, for large $p$ and small $\beta$, the numerical time-series are not expected to represent reliably the true superstable orbits [42]. Nevertheless, the existence of these orbits is established and their scaling properties are known, without relying on the time-series.

## 4.7 Width of the periodic window
464

465

466

467

468

469

The half-width $u_s$ of the superstable orbit is closely related to the width $\Delta\alpha_p$ of the periodic orbit which is the domain from $\alpha_p^{\min}$ to $\alpha_p^{\max}$, around the superstable orbit of parameter $\alpha_p$, where the periodic series is stable. Indeed, for $u$ small, the series $y_k^p$ initiated with $y_0^p = \frac{1}{2} + u$ is periodic for a given $\alpha=\alpha_p+\delta$, with $\delta$ small and adequately chosen, and we have:

$$y_p^p = y_0^p = \frac{1}{2} + u = f^{(p)}\left(\frac{1}{2} + u\right). \qquad (20)$$

470

471

We have $y_1^p = \alpha(\frac{1}{4} - u^2)$, thus all $y_k^p$ with $k \ge 1$ differ from $Q_p(\alpha)$ by terms of order larger than $u^2$, hence, keeping only the terms in first order in $u$:

472

$$\frac{1}{2} + u = f^{(p)}\left(\frac{1}{2} + u\right) = Q_p(\alpha) = Q_p(\alpha_p) + \left.\frac{dQ_p}{d\alpha}\right|_{\alpha_p} \delta. \qquad (21)$$

473

hence:

474

$$\delta(u) = \left.\frac{u}{\frac{dQ_p}{d\alpha}}\right|_{\alpha_p}. \qquad (22)$$

475

For this periodic solution to be stable, however, $u$ has to be smaller than $u_0$, hence giving:

476

$$|\alpha^p \prod_{j=0}^{p-1} (1 - 2y_j^p)| = 1. \qquad (23)$$

477

In first order in $u_0$, we have:

478

$$\prod_{j=0}^{p-1}(1 - 2y_j^p) = (1 - 2y_0^p) \prod_{j=1}^{p-1}(1 - 2y_j^p) \cong 2u_0 \prod_{j=1}^{p-1}(1 - 2x_j^p) = 2 \frac{u_0}{u_s} \alpha^{-p}, \qquad (24)$$

479

480

where we have recognized the partial product from the expression of $u_s$ of Eq. (17). Incorporating into Eq. (22) and Eq. (23), $\Delta\alpha_p=2\delta(u_0)$ has the following expression:

481

$$\Delta\alpha_p = \left.\frac{u_s}{\left.\frac{dQ_p}{d\alpha}\right|_{\alpha_p}}\right|_{\alpha_p}. \qquad (25)$$

482

483

484

For $p=3$, for example, we have $u_s \cong 2.84 \times 10^{-2}$ and $\left.\frac{dQ_p}{d\alpha}\right|_{\alpha_p} \cong 2.08$, giving $\Delta\alpha_p \cong 0.0136$, while the exact calculation gives 3.841499–3.828427 $\cong 0.0131$ (see section 1 above and [28]). The width of the periodic windows given by Eq. (25) is scaling faster than $u_s$. For example, for the ultimate
485

periodic window of period 13 ($\beta_p \approx 1.47 \times 10^{-7}$; Table 1), we get $\Delta\alpha_p \approx 6.3 \times 10^{-15}$. The periodic windows get narrower than the numerical accuracy [32,34,35].

486

487

488

## 5. Conclusion

489

490

In this paper, we have investigated numerically the values of the parameter $\alpha$ of the logistic map leading to superstable periodic orbits, and determined approximate values, reliable to 14 digits, comprehensively for period $p \le 21$ (Supp. Mat. Table S1), and for the first two superstable orbits for $p \le 1000$, with the addition of $p=1024$ and $p=2048$ (Supp. Mat. Table S2). Our work complements the results already obtained recently by other workers [32-36] by focussing on the role of scaling in the fine and global structure of superstable orbits.

491

492

493

494

495

496

The superstability condition is equivalent to the search of the roots of a degree $D_p=2^p-1$ superstability polynomial $P_p(\beta)$, a task performed numerically. While the number of superstable points $Z_p$ increases fast with $p$, like $2^{(p-1)}/p$, and is exactly $(2^{(p-1)}-1)/p$ for prime numbers [34-35,40], our numerical algorithm provides nearly the expected number of zeroes for most values of $p$, and exactly for all prime numbers until $p=19$ included.

497

498

499

500

501

The superstable parameter values display various fine structures with scaling properties, sometimes slightly different for $p$ even and odd, and, aside from odd-even and prime numbers, some values of the period appear to play a particular role. For example, when considering the difference between the first and the second superstable parameters, an upper bound scaling with period is observed when $p$ is a power of 2, a property probably linked to the structuring role of the bifurcations by period doubling. When a fine structure is established for moderate values of period $p$, predictions to larger values of $p$ can be attempted using scaling relations. This also suggests that expanding the determination of superstable orbits to larger values of the period, a task which is increasingly difficult because of the challenging numerical precision issues, may not lead to any novel features, the basic structures being already revealed from periods around $p=15$ to 17.

502

503

504

505

506

507

508

509

510
511

512

513

514

515

516

517

518

519

520

Intertwined with the fine structures, global scaling properties are also observed. Indeed, the
superstability points are distributed following a gamma distribution in terms of $\beta=4-\alpha$, with
parameter 0.5. The last superstability orbit, characterized by $\beta=\epsilon_p$, is in remarkable agreement with
the scaling equation $\pi^2/4p$ predicted by Lorenz in 1964 [7]. The following large superstability
parameters $4-\beta_{p,j}$ also follow a more general scaling relation established by Simó and Tatjer in 1991
[22], which, to a first approximation, simply adds a $(1+2j)^2$ multiplier to the Lorenz scaling, hence
providing an additional verification of the accuracy of the numerical determination of superstable
parameters. Other parameters, such as the width of the superstable orbits or the related width of the
periodic window also display remarkable scaling properties for large values of $p$, observed
numerically and understood theoretically [22].

521

522

523

524

525

526

527

528

529

After the superstable parameters are determined, the periodicity of the time-series can be
checked, and, while the time-series are precisely periodic, as expected, for all the superstable orbits
till $p=14$, instabilities appear for the two larger superstable parameters found for $p=15$ and for an
increasingly larger fraction of parameters for $p \ge 16$. This fact is due to the limitations in the
representation of real numbers in computers, even in double precision. This inherent feature of the
numerical logistic map, which differs essentially from the ideal theoretical function, may not only
be a limitation, but also an interesting feature to concentrate on practical properties preserved by
the numerical model, properties which are potentially more able to bear some relevance in
experimental situations and natural systems.

530

531

532

533

534

535

536

Numerous properties of these numerical superstable orbits remain poorly known, in particular
their attraction basins. Theoretical topological properties of the attraction sets of the periodic orbits
have been established in the logistic map since decades [20-22], but the systematic study in the case
of explicit superstable orbits is not an easy task. First results, using the numerical superstable
parameters (see Supp. Mat. Figs. S9 and S10), indicate that, for low values of $p$, or, for larger values
of $p$, for the smallest values of $\alpha$, the time-series always converge to the periodic time-series for all
choices of the initial point $x_0$. Large dispersion of the relaxation time, however, can be observed. In
other cases, in particular for large values of *p* and small values of *β*, no relaxation to the periodic orbit seems possible, even after large number of iterations. This seems to be the general case for the ultimate orbit with $β=ε_p$. More work will be needed to establish the properties of these numerical attraction basins, their possible fractal structure, whether they confirm the theoretical topological properties [19-20], and whether the assessment of these properties is affected by digital precision.

The dominant presence of superstable orbits, around the peak of the distribution of superstable points, or the total absence of superstable orbits of period *p* or smaller for $β<β_p$, must have consequences on the nature of chaotic regimes that can compete with periodic orbits. While some consequences can be reasonably expected, for example in the case of transient intermittent bursts in the time-series, they remain speculative at this stage and the practical consequences of the established theoretical framework on the interaction between SIC properties, intermittent bursts and superstable periodicity remains to be worked out [9]. Similarly, the relevance and actual meaning of scaling properties must be evaluated. These issues are not only relevant to legitimately refine the mathematical understanding of periodicity in chaotic models, including quadratic recurrence. Indeed, windows of periodic regimes, stable or transient, have been observed in natural systems [43-44], and might be an unavoidable ingredient of natural non-linear systems [9-10,14], intimately linked to intermittency regimes [45,46] to which an increasing attention is being dedicated [i.e., 47]. Superstable periodic orbits, which can be studied analytically, may reflect some important features of chaotic regimes that might be more difficult to handle theoretically. Indeed, the interplay between periodic and chaotic states deserves more attention to understand the temporal evolution of complex systems, which will be essential to progress in the understanding of the predictability of natural systems affecting our life and society [10-14].

**Declaration of Competing Interest**

The authors declare no competing interests.
563

564

CRediT authorship contribution statement

565

566

**Frédéric Perrier:** Initiation, numerical solutions, software, analytical results and figures, writing. **Frédéric Girault:** Critical analysis, verification, validation, writing.

567

568

569

Acknowledgements

570

571

572

573

574

575

576

The authors thank Pierrick Chanteux, Johann Boudot and Josef Schomburg for their availability in all circumstances and their expertise in computer technology and software installation. The original version of this manuscript was substantially modified thanks to the detailed work of two anonymous reviewers and the communication of important papers by Joan Carles Tatjer and Zbigniew Galias. This work contributes to the IdEx Université de Paris ANR-18- IDEX-0001.

577

578

References

579

580

[1] Verhulst P-F. Notice sur la loi que la population suit dans son accroissement. Corresp Math Phys 1838;10:113-126.

581

582

[2] Pelinovky E, Kurkin A, Kurkina O, Kokoulina M, Epifanova A. Logistic equation and COVID-19. Chaos Soliton Fract 2020;140:110241.

583

584

585

[3] Wu K, Darcet D, Wang Q, Sornette D. Generalized logistic growth modeling of the COVID-19 outbreak: comparing the dynamics in the 29 provinces in China and in the rest of the world. Nonlin Dyn 2020;101:1561-1581.

586

587

[4] Schröder E. Über iterirte Functionen. Math Annal 1870;3:296-322.

588

589

590

[5] Julia G. Mémoire sur l'itération des fonctions rationnelles. J Math Pures App 1918;8:47-246.

591

592

[6] Fatou P. Sur les équations fonctionnelles. Bull. Soc Math France 1919;47:161-271.

593

594

[7] Lorenz E.N. The problem of deducing the climate from the governing equations. Tellus 1964: 16;1-11.

593

594

[8] May RM. Simple mathematical models with very complicated dynamics. Nature 1976;261:459-467.

595

596

597

[9] Strogatz SH. Nonlinear dynamics and chaos: With applications to physics, biology, chemistry, and engineering. Westview Press Boulder Colorado. Second Edition; 2015.

[10] Ghil M. A century of nonlinearity in the geosciences. Earth Space Sci 2019;6:1007-1042.

[11] Allègre CJ, Le Mouël J-L, Provost A. Scaling rules in rock fracture and possible implications for earthquake prediction. Nature 1982;297:47-49.
598

599

[12] Faybishenko B. Nonlinear dynamics in flow through unsaturated fractured porous media: Status and perspectives. *Rev Geophys* 2004;42:RG2003.

600

601

[13] Phillips JD. Deterministic chaos and historical geomorphology: A review and look forward. *Geomorphology* 2006;76:109-121.

602

603

[14] Turcotte, DL. Fractals and chaos in geology and geophysics. Cambridge University Press, New York ; 1992.

604

605

[15] Feigenbaum MJ. Quantitative universality for a class of nonlinear transformations. *J Stat Phys* 1978;19:25-52.

606

[16] Li T-Y, Yorke JA. Period three implies chaos. *Am Math Monthly* 1975;82: 985-992.

607

608

[17] Sharkovskii AN. Coexistence of cycles of a continuous map on the line into itself. *Ukrain. Math. Zh.* 1964;16:61-71.

609

610

[18] Metropolis N, Stein ML, Stein PR. On finite limit sets for transformations on the unit interval. *J Combin Theor* 1973;15:25-44.

611

612

[19] Tsuchiya T, Yamagishi D. The complete bifurcation diagram for the logistic map. *Zeitschr Naturforsch* 1997;52a:513-516.

613

614

[20] Douady A, Hubbard J. Etude dynamique des polynômes complexes. *Documents mathématiques, Société Mathématique de France* ; 2007.

615

616

[21] de Melo W., van Strien S. One-dimensioanl dynamics. Ergebnisse der Mathematik and ihrer Grenzgebiete (3) V25 ; 1993.

617

618

[22] Simó C., Tatjer JC. Windows of attraction of the logistic map. *European Conference on Iteration Theory* (Batschuns, 1989) 1991:335-342.

619

620

[23] Lyubich M. Almost every real quadratic map is either regular or stochastic. *Ann of Math* 2002;156:1-78.

621

[24] Saha P, Strogatz SH. The birth of period three. *Math Mag* 1995;68:42-47.

622

[25] Gordon WB. Period three trajectories of the logistic map. *Math Mag* 1996;69:118-120.

623

624

[26] Lee H. Analytical study of the superstable 3-cycle in the logistic map. *J Math Phys* 2009;50:122702.

625

[27] Zhang C. Period three begins. *Math Mag* 2010;83:295-297.

626

627

[28] Fidélis AJ, Martins LC. Algebraic orbits on period-3 window for the logistic map. *Nonlin Dyn* 2015;79:1015-1021.

628

629

[29] Hoppenstreadt F, Hyman JM. Periodic solutions of a logistic difference equation. *Siam J Appl Math* 1977;32:73-81.

630

631

[30] Brown A. Equations for periodic solutions of a logistic difference equation. *J Austral Math Soc* 1981;23:78-94.

632

633

[31] Brown A. Solutions of period seven for a logistic difference equation. *Bull Austral Math Soc* 1982;26:263-284.

634

635

[32] Tucker W, Wilczak D. A rigorous lower bound for the stability regions of the quadratic map, *Physica D* 2009;238:1923-1936.

636

637

[33] Golmakani A, Koudjinan CE, Luzzatto S, Pilarczyk P. Rigorous numerics for critical orbits in the quadratic family. *Chaos* 2020;30:073143.

638

639

[34] Galias Z. Rigorous numerical study of low-period windows for the quadratic map. *Int J Bifurc Chaos* 2015;25:1550139.

640

641

[35] Galias Z. Systematic search for wide periodic windows and bounds for the set of regular parameters for the quadratic map. *Chaos* 2017;27:053106.

642

643

[36] Provatas N, Mackey MC. Asymptotic periodicity and banded chaos. *Physica D* 1991;53:295-318.

644

645

646

[37] Diaz-Ruelas A, Baldovin F, Robledo, A. Logistic map trajectory distributions: Renormalization-group, entropy and criticality at the transition to chaos. *Chaos* 2021;31:033112.

647

648

[38] Galias Z. Periodic orbits of the logistic map in single and double precision implementations. *IEEE Trans Circ Syst II* 2021;68:3471-3475.
649

650

[39] Jenkins MA, Traub JF. A three-stage algorithm for real polynomials using quadratic iteration. SIAM J Num Anal 1970;7:545-566.

651

652

[40] Guckenheimer J. The bifurcation of quadratic functions. Ann New York Acad Sci 1979;316:75-85.

653

[41] van der Pol B. On relaxation oscillations. London Edinb Dublin Phil Mag 1926;2:978-992.

654

655

[42] Hammel SM, Yorke JA, Grebogi C. Do numerical orbits of chaotic dynamical processes represent true orbits? J Complex 1987;3:136-145.

656

657

[43] Morat P, Le Mouël J-L. Electrical signals generated by stress variations in porous non saturated rocks. C R Acad Sci Paris 1992;315:955-963.

658

659

[44] Perrier F, Le Mouël J-L. Stationary and transient thermal states of barometric pumping in the access pit of an underground quarry. Sci Total Environ 2016;550:1044-1056.

660

[45] Manneville P., Intermittency in dissipative dynamical systems. Phys Lett A 1980;79:33-35.

661

662

[46] Constantin P, Procaccia I. Dimension of the carrier of turbulence in fluid mechanics. Phys Rev A 1992;46:4736-4741.

663

664

665

[47] Diaz-Rueslas A, Jensen HJ, Piovani D, Robledo A. Tangent map intermittency as an approximate analysis of intermittency in a high dimensional fully stochastic dynamic system: The Tangled Nature model. Chaos 2016;26:123105.

666

667
<table><caption>**Table 1.** Overview of results obtained for the superstable periodic orbits with period $p \le 21$. $D_p$ is the degree of the superstability polynomial, $X_p$ is the expected exact number of superstable values [35-36,40], and $Z_p$ is the number of zeroes found for $\alpha$ between 3.1 and 4. The first non-trivial superstable orbit is found for $\alpha=\gamma_p$ (smallest superstable value which is not a superstable value for a period $q$ divisor of $p$). The ultimate (largest) superstable orbit is found for $\alpha=4-\epsilon_p$. For $p \ge 15$, $B_p$ and $S_p$ are the mean value and the standard deviation, respectively, of the superstable $\beta$ distribution (Fig. 6).</caption><thead><tr><th>p</th><th>D<sub>p</sub></th><th>X<sub>p</sub></th><th>Z<sub>p</sub></th><th>γ</th><th>ε<sub>p</sub></th><th>B<sub>p</sub></th><th>S<sub>p</sub></th></tr></thead><tbody><tr><td>2</td><td>3</td><td>1</td><td>1</td><td>3.2360679774997898</td><td>0.7639320225002102</td><td></td><td></td></tr><tr><td>3</td><td>7</td><td>1</td><td>1</td><td>3.8318740552833157</td><td>0.1681259447166843</td><td></td><td></td></tr><tr><td>4</td><td>15</td><td>3</td><td>3</td><td>3.4985616993277020</td><td>0.0397298727788473</td><td></td><td></td></tr><tr><td>5</td><td>31</td><td>3</td><td>3</td><td>3.7389149129706847</td><td>0.0097329530262988</td><td></td><td></td></tr><tr><td>6</td><td>63</td><td>7</td><td>7</td><td>3.6275575295155229</td><td>0.0024168817454324</td><td></td><td></td></tr><tr><td>7</td><td>127</td><td>9</td><td>9</td><td>3.7017691535379562</td><td>0.0006029390379023</td><td></td><td></td></tr><tr><td>8</td><td>255</td><td>19</td><td>19</td><td>3.5546408627688257</td><td>0.0001506379861498</td><td></td><td></td></tr><tr><td>9</td><td>511</td><td>29</td><td>29</td><td>3.6872161809341524</td><td>0.0000376523925003</td><td></td><td></td></tr><tr><td>10</td><td>1023</td><td>55</td><td>57</td><td>3.6470488015219011</td><td>0.0000094125878389</td><td></td><td></td></tr><tr><td>11</td><td>2047</td><td>93</td><td>93</td><td>3.6817186741371324</td><td>0.0000023531109172</td><td></td><td></td></tr><tr><td>12</td><td>4095</td><td>179</td><td>181</td><td>3.5822298358203599</td><td>0.0000005882752174</td><td></td><td></td></tr><tr><td>13</td><td>8191</td><td>315</td><td>315</td><td>3.6797028056802539</td><td>0.0000001470686315</td><td></td><td></td></tr><tr><td>14</td><td>16383</td><td>595</td><td>598</td><td>3.5972381983725583</td><td>0.0000000367671458</td><td></td><td></td></tr><tr><td>15</td><td>32767</td><td>1095</td><td>1095</td><td>3.6789763419033976</td><td>0.0000000091917856</td><td>0.0367 ± 0.0017</td><td>0.0555 ± 0.0012</td></tr><tr><td>16</td><td>65353</td><td>2067</td><td>2071</td><td>3.5666673798562676</td><td>0.0000000022979467</td><td>0.0352 ± 0.0013</td><td>0.06023 ± 0.00094</td></tr><tr><td>17</td><td>131071</td><td>3855</td><td>3855</td><td>3.6787167827358758</td><td>0.0000000005744867</td><td>0.03020 ± 0.00076</td><td>0.04724 ± 0.00054</td></tr><tr><td>18</td><td>262143</td><td>7315</td><td>7319</td><td>3.5942221098256288</td><td>0.0000000001436220</td><td>0.02823 ± 0.00054</td><td>0.04638 ± 0.00038</td></tr><tr><td>19</td><td>524287</td><td>13797</td><td>13797</td><td>3.6786244032684179</td><td>0.0000000000359046</td><td>0.02547 ± 0.00035</td><td>0.04084 ± 0.00025</td></tr><tr><td>20</td><td>1048575</td><td>26271</td><td>26257</td><td>3.5775498113692299</td><td>0.0000000000089764</td><td>0.02375 ± 0.00024</td><td>0.03928 ± 0.00024</td></tr><tr><td>21</td><td>2097151</td><td>49939</td><td>49941</td><td>3.6785915791011803</td><td>0.0000000000022440</td><td>0.02180 ± 0.00016</td><td>0.03564 ± 0.00011</td></tr></tbody></table>
677 **Figure captions**

678

679 **Fig. 1.** Stationary distribution of points for the logistic map $x_{n+1}=\alpha x_n(1-x_n)$ for $5000 \le n \le 10000$ with $x_0=0.5$, shown as a function of $\beta=4-\alpha$ with a logarithmic scale. The dash line represents the fixed point (period 1). The dotted lines represent the two points of the period 2 cycle. The dash-dot lines indicate the minimum and the maximum accessible values.

683

684 **Fig. 2.** Superstability polynomial of the logistic map $x_{n+1}=\alpha x_n(1-x_n)$ shown as a function of $\beta=4-\alpha$ with a logarithmic scale and various values of the period $p$. The vertical dotted lines correspond to zeroes of $P_p(\beta)$, in green, blue and red for $p=7, 8$ and $9$ respectively.

687

688 **Fig. 3.** First non-trivial superstable orbit $\gamma_p$ versus period $p$. It is the smallest superstable value which is not a superstable value for a period $q$ divisor of $p$. Periods which are a power of 2, multiple of 2 or odd numbers are shown with different colours. Prime numbers are shown in red. Lines refer to power-law fits (see text).

692

693 **Fig. 4.** Difference $\Delta\gamma_p$ between the first and the second non-trivial superstable orbit versus period $p$. Periods which are a power of 2, multiple of 2 or odd numbers are shown with different colours. Prime numbers are shown in red. The red line ($\Delta\gamma_p=11.92p^{-2.245}$) refers to a power-law fit of the points with power of 2 periods.

697

698 **Fig. 5.** Number of superstable points found in two intervals as a function of period $p$. Periods which are a power of 2, multiple of 2 or odd numbers are shown with different colours. Prime numbers are shown in red. In the first interval, only superstable orbits with even periods exist. The lines refer to power laws of the form $K \times 10^{bp}$ with $b=0.13, 0.15$ and $0.18$ for the orange, grey and blue lines, respectively.

699

700

701

702
703

704
705
**Fig. 6.** Distribution of the superstability points of the logistic map $x_{n+1}=\alpha x_n(1-x_n)$ shown as a function of $\beta=4-\alpha$ with a logarithmic scale and for values of the period $p$ between 15 and 21.

706

707
708
709
**Fig. 7.** Maximum interval $\Delta\alpha_{max}$ between consecutive superstable points as a function of period $p$. Periods which are a power of 2, multiple of 2 or odd numbers are shown with different colours. Prime numbers are shown in red.

710

711
712
713
714
715
716
**Fig. 8.** Numerical results for the superstability points of the logistic map $x_{n+1}=\alpha x_n(1-x_n)$ as a function of period $p$: a) Position $\epsilon_p$ of the last superstable points. The grey line represents Lorenz scaling Eq. (9) [7]; b) Difference (expressed in %) between the numerical result $\epsilon_p$ and Lorenz scaling Eq. (9); c) Normalized difference $\Omega_j$ given by Eq. (11) (expressed in %) between the numerical results and the Simó and Tatjer scaling Eq. (10) [22] for the smallest $\beta_{p,j}$ values for $j$ values from 1 to 4.

717
718
719

720
721
722
723
724
**Fig. 9.** Time-series for several examples of superstable orbits. The values of $\alpha$, taken from Table S1, are shown with the corresponding colour on the right-hand side.

725
726

720
721
722
723
724
**Fig. 10.** Numerical instability of periodic time-series for the ultimate superstable orbit for $p=15$. a) First 100 points of the time-series obtained from $x_0=1/2$ (black points) and last 100 points (blue points) for a total number of $2^{18}$ iterations; b) Stacked amplitude spectrum of the obtained time-series. Spectra for portions of 2048 points are stacked.

725
726

725
726
**Fig. 11.** Half-width of the superstable orbits of the logistic map $x_{n+1}=\alpha x_n(1-x_n)$ shown versus $\beta=4-\alpha$ for $p=7$ and for values of the period $p$ between 14 and 21. The diamonds refer to numerical estimates,
727 averaged in logarithmic spaced intervals versus $\beta$. The full lines refer to the theoretical prediction

728 Eq. (18). The horizontal dash lines refer to the predicted half-width of the ultimate orbit Eq. (19).

729