{'fraction_non_alphanumeric': 0.04452106698999697, 'fraction_numerical': 0.014436192785692635, 'mean_word_length': 4.331919191919192, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Universal scaling of Lyapunov-exponent fluctuations in space-time chaos

Finite-time Lyapunov exponents of generic chaotic dynamical systems fluctuate in time. These fluctuations are due to the different degree of stability across the accessible phase-space. A recent numerical study of spatially-extended systems has revealed that the diffusion coefficient D of the Lyapunov exponents (LEs) exhibits a non-trivial scaling behavior, D(L) ~ L^{-\gamma}, with the system size L. Here, we show that the wandering exponent \gamma can be expressed in terms of the roughening exponents associated with the corresponding"Lyapunov-surface". Our theoretical predictions are supported by the numerical analysis of several spatially-extended systems. In particular, we find that the wandering exponent of the first LE is universal: in view of the known relationship with the Kardar-Parisi-Zhang equation, \gamma can be expressed in terms of known critical exponents. Furthermore, our simulations reveal that the bulk of the spectrum exhibits a clearly different behavior and suggest that it belongs to a possibly unique universality class, which has, however, yet to be identified.


I. INTRODUCTION
Lyapunov exponents (LEs) λ i are the most powerful tool for a detailed characterization of chaotic dynamics [1]: they allow determining the number of unstable directions, the Kolmogorov-Sinai entropy, the fractal dimension and other dynamical invariants. The number of LEs equals the number of degrees of freedom N . The LEs, ordered from the largest to the smallest one, define the so-called spectrum of LEs. In spatially extended systems, N ∝ L d (where L is the system linear size and d is the dimensionality of the space); in the thermodynamic limit (L → ∞) the Lyapunov spectrum converges to an asymptotic curve, that is specific to each system and depends only on ρ = i/L d . The existence of such a limit shape is taken as the proof of extensivity, i.e. that quantities like the fractal dimension or the Kolmogorov-Sinai entropy are proportional to the system volume.
The LEs quantify the exponential expansion (contraction) growth rates along the N covariant or characteristic directions in the infinite-time limit. In fact, tangent space can be decomposed into covariant subspaces and the corresponding vectorbase provides relevant information that is encoded in the socalled characteristic/covariant Lyapunov vectors (LVs) [2][3][4]. These vectors actually contribute to identify the local structure of the invariant measure, to uncover the possible presence of collective phenomena [5], to determine the dimension of the inertial manifold [6], or to detect spurious LEs observed in embedded time series [7], to cite a few applications.
In a time interval t, an infinitesimal perturbation pointing along the i-th LV is expanded/contracted in tangent space by a factor e Λi(t) , where Λ i (t) is the so-called expansion rate. As a result of the variable degree of instability in phase-space, Λ i (t) fluctuates along the trajectory (or, equivalently, across phase-space). Nevertheless, the finite-time Lyapunov exponent (FTLE) λ i (t) = Λ i (t)/t converges, in the infinite-time limit, to the i-th LE, lim t→∞ λ i (t) = λ i , where brackets indicate an average over trajectories [8] (hereafter, phase-space average and time average are assumed to coincide).
It is natural to ask how the FTLE fluctuations scale with both time and system size. As already argued in [9], this question is not only connected with the convergence to the thermodynamic limit, but also with the extensivity of space-time chaos. The best way to approach the problem is by introducing the (time-dependent) variances, and the corresponding diffusion coefficients, It is expected that for large-enough times, the distribution P(λ, t) of FTLEs [λ(t) = (λ 1 (t), . . . , λ N (t))] is described by a suitable large-deviation function, P(λ, t) ∝ exp[−S(λ)t], and S( λ ) = 0. Under the fairly general assumption that the LEs fluctuations are short-range in time, (this is typically true away from bifurcations and phase transitions), the central limit theorem implies that S(λ) has a quadratic structure around its minimum, i.e. that P(λ, t) is Gaussian, where † denotes the transpose, while the matrix Q is the inverse of the symmetric diffusion matrix, i.e., Q = D −1 . Thus, within the Gaussian approximation, D describes how strong FTLE fluctuations are for a given system size. In particular, the diagonal elements D ii correspond to the diffusion coefficient of the expansion rates Λ i around the average growth λ i t.
Chaos extensivity would naively suggest that the diffusion coefficients should scale as L −d with the system size. Based on the numerical simulation of a variety of systems in d = 1 Kuptsov and Politi [9], have however found that D ii (L) ∼ L −γ with a wandering exponent that is smaller than 1: γ ≃ 0.85 for i > 1, while γ = 1/2 for i = 1. At the same time, it was found that the off-diagonal terms decay as L −1 and so do the eigenvalues of the matrix D, in agreement with the expectation for an extensive chaotic dynamics. It is rather intriguing that although the matrix D becomes increasingly diagonal in the thermodynamic limit, its diagonal terms scale differently from the matrix eigenvalues [9]. In the absence of theoretical arguments, one cannot a priori exclude that the scaling behavior of the bulk of the spectrum is affected by strong finite-size corrections, so that γ → 1 as L → ∞. We shed some light on this problem by resorting to the wellknown connection between LV dynamics and the kinetics of rough surfaces. This allows unveiling a theoretical connection between the scaling properties of the FTLEs fluctuations and the velocity fluctuations of rough surfaces subject to stochastic forces. This mapping leads to a scaling relation between the universal roughening exponents and the FTLE wandering exponent γ, which thus turns out to be a true critical exponent.
More specifically, in Section II, we develop the theoretical scaling arguments which lead us to derive the relevant mathematical relationships. In Section III, our systematic investigation of the maximal Lyapunov exponent confirms the theoretical expectations as well as the relationship with Kardar-Parisi-Zhang (KPZ) dynamics [18]. Section IV is devoted to the entirely new analysis of the bulk of the spectrum which suggests the correspondence with some yet unknown stochastic field theory. Finally, in Section V, we summarize the main results and briefly discuss the open problems.

II. THEORY
In this section we focus on the scaling behavior of the diagonal elements D ii of the matrix D. For the sake of simplicity, from now on, we drop the index i in the formulas and reintroduce it in the next sections, when it will be necessary to distinguish between different values of the integrated density ρ = i/L d .
In the following we show that D(L) ∼ L −γ , where the exponent γ can be expressed in terms of the scaling properties of the corresponding LV, v(x, t). Our scaling theory is based upon the well-known interpretation of the dynamics of a LV as the statistical evolution of a rough surface [10][11][12][13][14].
For each given LV v(x, t) we define an associated "surface" field h(x, t) through the logarithmic transformation h(x, t) = ln |v(x, t)|. The LV surface so defined is known to be generically rough and scale-invariant [10][11][12][13][14]. As we are interested in the expansion factor of the LV over a time interval t, it is useful to introduce the field which corresponds to the logarithm of the finite-time expansion factor (over a time t). The fluctuations around the average surface position at any given time are quantified by the surface width W , where the angular brackets denote an average over an ensemble of different trajectories, while the overline is a spatial average (here and in the following). Scale invariance generically leads to finite-size scaling of surface fluctuations that can be cast in the typical scaling form [15] where F (u) is a dynamical scaling function, which reaches asymptotically (u → ∞) a constant value and grows as u 2α/z for small values of u. All of the above means that at short times W 2 grows as t 2α/z , until t ∼ L z when W 2 saturates to a size-dependent value L 2α . The roughness exponent α and the dynamic exponent z quantify space-time correlations and fully characterize the statistical and dynamical behavior of the surface (equivalently, the LV.)

A. Main result
Before proceeding with the details of the theoretical derivation, we anticipate our main result, i.e. that for any FTLE its variance χ 2 (t, L) [defined by (1) with i = j] scales as where G(u ≫ 1) ≃ const. By then comparing this formula with Eq. (2), this implies that the LE diffusion coefficient scales as so that the wandering exponent γ heuristically observed in Ref. [9] is a truly critical exponent, connected to the LV surface roughening exponents, As a result, γ can be determined, once the roughening exponents of the corresponding surfaces are known. This is valid for any spatial dimension d.

B. Derivation of the scaling function G
We now derive Eq. (8) and give further details on the scaling of the diagonal elements of D in the intermediate regime before saturation. First, let us notice that the expansion rate Λ(t) necessarily refers to some norm in tangent space. The norm selection is irrelevant for the calculation of dynamical invariant quantities like the LEs or the diffusion coefficient D, as they involve an infinite-time limit. However, the finite-time expansion rates Λ(t) depend explicitly on the norm. Given a perturbation v(x, t), a rather broad family of q-norms can be defined as follows The standard Euclidean norm corresponds to q = 2. In the following, our theoretical arguments will be developed with reference to the 0-norm: v(x, t) 0 = x |v(x, t)| 1/L d , unless otherwise specified. The 0-norm is the most convenient and natural choice because, in this framework, computing the LV norm corresponds to determining the average height of the surface h(x, t) = ln |v(x, t)|, while the FTLE corresponds to the surface velocity, so that which, by definition (4), corresponds to the net displacement of the average surface position in a time interval t. Now, by combining Eqs. (1) and (11) one obtains This expression, like Eq. (5) for W 2 , is a quadratic correlation function of φ(x, t) and, therefore, it should scale as Consistency with Eq. (2) requires that G χ (u) diverges linearly with time. This information can be included in the above equation, by writing the scaling function as a product, Altogether, this leads to Eq. (7) and the main relation (9) for the scaling behavior of the asymptotic diffusion coefficient. Although our theoretical arguments are based on the use of the 0-norm, this does not affect the validity of our main scaling relation (9) as, in the infinite-time limit, the diffusion coefficient D is actually a dynamical invariant, independent of the norm used.
Before saturation sets in (t ≪ L z ), the effective diffusion coefficient is time-dependent and exhibits long-range temporal correlations. These correlations are important, since they correspond to an anomalous diffusion of the expansion rate Λ, as it can be inferred by looking at the "short-time" behavior of the dynamic scaling function G in Eq. (7), i.e. for 1 ≪ t ≪ L z . We indeed find that, still in the 0-norm framework, Let us now show how this result arises from a simple scaling argument. At short times, correlations only extend over a linear distance of order ∼ t 1/z so that the system can be considered as formed by a number is the local (intrablock) fluctuation. From Eq. (12), we can now estimate χ 2 so that for t ≪ L z . By then recalling that W 2 (t) ∼ t 2α/z in this time regime, we finally obtain Eq. (14). This concludes our scaling analysis.
It is worth to remark that, contrary to the asymptotic behavior, Eq. (14) is valid only with reference to the 0-norm. This is because, while the stationary diffusion coefficient implies an infinite-time limit, the time-dependent effective diffusion coefficient is defined for t < L z . This will become evident later on from the comparison with numerical calculations with both the 0-norm and the Euclidean norm.
Let us finally mention that the scaling function Eq. (14) arises also in the context of the nonequilibrium roughening, where it describes the velocity-fluctuations [16,17] of a driven interface.

C. Universality
It is well known that for a wide class of extended dynamical systems, which include coupled-map lattices, Lorenz-96 model, Kuramoto-Sivashinsky equation, and many others [11]-the LV surface h(x, t) associated with the first LV belongs to the universality class of KPZ [18]. The universality class of KPZ in extended dynamical systems is, therefore, very large and includes models which, in spite of relevant differences in the microscopic details, share a universal behavior of the first LV. This universality class seems to include all dissipative models with short-range interactions as well as some symplectic models. Therefore, we expect D 11 to be characterized by the same wandering exponent, for all models in the KPZ universality class. For instance, in d = 1 one can plug the exact values α = 1/2 and z = 3/2 in Eq. (9) to obtain γ = 1/2.
Regarding the wandering exponent in the bulk LEs (i ≫ 1), it is known that the corresponding LV surfaces are characterized by a different set of scaling exponents. This issue has been much investigated in the last few years. Numerical simulations in d = 1 suggest that the dynamic exponent is z ≃ 1 [13], and this result is also supported by theoretical arguments [14]. Moreover, 2α has been found to be much smaller than that in KPZ, lying in the range 0.15-0.2 [13,14]. Our main relation in Eq. (9) suggests a wandering exponent γ = 0.8-0.85 for the diffusion of the bulk LEs that is consistent with the numerical observations in Ref. [9].

III. THE LARGEST LYAPUNOV EXPONENT: NUMERICAL RESULTS
In this and the next section we compare our theoretical predictions with detailed numerical simulations of several systems. Here we focus on the diffusion coefficient of the first LE.

A. One-dimensional systems (d = 1)
The first model we analyze is a chain of Hénon maps, whose LE fluctuations have been recently studied in Ref. [9]. The model writes where Dx n ≡ x n−1 − 2x n + x n+1 is the discrete Laplacian operator, n = 1, . . . , L, and we select a = 1.4, b = 0.3, ǫ = 0.025. In Fig. 1 we plot L γ χ 2 /t , obtaining that for different system sizes the data indeed collapse for γ = 1/2 onto a dynamic scaling function G that follows Eq. (7) and the predicted asymptotes both above (t ≫ L z ) and below (1 ≪ t ≪ L z ) the crossover time. In particular, the very good data collapse observed at long times validates Eq. (9). γ = 1/2 was also observed in Ref. [9] where the Euclidean norm was used to measure vector metrics in tangent space, instead of the 0-norm considered in this paper. This confirms that the norm choice does not affect the stationary behavior. At shorter times, finite-size corrections are more sizable, but one can nevertheless appreciate an increasing quality of the data collapse with L. The initial slope increases with L and approaches the theoretical prediction ν = 1/3. Notice that this means that χ 2 initially grows as t 4/3 , i.e. the Lyapunov dynamics is superdiffusive in the intermediate regime before saturation of fluctuations.
We have also studied numerically the Lorenz-96 model [19] dy a time-continuous toy-model of the atmosphere that represents the value of a scalar variable on a mid-latitude. The data collapse shown in Fig. 2 confirms that the diffusion of the largest LE is well described by Eq. (7). The only difference with re-  spect to the previous model is that, as the arrow indicates, the curves converge towards the asymptotic shape from below. As a last example of a (pseudo) 1d system, we present our numerical results for a very different type of chaotic system: a model with delayed feedback. For many models of this type, the main LV scales as in typical 1d spatio-temporal chaotic systems [11,20], after identifying the delay T with the system size L. Note that one also must rescale the time axis t by a factor 1/T , as it is so for the Lyapunov spectrum [21]. More specifically, we have carried out simulations of the Mackey-Glass model [22], The results are plotted in Fig. 3. Apart from a relatively slow convergence they are again in agreement with our theoretical predictions.

B. Two-dimensional systems (d = 2)
It is very instructive to check the theoretical predictions in two-dimensional models since the KPZ critical exponents depend on the spatial dimension. Current numerical capabilities allow us to study a two-dimensional lattice composed by L 2 coupled logistic maps (on a torus geometry) (20) where f [u(t)] ≡ 4u(t) [1 − u(t)] and the sum is over the set N of nearest neighbors.
In two dimensions, only numerical estimates of the KPZ scaling exponents α and z are available. The best estimations [23] are α ≃ 0.387, and z ≃ 1.613, so that we predict γ ≃ 0.839 and ν ≃ 0.720 for the wandering exponent in (9) and the time exponent in (14), respectively. Numerical results for the dynamic scaling function are plotted in Fig. 4. The data collapse is excellent, confirming again the validity of our theoretical arguments. Notice that the LE fluctuations decay with the system size faster in 2d than in 1d since γ is larger. At short times one observes again the presence of strong finitesize corrections but one can nevertheless appreciate the predicted scaling behavior ∼ t 0.72 as the system size is increased.

IV. THE BULK OF THE LYAPUNOV SPECTRUM: NUMERICAL RESULTS
In this section, we investigate numerically the scaling behavior of D ii for i > 1, including the intermediate time regime before saturation, by revisiting the chain of Hénon maps and studying the Lorenz-96 model. It is worth recalling that, in the limit of large system-sizes, the LEs depend on the integrated density ρ = (i − 0.5)/L (in this section we limit ourselves to studying one-dimensional systems). As a consequence, a meaningful comparison of LVs for different system sizes must be made by selecting the index i which corresponds to the same density ρ. Since i is, by definition, an integer variable in the following we interpolate between the two nearest integers that correspond to the given ρ-value.

A. Chain of Hénon maps
In Fig. 5 we plot the results of simulations performed with the chain of Hénon maps for ρ = 0.25. This ρ-value is (i) sufficiently distant from the singularity at ρ = 0 to avoid crossover problems, and (ii) small enough to be computationally achievable in large systems [24]. We monitor the evolution of dχ 2 /dt rather than χ 2 /t. The two quantities would be equally valid as both obey the same scaling relation (7), but we prefer to use the former one, since it converges faster (i.e. for smaller values of t) to the asymptotic value. The solid curves in Fig. 5 correspond to simulations performed with the 0-norm for different system sizes. Altogether, the good data collapse confirms our scaling analysis, with γ ≃ 0.865 (close to the numerical value 0.85 measured in [9]) and z = 1 (as determined from direct LV studies [13,14]). Moreover, the initial growth agrees with the theoretical prediction (see the dashed line, whose slope is ν = (d − γ)/z ≃ 0.135).
For comparison, in Fig. 5 we plot also the results obtained by computing the FTLE fluctuations obtained from the standard Gram-Schmidt orthogonalization procedure (see the symbols). They do not only correspond to a different way of computing the LEs, but also to a different norm (namely, the Euclidean or 2-norm). One can see that the asymptotic value of the diffusion coefficient fully agrees with the previous results: this is consistent with the expectations that longtime LEs are independent of the norm adopted and this extends to their fluctuations too. The shape of the corresponding dynamic scaling-function is, however, very different for both metrics, as expected. In fact, the LEs exhibit a sub-diffusive transient rather than super-diffusive behavior if the Euclidean metric is used.
Besides estimating the exponent γ, we have directly determined 2α from the covariant LVs for ρ-values below 1, in order to test the validity of our main relation (9). The best estimation of α is typically obtained from the structure factor (power spectral density) of h, which follows a power-law decay due to its self-affine character, whereĥ is the Fourier transform. A general representation is portrayed in Fig. 6, where the effective value of 2α is plotted versus ρ for different system sizes. In the bulk (i.e. for 0 < ρ < 1), the data reveals a clear tendency to flatten towards 2α ≈ 0.16 for increasing the system-size. This α-value corresponds to γ ≈ 0.84, to be compared with the direct estimate γ = 0.865. This agreement, besides validating relation (9), hints at a possible universal behavior of the LV structure in the bulk of the spectrum. A careful numerical analysis (analogous to that described in Fig. 5) for ρ = 0.75 (data not shown) further confirms that γ is independent of ρ. Appreciable deviations from the basal value 2α ≈ 0.16 are clearly visible in Fig. 6 in the vicinity of the two extreme values ρ → 0 and ρ → 1. For ρ ≪ 1, we know that the first LV follows KPZ scaling, 2α = 1, and the data must show a crossover towards such a different scaling, when ρ = 0.5/L ≈ 0 is approached. A similar behavior is found for ρ = 1. Note that this value does not correspond to the smallest LE (which is obtained for ρ = 2, as there are 2L exponents in a system of size L), but rather to the edge of the band of positive LEs. In fact, for these parameter values there exists a gap between the positive and negative bands of LEs: the singularity for ρ = 1 thus reinforces the idea that band edges are characterized by a quantitatively different behavior, although here we do not see a direct way to map the evolution onto, e.g., KPZ dynamics.

B. Lorenz-96 model
The most intriguing message that arises from the study of the Hénon maps is the possibly universal scaling behavior of the bulk LVs. Given the relevance of such an observation, we have studied also the Lorenz-96 model (18). Being a continuous-time system, simulations are heavier than in the previous case and for this reason we have been able to carry out extensive simulations only for ρ = 0.1 which is nevertheless far enough from ρ = 0 to draw meaningful conclusions. The results for L up to 256 for the 0-norm and covariant LVs, and up to 512 for the Gram-Schmidt LVs, are plotted  in Fig. 7. The good data collapse confirms that the dynamic exponent is z = 1 as in the previous model. As for γ, we find a slightly different value, namely γ = 0.897 (to be compared with γ = 0.865). So far it is not possible to determine whether this difference is significative or just due to strong model-dependent finite-size effects hiding a universal systemindependent value. However, the closeness of both numbers suggests that γ is universal in the bulk Lyapunov spectrum.
We have also determined the values of α from the struc- ture factors S(k) of the LV-surfaces. The values of 2α for the positive LEs are plotted in Fig. 8. As shown above for the chain of Hénon maps, here the curve also becomes increasingly flat as the system size grows except for the points close to ρ = 0 and ρ ≃ 0.33: (i) ρ = 0 corresponds to the first LE, where we know that 2α = 1; (ii) ρ ≃ 0.33 corresponds to the vanishing Lyapunov exponent [14], so that S(k) ∼ k 0 , i.e. 2α = −1. These results suggest the existence of a common α for the bulk in the thermodynamic limit. Our best estimation is 2α ≃ 0.13, not far from the estimation 0.15 in [14]. Moreover, assuming z = 1, the value of γ expected via Eq. (9) is 0.87, which is relatively close to the value 0.897 observed in Fig. 7.

V. CONCLUSIONS AND OPEN PROBLEMS
Altogether, in this paper we have shown that the analogy between roughening phenomena and LV dynamics in spatially extended systems is rather fruitful in that it allows relating the scaling behavior of the LE diffusion coefficients with the roughness exponents of the corresponding vectors. Numerical simulations of various models support the general scaling relation (7) and of its asymptotic behavior (8) and (14). With reference to the first LV, our analysis confirms the validity of the relationship with the KPZ equation in dissipative systems.
More intriguing is the question of the scaling behavior for the bulk of the Lyapunov spectrum. Our studies of Hénon maps and of the Lorenz-96 model reveal that in both cases, independently of ρ, z = 1 and α ≈ 0.07. Nonetheless, small but not-so-negligible differences for the latter exponent are observed.
Two questions are, therefore, still open in the case of the fluctuations of the bulk LEs: (i) whether α is strictly larger than zero (γ < 1) in the thermodynamic limit; (ii) the very existence of a single universality class. The γ < 1 issue arises from the comparison with the cross-correlations of different LEs-namely, the scaling of the off-diagonal terms D ij with i = j. In Ref. [9] such terms have been found to scale as 1/L, i.e. their γ-value is 1. This finding has been interpreted as the evidence of an extensive behavior of the large deviation function (which is proportional to the inverse of D). On the one hand, it is strange that the diagonal terms scaling differs from that of the off-diagonal ones. On the other hand, the very fact that γ < 1 implies that α > 0. This, in turn, indicates a strong localization of the covariant LV, a property that has been observed in different models by using different methods [13,25]. On the basis of the results derived in this paper for the Hénon maps (that are more reliable than those for the Lorenz-96 model) we feel confident in stating that γ < 1. In order to draw firmer conclusions, however, we believe that it is necessary to make some substantial progress on the theoretical side by either identifying a minimal stochastic model of the LV dynamics for i ≫ 1 (such as the KPZ equation for the first vector), or by resorting to new results within the field of random matrices.