**2. Methods**

While fractal processes are characterized by long-term correlations that are described by a single scaling exponent, multifractal time series subsets with small and large fluctuations can scale differently, and the analysis of long-term correlations results in a hierarchy of scaling exponents [46]. Multifractal analysis of temporal series can be performed using different methods, such as the wavelet transform modulus maxima (WTMM) method [55], multifractal detrended fluctuation analysis (MF-DFA) method [46] and multifractal detrending moving average method (MF-DMA) [56]. In this work, we employ MF-DFA, which has been found to produce reliable results [57] and has been widely used to analyze physiological signals [58–60], geophysical data [61], weather data [62], and financial time series [63].

The implementation of the MF-DFA algorithm can be described as follows [46]:

i The first step is the integration of the original series *<sup>x</sup>*(*i*), *i* = 1, . . . , *N* to produce

$$X(k) = \sum\_{i=1}^{k} [x(i) - \langle x \rangle], \quad k = 1, \dots, N,\tag{1}$$

where *x* = 1 *N* ∑*k i*=1 *x*(*i*) is the average.

ii Next, the integrated series *X*(*k*) is divided into *Nn* = *int*(*N*/*n*) non-overlapping segments of a length *n*, and in each segmen<sup>t</sup> *ν* = 1, ... , *Nn*, the local trend *Xn*,*<sup>ν</sup>*(*k*) is estimated as a linear or higher order polynomial least square fit and subtracted from *<sup>X</sup>*(*k*).

iii The detrended variance

iv

$$F^2(n, \nu) = \frac{1}{n} \sum\_{k=(\nu-1)n+1}^{\nu n} \left[ X(k) - X\_{n,\nu}(k) \right]^2 \tag{2}$$

is calculated for each segmen<sup>t</sup> and then averaged over all segments to obtain the *q*th order fluctuation function:

$$F\_{\emptyset}(n) = \left\{ \frac{1}{N\_n} \sum\_{\nu=1}^{N\_n} [F^2(n,\nu)]^{q/2} \right\}^{1/q},\tag{3}$$

where, in general, *q* can take on any real value except zero.

Repeating this calculation for all box sizes provides the relationship between the fluctuation function *Fq*(*n*) and box size *n*. *Fq*(*n*) increases with *n* according to a power law *Fq*(*n*) ∼ *nh*(*q*) if long-term correlations are present. The scaling exponent *h*(*q*) is obtained as the slope of the linear regression of log *Fq*(*n*) versus log *n*.

The power law exponent *h*(*q*) is called the generalized Hurst exponent, where for stationary time series, *h*(2) is identical to the well-known Hurst exponent *H*. For positive *q* values, *h*(*q*) describes the scaling behavior of large fluctuations, while for negative *q* values, *h*(*q*) describes the scaling behavior of small fluctuations, while *h*(*q*) is independent of *q* for monofractal time series and a decreasing function of *q* for multifractal time series.

The generalized Hurst exponents are related to the Renyi exponents *τ*(*q*) defined by the standard partition function-based multifractal formalism *τ*(*q*) = *qh*(*q*) − 1. For the monofractal signals, *τ*(*q*) is a linear function of *q* (as *h*(*q*) = *const*.) and for multifractal signals *τ*(*q*) is a nonlinear function of *q*. A multifractal process can also be characterized by the singularity spectrum *f*(*α*), which is related to *τ*(*q*) through the Legendre transform:

$$
\alpha(q) = \frac{d\tau(q)}{dq},\tag{4}
$$

$$f(\alpha(q)) = q\alpha(q) - \tau(q),\tag{5}$$

where *f*(*α*) is the fractal dimension of the support of singularities in the measure with Lipschitz–Holder exponent *α*. The singularity spectrum of the monofractal signal is represented by a single point in the *f*(*α*) plane, whereas the multifractal process yields a single humped function.

Multifractal spectra reflect the level of complexity of the underlying stochastic process and can be characterized by a set of three parameters, which are determined as follows. The singularity spectra are fitted to a fourth degree polynomial:

$$f(\mathfrak{a}) = A + B(\mathfrak{a} - \mathfrak{a}\_0) + C(\mathfrak{a} - \mathfrak{a}\_0)^2 + D(\mathfrak{a} - \mathfrak{a}\_0)^3 + E(\mathfrak{a} - \mathfrak{a}\_0)^4 \tag{6}$$

The multifractal spectrum parameters are found as the position of the maximum *α*0 = arg max*α f*(*α*), the width of the spectrum *W* = *αmax* − *αmin* obtained from extrapolating the fitted curve to zero, and the skew parameter *r* = (*<sup>α</sup>max* − *<sup>α</sup>*0)/(*<sup>α</sup>*0 − *<sup>α</sup>min*), where *r* = 1 for symmetric shapes, *r* > 1 for right-skewed shapes and *r* < 1 for left-skewed shapes. These three parameters can be used to evaluate the complexity of the underlying process. A small value of *α*0 means that the process is correlated and more regular in appearance. The width *W* of the spectrum measures the degree of multifractality of the process, where a wider range of fractal exponents leads to "richer" structures. The skew parameter *r* indicates which fractal exponents are dominant: the *f*(*α*) spectrum is rightskewed (*r* > 1), and the process is characterized by a "fine structure" (small fluctuations) if high fractal exponents are dominant, whereas the process is more regular or smooth, the *f*(*α*) spectrum is left-skewed (*r* < 1), and the fractal exponents describe the scaling of large

fluctuations if the low fractal exponents are dominant. In summary, a signal with a high value of *α*0, a wide range *W* of fractal exponents (higher degree of multifractality) and a right-skewed shape (*r* > 1) may be considered more complex than one with the opposite characteristics [60].

The two sources of multifratality in a time series are (1) a broad probability density function for the values of the time series and (2) different long-term correlations for small and large fluctuations. The type of multifractal can be found by randomly shuffling the series and analyzing its behavior. For multifractals of the second type, the shuffled series exhibits simple random behavior (since long-term correlations are destroyed), and the width of the *f*(*α*) spectrum is reduced to a single point. For multifractals of the first type, the width of the *f*(*α*) spectrum remains the same (since the probability density cannot be removed), and for multifractals of types 1 and 2, the shuffled series shows weaker multifractality than the original series [46].

The time-dependent MF-DFA algorithm is based on the sliding window technique and yields a temporal evolution of multifractality in the system. Given a time series *x* = *x*1, ... , *xN*, many sliding windows *zt* = *<sup>x</sup>*1+*t*Δ, ... , *xw*+*t*Δ, *t* = 0, 1, ... , *N*−*w*Δ are constructed, where *w* ≤ *N* is the window size, Δ ≤ *w* is the sliding step and the operator [.] denotes taking the integer part of the argument. The values of the time series in each window *zt* are then used to calculate the multifractal spectrum at a given time *t* using the method described above. This allowed us to obtain time evolutions for the three complexity parameters.
