A Novel Intermittent Jumping Coupled Map Lattice Based on Multiple Chaotic Maps

: Coupled Map Lattice (CML) usually serves as a pseudo-random number generator for encrypting digital images. Based on our analysis, the existing CML-based systems still suffer from problems like limited parameter space and local chaotic behavior. In this paper, we propose a novel intermittent jumping CML system based on multiple chaotic maps. The intermittent jumping mechanism seeks to incorporate the multi-chaos, and to dynamically switch coupling states and coupling relations, varying with spatiotemporal indices. Extensive numerical simulations and comparative studies demonstrate that, compared with the existing CML-based systems, the proposed system has a larger parameter space, better chaotic behavior, and comparable computational complexity. These results highlight the potential of our proposal for deployment into an image cryptosystem. Author Contributions: Conceptualization, and methodology, and project administration, and X.L.; funding acquisition, R.H., and X.L.; software, A.D.; writing, validation,


Introduction
In 1985, Kaneko [1][2][3][4][5] formally proposed the Coupled Map Lattice (CML) consisting of discrete time steps, discrete spatial positions, and continuous state variables, and then deeply investigated its chaotic behavior. The series of studies [1][2][3][4][5] reveal that the CML system can create spatiotemporal chaos by combining local nonlinear dynamics and spatial diffusion. On the other hand, Fridrich [6] set forth that chaotic systems have remarkable cryptography properties, including simple structure, ergodicity, and high sensitivity to initial conditions and control parameters, and then constructed a chaos-based permutationand-diffusion cipher architecture for encrypting images. These pioneering achievements have been continuous sources that inspire subsequent methods  and have enriched the field of chaotic cryptography.
Zhang and Wang [7,8] proposed a Non-adjacent Coupled Map Lattice (NCML) in which the spatial positions of coupled lattices are dynamically determined through a nonlinear Arnold cat map. Then, the spatiotemporal chaotic sequences were used to drive a bit-level group permutation phase. Guo et al. [9] redesigned the NCML system by replacing a Logistic map with a piecewise-linear chaotic map, in the hope of boosting the turbulence of local reaction. To increase the spatial diffusion, reference [10] provided a hybrid-coupling mechanism that considers the adjacent and non-adjacent interactions simultaneously. Later, Zhang et al. [11] employed a 3D Arnold cat map to establish coupling relations between lattices, which enlarges the chaotic regime over a wide range of parameters. In [12], Liu et al. developed a NCML-based S-box shuffling method, which strengthens the ability to resist linear password attacks. Wang et al. [13] put forward a Dynamically Coupled Map Lattice (DCML), in which the coupling coefficient dynamically varies with the spatiotemporal indices. Then, Wang et al. [14] combined the ideas behind NCML and DCML, and invented a Non-adjacent Dynamically Coupled Map Lattice (NDCML). In [15], the spatiotemporal chaotic sequences generated by the NDCML system were used to govern a bidirectional image diffusion phase. Along this line of thought, Tao et al. [16] constructed a tailor-made dynamical coupling architecture, and claimed that it can evenly spread diffusion energies over lattices. Reference [17] described a customized globally CML whose coupling term is a superposition of multiple chaotic maps. Zhang et al. introduced a fractional order Logistic map into the CML system, and investigated the spatiotemporal dynamics conditioned on non-adjacent coupling [18] or time delay [19], respectively. Inspired by [19], Lv et al. [20] defined a dedicated delay function, and developed a new CML system with time delay. This CML system [20] serves as a pseudo-random number generator for encrypting RGB images in a channel-crossed manner [21]. Wang et al. [22,23] focused on a kind of cross CML whose diffusion mechanism proceeds along the temporal and spatial dimensions simultaneously. A dynamical coupling coefficient controlled by Tent map [22] or a pair of module operations [23] was introduced into the cross CML in order to enhance the chaotic behaviors. Reference [24] presented a mixed CML, which organizes three chaotic maps together guided by the spatial positions of lattices. In addition, some works [25][26][27][28][29] extended the original CML system to two dimensions, and realized the deployment in some real-world applications like image encryption or hiding.
In summary, the ideas of the above works [7][8][9][10][11]13,14,16,17,[22][23][24] mainly concentrate on three aspects: non-adjacent coupling [7][8][9][10][11], dynamical coupling [13,14,16,22,23], and multi-chaos mixing [17,24]. Based on our analysis, there still exist some problems in the three aspects. First, for the non-adjacent coupling, most of the previous works [7][8][9][10][11] merely adopted the Arnold cat map to determine the spatial positions of coupled lattices. However, we find that the Arnold cat map will repeatedly sample the same spatial position at a regular period. The repeated sampling usually tends to accumulate the diffusion energies heavily for a part of lattices, and thus incurs a non-uniform distribution of diffusion energies over lattices. Second, for the dynamical coupling, the previous works [13,14,16,22,23] simply endow the coupling coefficient with dynamics, but neglect the influence of coupling states on the spatiotemporal chaos. Third, for the multi-chaos mixing, the current works [17,24] fulfilled the mixing in a deterministic manner so that the mixing behavior can be clearly inferred according to the spatiotemporal indices. These problems not only limit the parameter space of spatiotemporal chaos, but also result in local chaotic behaviors, implying that some of the lattices may always lie in a non-chaotic regime. In such circumstances, extra computation or human labor may be required to carefully screen out the chaotic lattices before performing image encryption algorithms.
To alleviate the above problems, in this paper, we propose a novel intermittent jumping CML, abbreviated as IJCML, based on multiple chaotic maps. In our proposal, the pseudorandom information generated from multiple chaotic maps is integrated into the lattices in a complex nonlinear manner. The intermittent jumping mechanism is twofold. First, the intermittence refers to dynamically switching the coupling states between "off" and "on", which increases the complexity of diffusion. Second, the jumping establishes new non-adjacent coupling relations between lattices, which avoids the problem of repeated sampling so as to equalize the distribution of diffusion energies over lattices. We will describe the IJCML system at length in Section 2.
Our work studies the chaotic behaviors and nonlinear phenomena induced by the intermittent jumping mechanism, and observes the spatiotemporal chaos of the IJCML system from the following perspectives: diffusion energy analysis, power spectrum analysis, equilibrium degree analysis, information entropy analysis, inter-lattice independence analysis, Lyapunov exponent analysis, bifurcation diagram analysis, spatiotemporal behavior analysis, and computational complexity analysis. Note that these perspectives are widely chosen by [7][8][9][10][11]13,14,16,17,[22][23][24] for analyzing a dynamical system. We will conduct extensive numerical simulations and comparative studies in Section 3. The simulation results demonstrate that the proposed IJCML system is superior to or comparable to the existing ones [1,7,13,14] in terms of the aforementioned analyses, and thus highlight the potential for deployment into a practical image cryptosystem.

Intermittent Jumping Coupled Map Lattice
There are three essential elements in the CML system [1]: coupling mechanism, nonlinear mapping function, and multiple processing units (lattices), which forms a reactiondiffusion architecture. In each time step, the reaction procedure updates the state variable of a lattice via the nonlinear mapping function. At the same time, the diffusion procedure establishes a coupling relation between the current lattice and a set of other lattices. The traditional CML system can be expressed as: where n = 1, 2, · · · , N represent the time steps. The indices l = 1, 2, · · · , L stand for the spatial positions of lattices with a periodic boundary condition, and L is the size of the CML system. The symbol e denotes the coupling coefficient, whose value lies in the interval [0, 1]. In general, the Logistic map x k+1 = f (x k ) = µ · x k · (1 − x k ) serves as the nonlinear mapping function, where the parameter µ ∈ (0, 4], and x k represents the state variable after k iterations. Remarkably, when µ ∈ (3.57, 4], the Logistic map exhibits chaotic behaviors, including ergodicity and high sensitivity to initial conditions. See its bifurcation diagram and Lyapunov exponent in Figure 1a. Clearly, the coupling mechanism in Equation (1) is static, due to the constant coupling coefficient and the regularity of the coupling relations. To alleviate this problem, the NCML system [7] defines a non-adjacent coupling mechanism as follows: where the two spatial positions u(l) and v(l) are obtained from the Arnold cat map. That is: where s and t are the parameters of the Arnold cat map. However, as discussed before, the Arnold cat map used in Equation (3) belongs to a time-invariant transformation, and tends to sample the same spatial position repeatedly within a single time step.
The DCML system [13] provides another way to break the static defect, which can be formulated as: where the coupling coefficient is no longer a constant, but varies with the spatiotemporal indices l and n. To be specific, the DCML system [13] leverages an auxiliary Logistic map of the form: where k = 1, 2, · · · , LN, to generate the sequence of dynamical coupling coefficients. In practical use, this sequence is reshaped into a matrix of size L × N so as to gear towards the spatiotemporal indices in Equation (4). The auxiliary parameter µ aux is set to 3.99 for achieving outstanding dynamics, and the initial value of ε k is set to e. The NDCML system [14], which combines the non-adjacent coupling mechanism and the dynamical coupling coefficient together, can be viewed as a natural extension of the NCML [7] and DCML [13] systems. That is: x n+1 (l) = (1 − ε n (l)) · f (x n (l)) + (ε n (l)/ 2)· [ f (x n (u(l))) + f (x n (v(l)))], (6) where u(l) and v(l) are determined by Equation (3), and ε n (l) is obtained from the auxiliary Logistic map defined in Equation (5).
In this paper, we propose a novel Intermittent Jumping Coupled Map Lattice (IJCML), in which the pseudo-random information generated from the Logistic map, Sine map, and Chebyshev map will be integrated together in a dynamical manner. The definitions of the Sine map and the Chebyshev map can be described by the following equations: where γ and φ are the systems' parameters. As shown in Figure 1b, the Sine map shares similar chaotic behaviors with the Logistic map. As shown in Figure 1c, z k fills in the range [− 1, 1] , and the Lyapunov exponent increases monotonically with φ. The Chebyshev map starts to enter the chaotic regime when φ ≥ 1.
With these preparations, we formulate the IJCML system as: x n+1 (l) = (1 − e) · f (x n (l)) + (e/2)·[g(x n (u n (l))) + g(x n (v n (l)))], if w n (l) ≥ 0.5, f (x n (l)), otherwise, where w n (l) determines the couple state of the lth lattice at the nth time step. In this paper, the Chebyshev map is used to generate w n (l), taking the form w n (l) = (z n (l) + 1)/2, where z n (l) is the reshaped version of z k . The parameter φ is set to 3.999 so that w n (l) varies with the spatiotemporal indices dynamically. The intermittence is reflected in the phenomenon that a lattice sometimes interacts with other ones (when w n (l) ≥ 0.5), and other times updates alone (when w n (l) < 0.5). Clearly, the coupling states are mutually different between lattices.
When w n (l) ≥ 0.5, the spatial positions of coupled lattices are u n (l) and v n (l), respectively. In this paper, we abandon the Arnold cat map, and resort to a chaotic map to determine u n (l) and v n (l). Specifically, the pseudo-random information w n (l) is reused here in the following form: where the floor sign · rounds down to the nearest integer of the number enclosed within the sign. The sign taking its right operand, namely w n (l) · 10 7 in Equation (10), as input represents a compound S-box substitution operation. Specifically, the right operand is first converted to a six-digit hexadecimal integer. For example, if w n (l) · 10 7 = 1234567, its hexadecimal representation equals '12D687'. Second, the six-digit hexadecimal integer is divided into three groups, each of which contains two digits. The three groups of '12D687' are '12', 'D6', and '87', respectively. Third, apply substitution operation to each group in turn, where we use the S-box recommended in AES [30]. See Figure 2 for details. By doing so, we obtain 'C9', 'F6', and '17', respectively. Forth, concatenate the substitution values, yielding 'C9F617'. Convert this six-digit hexadecimal integer into decimal representation, yielding 13235735. As shown in Equation (10), the compound S-box substitution operation is only used in the course of calculating v n (l), and the nonlinearity of the S-box ensures that u n (l) = v n (l) holds. This is a straightforward way to avoid sampling the same spatial position for a specific spatiotemporal index. The pseudo-randomness of u n (l) and v n (l) inherits from w n (l), so that they vary with the spatiotemporal indices as well. The new non-adjacent coupling relations are time-varying, compared with the one in NCML [13] or NDCML [14]. In this paper, we use the term "jumping" to characterize the pseudorandomness and the spatiotemporal variability of the coupling relation.  (9) and (10), the proposed IJCML system combines the Logistic map, the Sine map, and the Chebyshev map together in a nonlinear complex way. First, the Logistic map serves as the nonlinear mapping function as usual. Second, the Sine map plays a key role in the coupling term. In this paper, the parameter of the Sine map is set to 0.999, so that, in each time step, the state variables of the lattices at the spatial positions u n (l) and v n (l) will be updated along chaotic trajectories. Third, the Chebyshev map generates the pseudo-random signal that is used to ensure the dynamic switching of coupling states and coupling relations. Additionally, as long as we protect the initial state variables of the Sine map and the Chebyshev map (e.g., treating them as keys), it is virtually impossible for an adversary to infer the mixing behaviors. This alleviates the problem of deterministic mixing in the existing works [17,24]. This work provides an intermittent jumping mechanism, which can effectively promote the turbulence evolution amongst lattices.

Analysis of Dynamic Behaviors
In this section, we analyze the chaotic behaviors of the IJCML system from nine aspects, and make comparisons with the baseline systems, including CML [1], NCML [7], DCML [13], and NDCML [14]. Throughout our analysis, the number of lattices L is set to 100, while the maximum number of iterations, namely N, equals 1000. Unless explicitly stated, the auxiliary parameters, including s, t, µ aux , are set to 23, 12, 3.99, respectively.

Diffusion Energy Analysis
The goal of this analysis is to check whether the distribution of diffusion energies between lattices is uniform or not. In this paper, the diffusion energy of the lth lattice is defined as: where δ(·) denotes the Dirac delta function. Equation (11) reveals that there are two factors determining the diffusion energy. One is the total number of times that the lth lattice is sampled after N iterations. Another is the square of the coupling coefficient. Note that, for a specific CML-based system, the spatiotemporal variability of u n (l ), v n (l ), and ε n (l ) in Equation (11) may be absent. For example, for the original CML system, we shall set that u n (l ) = l − 1, v n (l ) = l + 1, and ε n (l ) = e, respectively. Since this analysis is primarily concerned with whether the distribution is uniform or not, rather than the magnitude of E(l), we normalize the diffusion energy as follows: where l = 1, 2, · · · , L. It is desirable to equalize this distribution so that each lattice contributes equally to the diffusion of spatiotemporal dynamics. As demonstrated in [7,8], the non-adjacent coupling mechanism can strengthen the diffusion of spatiotemporal dynamics to some extent. Hence, in this paper, we define inter-lattice diffusion distances as follows: which measures how far the lth lattice's diffusion energy is spread to the others. Note that, for the IJCML system, the calculations of Equations (11) and (13) are triggered only when w n (l ) ≥ 0.5. It is desirable for larger inter-lattice diffusion distances with uniform distribution. Figure 3 exhibits results for the diffusion energy analysis. In this simulation study, we set e to 0.99. The blue histogram in the upper panel is the normalized distribution of diffusion energies, while the bottom panel shows the histogram of the inter-lattice diffusion distances. In addition, we print the standard deviation (std − ) of the blue histogram, and print the average (avg + ) and the standard deviation (std − ) of the orange one, where the superscript "+" (or "−") is intended to indicate that a higher (or lower) value is better. We find that although the CML system achieves the best uniformity with zero-valued standard deviations, D(l) equals 1 for all l, reflecting that each lattice only interacts with its adjacent counterpart in a regular way. As shown in Figure 3b, the inter-lattice diffusion distances reach larger values, which demonstrates that the NCML system indeed enlarges the range of coupling. However, there are several peaks that appear periodically in the histograms. This is due to the use of the Arnold cat map, which repeatedly samples the same spatial position at a regular period. Such non-uniform distribution implies that only several lattices may dominate the diffusion of spatiotemporal dynamics, leading to the local chaotic behaviors. As shown in Figure 3c, the dynamical coupling coefficient can slightly modulate the diffusion energies. However, it fails to equalize NDCML's histograms. In contrast, the proposed IJCML system equalizes the distribution of diffusion energies, and achieves larger inter-lattice diffusion distances at the same time. This validates that the IJCML system can uniformly spread the diffusion energies to non-adjacent lattices, so as to encourage the diffusion of spatiotemporal dynamics.

Power Spectrum Analysis
Power spectrum is a significant tool to visualize whether a sequence has chaotic characteristics or not. The Fourier theorem sets forth that any periodic signal can be expressed as Fourier series, corresponding to a discrete frequency spectrum. By contrast, an aperiodic signal is represented by Fourier integral with a continuous frequency spectrum. As such, if the power spectrum has one or more distinct peaks, the signal must be periodic or quasi-periodic. Conversely, if the power spectrum has a noise-like modality without obvious peaks, the signal is chaotic. Figure 4 shows results for the power spectrum analysis, in which we investigate twelve combinations of settings for the parameters µ and e. In this simulation study, the sequence with 1000 elements generated by 50th lattice is selected for analysis. In addition, for the sake of observation, we swap the left and right halves of a power spectrum, in the sense of shifting the zero-frequency component to the center. . Power spectrum analysis for the CML system, the NCML system, the DCML system, the NDCML system, and the proposed IJCML system. By comparison, we find that the parameter e has a relatively small influence on the power spectrum. Thus, the following discussion mainly focuses on the effect of different values of µ. As we see, for the baseline systems, almost all the power is concentrated at zero-frequency component when µ = 2.49 or 2.99. This means that there is almost no fluctuation in the corresponding sequence, so that it belongs to a direct-current signal. When µ = 3.49, the baseline systems distribute the power onto several frequency components. This implies that the corresponding sequence contains periodic transitions between several state variables, so that it belongs to a deterministic signal. When µ = 3.99, five power spectrums share the similar noise-like modality, meaning that all the systems have turned into the complete turbulence pattern [4]. Remarkably, even when µ = 2.49, 2.99, or 3.49, the IJCML system's power is evenly dispersed over all frequency components, in the sense that there exists only one small peak submerged in the noise-like fluctuations. This phenomenon demonstrates that the proposed system behaves better than the baseline ones, and indeed enlarges the parameter space of spatiotemporal chaos.

Equilibrium Degree Analysis
As the name implies, the equilibrium degree analysis is to evaluate whether the number of ones is close to the number of zeros in a binary sequence. Naturally, its definition can be written in the following form: where l = 1, 2, · · · , L. The notations P(l) and Q(l) record the number of ones and zeros, respectively, in the lth binary sequence, and we have N = P(l) + Q(l). The value of equilibrium degree lies in the interval [0, 1], and the higher the value, the better the equilibrium degree. For some bit-level image cryptosystems, the equilibrium degree analysis is helpful to check whether a CML-based system can generate pseudo-random binary sequences. Figure 5 shows results for the equilibrium degree analysis under the setting 2.8 ≤ µ < 4 and 0 ≤ e ≤ 1. In this simulation study, the sequence x n (l) generated by the lth lattice is first binarized through the operation mod( x n (l) · 10 7 , 2), where n = 1, 2, · · · , N. We calculate the average value of ED(l) over all lattices, where l = 1, 2, · · · , L. As we see, when µ < 2.99, the average values of equilibrium degrees are all equal to 0 for the baseline systems. When 2.99 ≤ µ ≤ 3.57, the baseline systems become trapped in an unstable transition zone, over which the average values are drastically varying. When µ > 3.57, the DCML system reaches a flat plateau with the average value close to 1 as shown in Figure 5c. For the CML system, however, there exists an obvious sunken area embedded into the flat plateau at the range 3.57 < µ < 3.99 and 0 < e < 0.18. For the NDCML system, two straight ravines occur at µ = 3.74 and 3.84, respectively, and directly go through the flat plateau. The NCML system has a wavelike plateau, which involves multiple sunken areas and ravines. This means that these baseline systems fail to achieve a stable equilibrium degree even when the nonlinear mapping function f (x) has entered into the chaotic regime (i.e., µ > 3.57). In contrast, due to the mixture of multiple chaotic maps, the proposed IJCML system possesses a more stable and larger flat plateau as shown in Figure 5e. If e is set to 0, the IJCML system's equilibrium degree degrades considerably, which validates that the intermittent jumping mechanism is helpful to enhance the chaotic behaviors. Based on these observations, we state that the IJCML system outperforms the baseline ones in terms of the equilibrium degree.

Information Entropy Analysis
Information entropy is a frequently-used criterion to measure the degree of confusion of a dynamical system. Following the quantization procedure previously applied in [13,14,16], we uniformly quantize a state variable x n (l) to an integer between 0 and 9. Specifically, the quantization procedure takes the form x n (l) · 10 . The information entropy of a specific lattice is defined by: where c = {c|c = 0, 1, · · · , 9} denotes the possible levels of quantization, and p(c) represents the probability of occurrence of c in a quantized sequence. Since we consider ten levels in the course of quantizing, the maximum value of H(c) equals log 2 10 ≈ 3.3219. It is desirable to achieve a higher value of H(c), which implies better pseudo-randomness of the quantized sequence. Figure 6 shows the results for the information entropy analysis. We conduct the simulation study under the setting 3 ≤ µ < 4 and 0 ≤ e ≤ 1, and calculate the average value of H(c) over all L lattices. We see that, for the baseline systems, the average values of H(c) form a stair-stepped upward trend as µ increases, and the two notable steps appear at µ = 3.02 and 3.48, respectively. When µ > 3.48, the DCML system has a smooth surface, in which the average values are independent of e. For the NDCML system, two straight ravines that appear at µ = 3.74 and 3.84, respectively, go through the surface. This phenomenon is consistent with the observation in Section 3.3, demonstrating that there still exist defects in the NDCML system even when µ is set to a higher value. For the CML system, a L-shaped valley with larger area is embedded into the surface. This analysis provides a warning that one should carefully select the parameters before using a CML-based image cryptosystem. Equipped with the non-adjacent coupling mechanism, the NCML system reduces the area of the valley, but fails to construct a smooth surface. For the IJCML system, the average value of H(c) increases along the axes of µ and e, and the leaf-like surface covers a larger parameter space of spatiotemporal chaos. These results verify that the IJCML system is superior to the baseline ones in terms of information entropy.

Inter-Lattice Independence Analysis
Some CML-based image cryptosystems may use multiple lattices simultaneously for driving the permutation and diffusion phases. In such circumstance, the higher the independence between two lattices, the stronger the ability to combat information leakage. In this paper, a Mutual-Information-based Entropy Distance (MIED) is defined to measure the independence between two quantized sequences. The mathematical formula is: where the superscript i (or j) stands for the lattice's index, while the symbol MAX is a constant, representing the maximum possible value of H(c i ) (or H(c j )). In this simulation study, we quantize a state variable into ten levels as before, so that MAX= log 2 10 ≈ 3.3219. In Equation (16), I(c i ; c j ) is the mutual information between c i and c j , and its definition is written as: The nonnegativity of mutual information ensures that 0 ≤ MIED(c i , c j ) ≤ MAX holds for any pair of lattices i and j. The upper and lower bounds, namely MAX and 0, correspond to the maximal and minimal degrees of independence, respectively. When µ = 3.49, the values of MIED, in general, form a middle-level plane containing fluctuations of various patterns. Interestingly, for the NCML system, the punctate leakages are expanded to box-like pits as e increases. This is because a greater coupling coefficient usually enhances the correlations between a pair of lattices i and j, thereby yielding a higher value of I(c i ; c j ). Moreover, as shown in the second column of Figure 7b, the box-like pits are arranged with a regular period. This arises from the fact that the Arnold cat map, which is used in the NCML system, forces some pairs of lattices to interact with each other repeatedly during the iterative procedure. Unfortunately, this defect is inherited by the NDCML system, as shown in Figure 7d.  When µ = 3.99, the values of MIED constitute a plateau lying close to MAX. The groove along the principal diagonal corresponds to the self-independence of a specific lattice, namely MIED(c i , c i ), where i = 1, 2, · · · , L. Clearly, a flat plateau with a narrow groove indicates good performance of the inter-lattice independence. For the CML system, increasing e reinforces the interactions between two lattices, and thus results in a wider groove and four sunken regions around the corners. A similar phenomenon also happens in the DCML system, implying that a dynamical coupling coefficient is insufficient to eliminate the coupling-caused correlations. For the NCML and NDCML systems, there exist multiple diagonal grooves embedded into the flat plateau due to the use of Arnold cat map. From the first two columns of Figure 7e, we find that the IJCML system can reach higher value of MIED when µ = 3.49 and e = 0.9. More significantly, the IJCML system forms the flattest plateau with a unit-width groove when µ = 3.99. These comparisons demonstrate that the IJCML system has a larger parameter space for generating the independent sequences randomly.
Further, we calculate the average values of MIED over all lattices via Equation (18) as follows: and investigate the inter-lattice independence under the setting 3.5 ≤ µ < 4 and 0 ≤ e ≤ 1.
The corresponding results are displayed in Figure 8. As we see, the CML and NCML systems share the same defect, in the sense that a cuboid-shaped valley appears at the range 0.12 < e < 0.18. This result suggests that the CML and NCML systems may be unsuitable for a multi-lattice-based image cryptosystem. It is clear that the DCML and NDCML systems are well-behaved when µ > 3.7, and are both insensitive to e. In most cases, the IJCML system achieves higher values than the baseline systems, especially when 3.5 < µ < 3.7 and 0.85 < e ≤ 1. This result reflects that the intermittent jumping mechanism can effectively compensate for the performance degeneration, and thus demonstrates that the IJCML system surpasses the baseline ones in terms of the inter-lattice independence.

Lyapunov Exponent Analysis
Lyapunov exponent [31] is a canonical metric in describing a dynamical system. It quantifies the average exponential rate of separation (or convergence) between two initially close trajectories along each direction of phase space. Mathematically, the definition of Lyapunov exponent is formulated as: where λ denotes the Lyapunov exponent of F(x), while i stands for the time step. Typically, a positive Lyapunov exponent indicates that F(x) has chaotic behaviors with the trajectories diverging exponentially. Kolmogorov-Sinai Entropy Density (KSED), which incorporates the Lyapunov exponents of all lattices, is devoted to measuring the overall dynamics of a CML-based system. Formally, its definition is written as: where λ + (l) represents the positive Lyapunov exponent of the lth lattice. In other words, Equation (20) overlooks all negative Lyapunov exponents during its calculation.
In addition, Kolmogorov-Sinai Entropy Breadth (KSEB) [32] is used, in this paper, to count the proportion of chaotic lattices. We calculate its value as follows: where L + denotes the number of lattices whose Lyapunov exponents are positive. Clearly, the higher the values of KSED and KSEB, the better chaotic behaviors of a dynamical system. Figure 9 exhibits the resulting values of KSED under the setting 3 ≤ µ < 4 and 0 ≤ e ≤ 1. Starting from µ = 3.57, the DCML system possesses a smooth surface with positive curvatures, validating that increasing µ can constantly strengthen the chaotic behaviors. For the NDCML system, however, there exist two objectionable ravines appearing at µ = 3.74 and 3.84, and straightly passing through the surface. The CML and NCML systems behave unstably when µ > 3.57, as the values of KSED are varying sharply and a distinct valley appears at the range 3.57 < µ < 3.99 and 0.11 < e < 0.16. The IJCML system achieves numerous high values of KSED that form a flat plateau. Noticeably, even when µ < 3.57, setting e to a higher value can still lead the system into the chaotic regime. Thus, the proposed system has a larger parameter space of spatiotemporal chaos, which benefits from the intermittent jumping mechanism. In Figure 10, we display the resulting values of KSEB. The maximum value of KSEB, being equal to 1.0, means that all the lattices in the system are activated into spatiotemporal chaos. At first glance, the IJCML system's flat plateau, as shown in Figure 10e, is of the largest area. We note that, among the baseline systems, only the DCML system has a stable flat plateau when µ > 3.57. For the other baseline systems, there exist deep valleys or straight ravines embedded into the flat plateau. In other words, even though µ is set to a higher value, some lattices are still far from spatiotemporal chaos. Therefore, one shall spend time inspecting the dynamics carefully for each lattice before deploying them in an image cryptosystem. Further, we design a statistical tool called Cumulative Percentage Function (CPF) for analyzing KSEB. The CPF is to reflect the percentage of parameter pairs {µ, e} making KSEB less than or equal to a specific threshold α. We show the plots of CPFs in Figure 11. The ideally optimal case is that all possible combinations of µ and e can lead a system into the fully chaotic regime, where all its lattices have the chaotic behaviors, namely KSEB = 1.0. Obviously, the ideally optimal case corresponds to an impulse function centered at 1.0. Among all CPFs, the pink one is closest to the impulse function, consisting of a wide flat region and a narrow step region. The higher the step occurs at α = 1.0 in Figure 11, the larger the area of flat plateau in Figure 10. Based on above observations, we state that the proposed IJCML system is better than the baseline ones in terms of the Lyapunov exponent analysis.

Bifurcation Diagram Analysis
The bifurcation diagram, which is a plot of the state variables versus the control parameter(s), is commonly used to analyze a dynamical system. Let all the systems, in this simulation study, work under the setting 3.4 ≤ µ < 4 and e = {0.1 , 0.5, 0.9}. In Figure 12, we depict the bifurcation diagrams, in which the state variables are harvested from the 1st, 50th, and 100th lattices. For convenient comparison, the bifurcation diagrams belonging to different systems are coated with different colors.   For the baseline systems, we can clearly observe that period-doubling bifurcations form a route to spatiotemporal chaos. For example, the NDCML system gets trapped in period-two oscillations when µ = 3.4. Two bifurcation points at µ = 3.45 double the periods of the orbits and give rise to pitchfork bifurcations. This process is repeated as µ increases. Starting from µ = 3.57, it is virtually impossible to observe the pitchfork bifurcations because the number of bifurcation points becomes large and the gaps between the bifurcation points are negligible. Unfortunately, the chaotic oscillations of the baseline systems may occasionally be interspersed with periodic windows in some cases. For example, when setting µ to 3.63, we can observe the periodic windows in the NDCML system, where the oscillations suddenly degenerate into several periods. Comparing the baseline systems' bifurcation diagrams, we find that they roughly share the same phe- For the baseline systems, we can clearly observe that period-doubling bifurcations form a route to spatiotemporal chaos. For example, the NDCML system gets trapped in period-two oscillations when µ = 3.4. Two bifurcation points at µ = 3.45 double the periods of the orbits and give rise to pitchfork bifurcations. This process is repeated as µ increases. Starting from µ = 3.57, it is virtually impossible to observe the pitchfork bifurcations because the number of bifurcation points becomes large and the gaps between the bifurcation points are negligible. Unfortunately, the chaotic oscillations of the baseline systems may occasionally be interspersed with periodic windows in some cases. For example, when setting µ to 3.63, we can observe the periodic windows in the NDCML system, where the oscillations suddenly degenerate into several periods. Comparing the baseline systems' bifurcation diagrams, we find that they roughly share the same phenomenon. That is, the period-doubling bifurcations trigger chaotic oscillations interspersed with periodic windows. By contract, the proposed IJCML system has a large number of bifurcation points with tiny gaps between them even when µ = 3.4. This demonstrates that the IJCML system enters the chaotic regime early. Additionally, the IJCML system greatly reduces the periodic windows compared with the baseline ones. As shown in Figure 12e, a higher coupling coefficient introduces more pseudo-random motions from the multiple chaotic maps, so as to further enhance the chaotic oscillations. Based on these comparisons, we claim that the IJCML system behaves better than the baseline ones in terms of the bifurcation diagram analysis.

Spatiotemporal Behavior Analysis
Spatiotemporal behavior is a means of testing diffusions between lattices. In this simulation study, e is fixed at 0.5 for simplicity, while µ is set to 3.15, 3.57, and 3.99, respectively, for adjusting the spatiotemporal behaviors.
In Figure 13, we show space-amplitude plots, which are used to analyze the spatiotemporal behaviors. Similarly, the space-amplitude plots belonging to different systems are coated with different colors. In Figure 13, we show space-amplitude plots, which are used to analyze the spatiotemporal behaviors. Similarly, the space-amplitude plots belonging to different systems are coated with different colors.
(a) Results for the CML system (b) Results for the NCML system (c) Results for the DCML system   When µ = 3.15, the CML system simply creates period-two responses. The other baseline systems twist the period-two responses at a regular pace, resulting in a X-shaped pattern. The twisting points correspond to stable solutions of a system, whose positions depend on the initial state variables. In contrast, the X-shaped pattern is unrecognizable for the IJCML system. These results reveal that, when the baseline systems are restricted to the frozen random pattern [4], the IJCML system has advanced into the defect chaotic diffusion pattern [4].
When µ = 3.57, we can observe the period-doubling behaviors from the space-amplitude plots of the baseline systems, in the sense that the long-thin X-shaped pattern evolves into a compound version. Specifically, there exist multi-period responses twisted alternately, and the dynamic range of the solutions has been enlarged at the same time. In particular, some of the lattices in the DCML system, for example, the 40th one, have entered the defect chaotic diffusion pattern [4]. These results reveal that the baseline systems are transitioning from the frozen random pattern [4] to the defect chaotic diffusion pattern [4].
When µ = 3.99, the period-changing behaviors become complex and unstable, resulting in extremely twisted and superimposed responses. Remarkably, the dynamic range is further extended towards the upper and lower bounds. These results reveal that all the systems have entered the complete turbulence pattern [4], while the IJCML system has the strongest ability to approach the lower bound. Consequently, we state that the IJCML system has better spatiotemporal behaviors than the baseline ones. When µ = 3.15, the CML system simply creates period-two responses. The other baseline systems twist the period-two responses at a regular pace, resulting in a X-shaped pattern. The twisting points correspond to stable solutions of a system, whose positions depend on the initial state variables. In contrast, the X-shaped pattern is unrecognizable for the IJCML system. These results reveal that, when the baseline systems are restricted to the frozen random pattern [4], the IJCML system has advanced into the defect chaotic diffusion pattern [4].

Computational Complexity Analysis
When µ = 3.57, we can observe the period-doubling behaviors from the spaceamplitude plots of the baseline systems, in the sense that the long-thin X-shaped pattern evolves into a compound version. Specifically, there exist multi-period responses twisted alternately, and the dynamic range of the solutions has been enlarged at the same time. In particular, some of the lattices in the DCML system, for example, the 40th one, have entered the defect chaotic diffusion pattern [4]. These results reveal that the baseline systems are transitioning from the frozen random pattern [4] to the defect chaotic diffusion pattern [4].
When µ = 3.00, the period-changing behaviors become complex and unstable, resulting in extremely twisted and superimposed responses. Remarkably, the dynamic range is further extended towards the upper and lower bounds. These results reveal that all the systems have entered the complete turbulence pattern [4], while the IJCML system has the strongest ability to approach the lower bound. Consequently, we state that the IJCML system has better spatiotemporal behaviors than the baseline ones.

Computational Complexity Analysis
Computational complexity reflects the time and space consumptions while executing a target algorithm so as to provide a guidance for developing practical software.
First, we theoretically analyze the time complexity. Like the baseline systems [1,7,13,14], the proposed IJCML system also has two nested "for" loops in its calculation procedure. The outer "for" loop runs the dynamical system by N iterations, while the inner one traverses the L lattices. Hence, all the five CML-based systems investigated in this paper have the quadratic time complexity, namely O(NL).
Second, we theoretically analyze the space complexity. Clearly, it is necessary to save the resulting pseudo-random matrix of size N × L. This requires quadratic space complexity, namely O(NL). Compared with the original CML system [1], the variants, i.e., the NCML system [7], the DCML system [13], the NDCML system [14], and the proposed IJCML system need to occupy additional space to save the intermediate variables obtained from the auxiliary chaotic maps. See Equations (3), (5), and (8) for details. Fortunately, since these auxiliary chaotic maps can be separately executed out of the two nested "for" loops, the final space complexity remains the same, namely O(NL).
However, it is undeniable that the variants may require more time and space consumptions in practice. To make an intuitive comparison, we experimentally count the average time for running 1000 iterations over 100 lattices, and inspect the size of the additional space. The parameters µ and e are set to [3, 3.99] and [0, 0.99], respectively, and both of them take 0.01 as step size. Our computing device is a desktop computer with a 2.90 GHz Intel i7-10700 central processing unit, 16.00 GB memory. Our programming environment is Matlab R2017a installed on the Window 10 operation system. The intermediate variables are saved as ".mat" format. Table 1 lists the simulation results. The IJCML system takes 0.5400 s, on average, to generate a pseudo-random matrix of size 1000 × 100, which is faster than the NCML system [7], the DCML system [13], and the NDCML system [14], and is comparable to the original CML system [1]. High efficiency of the IJCML system is due to the intermittent jumping mechanism because the coupling operation, as shown in Equation (9), will be skipped with a probability of about 50%. On the other hand, the NCML system [7] only requires 3.42 KB space to additionally save the time-invariant spatial positions. In contrast, the remaining three variants prepare about 746.00 KB space to accommodate the dynamical coupling coefficients [13,14] or the intermediate variables w n (l) obtained from the Chebyshev map. Fortunately, a space of size 746.00 KB is almost negligible for modern hardware devices. Consequently, the computational complexity of the proposed IJCML is better than or comparable to the baseline systems [1,7,13,14].

Conclusions
In this paper, we propose a novel IJCML system based on multiple chaotic maps. The intermittent jumping mechanism establishes a new coupling mode, in which not only the coupling states but also the coupling relations dynamically vary with the spatiotemporal indices. We conduct extensive numerical simulations and comparative studies to analyze the proposed system from the following aspects: the diffusion energy analysis, the power spectrum analysis, the equilibrium degree analysis, the information entropy analysis, the inter-lattice independence analysis, the Lyapunov exponent analysis, the bifurcation diagram analysis, the spatiotemporal behavior analysis, and the computational complexity analysis. The simulation results adequately demonstrate that, compared with the baseline systems [1,7,13,14], the proposed IJCML system has better chaotic behaviors, which brings stronger spatiotemporal chaos for a single lattice and higher independency between lattices. Our future work is to study the effective way for discretizing the IJCML system, and to consider the issues of practical realization and circuit implementation [33].