Intermittent chaos in the Lennard-Jones oscillator

Abstract

The analysis of the classical Lennard-Jones (LJ) oscillator under harmonic excitation reveals a fascinating phenomenon—a rich interplay between quasi-periodic motion and intermittent chaos. It is shown that chaotic oscillations can emerge even at low forcing amplitudes, while the system remains largely non-chaotic across a broad parameter range. The emergence of strange attractors composed of densely organized regular patterns represents another distinctive and persistent feature of the LJ oscillator. The analysis employs bifurcation diagrams, the largest Lyapunov exponents (LLEs), correlation dimension (CD), Shannon entropy (SE), and other tools. A comparison with the Duffing oscillator reveals a substantial discrepancy in oscillation regimes across the whole range of forcing amplitudes. To ensure numerical stability, the variable-step Adams–Bashforth–Moulton (ABM) method in an explicit–implicit formulation was used.

Keywords

Lennard-Jones oscillator Duffing oscillator bifurcation diagram largest Lyapunov exponent fractal measure chaos

1. Introduction

1.1. Overview

1.1.1. Lennard-Jones potential

The most general form of the Lennard-Jones (LJ) potential W introduced by Jones [1], can be written as

W (r) = \frac{A}{r^{n}} - \frac{B}{r^{m}},

(1)

where r is the distance between neighboring atoms or particles; A and B are the dimensional coefficients; and $n > m > 0$ are the integers. The repulsion force at smaller distances and attraction force at larger distances are given by the following equation [2,3]

F (r) \equiv - \partial_{r} W (r) = - (\frac{nA}{r^{n + 1}} - \frac{mB}{r^{m + 1}}),

(2)

from where, the equilibrium distance at which $F (r) = 0$ becomes

r_{eq} = {(\frac{nA}{mB})}^{n - m} .

(3)

Herein, positive force corresponds to repulsion and negative to attraction. Potential (1) is also known as the Mie potential [4], being introduced in the context of materials science [5].

Another form of potential (1) proposed in [6], corresponds to $n = 12; m = 6$ and condition $nA = mB$ , yielding

W (r) = 4 e ({(\frac{d}{r})}^{12} - {(\frac{d}{r})}^{6}),

(4)

where e determines the depth of the potential well, and d is the dimensional distance at which the potential vanishes. Both $n = 12$ and $m = 6$ in equation (4) were chosen to match the experimental data on the intermolecular interactions [5,7]. Analogously to equation (2), the repulsion force for potential (4) becomes

F (r) = - 24 e d^{6} \frac{2 d^{6} - r^{6}}{r^{13}},

(5)

yielding the equilibrium distance

r_{eq} = 2^{1 / 6} d,

(6)

at which $F (r_{eq}) = 0$ , and hence, the potential well occurs.

1.1.2. LJ oscillator

Quite a large number of works relate to the analysis of oscillations in dynamical systems with restoring force defined by equation (5); see [8 –15]. In the one-dimensional case, the deflection x from the equilibrium interatomic or interparticle distance $r_{eq}$ is related to the current distance, as

r = r_{eq} + x,

(7)

Introducing dimensionless deflection

x^{'} = \frac{x}{r_{eq}}

(8)

and taking into account relation (6), the restoring force becomes [16]

F (x^{'}) = - α \frac{1 - {(1 + x^{'})}^{6}}{{(1 + x^{'})}^{13}},

(9)

where

α = \frac{3 e}{2^{1 / 6} d^{13}} .

(10)

In the latter form, the restoring force for the LJ potential can be determined by a single constant α [16]). Bearing in mind equation (9), the harmonic LJ oscillator can be written as [11,16,17]

r_{eq} m \partial_{tt}^{2} x^{'} = - F (x^{'}) - r_{eq} η \partial_{t} x^{'} + f \cos (ω t) .

(11)

where m is the oscillating mass; η is the linear viscous friction coefficient; f is the forcing amplitude; ω is the angular frequency; and t is the time. Hereafter, the prime notation on the dimensionless variable $x^{'}$ will be omitted for clarity.

Following Ueda [18 –21], by introducing a dimensionless time variable

τ = ω t,

(12)

Equation (11) can be rewritten:

\partial_{τ τ}^{2} x^{'} = - F^{'} (x^{'}) - η^{'} \partial_{τ} x + f^{'} \cos (τ),

(13)

where

F^{'} (x^{'}) = - α^{'} \frac{1 - {(1 + x^{'})}^{6}}{{(1 + x^{'})}^{13}}; η^{'} = \frac{η}{m ω}; f^{'} = \frac{f}{r_{eq} m ω^{2}}; α^{'} = \frac{α}{r_{eq} m ω^{2}} .

(14)

Equation (13) serves as the starting point for further analysis. Hereafter, the prime notation and subscript τ -indices, denoting differentiation with respect to τ, will be omitted for clarity.

Using molecular dynamics methods to analyze multimass chains of oscillators (11) and similar ones, it was found that chaotic oscillations can occur for specific combinations of elastic, thermodynamic, or elastic-magnetic parameters [22 –26]. In most of this research, the occurrence of strange attractors and Feigenbaum cascades, associated with period doubling in the vicinity of chaotic regimes, was analyzed, primarily in dynamical systems comprising multiple LJ oscillators [9,22,25,27 –32].

1.2. Problem statement

The analysis of LJ harmonic oscillator (13) reveals several notable phenomena associated with the occurrence of intermittent chaos interlacing with quasi-periodic motion. The principal finding of oscillations in the considered dynamical system reveals a fascinating phenomenon: the absence of chaos across a sufficiently wide range of forcing amplitudes. Instead, chaos emerges at low amplitudes and displays strange attractors consisting of dense sets of regular patterns. Notably, these observations indicate that even small (dimensionless) forcing amplitudes can induce chaotic oscillations.

The analysis employs bifurcation diagrams, the largest Lyapunov exponents (LLEs), correlation dimension (CD) measures, Shannon entropy (SE), and other tools used to identify chaos. To ensure numerical stability, the variable-step Adams–Bashforth–Moulton (ABM) method in an explicit–implicit formulation was used [33,34].

The emergence of chaotic oscillations and the corresponding strange attractors observed in the LJ oscillator at small forcing amplitudes are crucial for understanding the chaotic nature of interatomic and intermolecular interactions in systems governed by the LJ potential. These phenomena could also find applications in studying energy transfer mechanisms in chemical and physical systems.

2. Governing equations

2.1. Equations of motion

2.1.1. Non-autonomous system

Equation (13) can be expressed in the form of a first-order non-autonomous dynamical system:

\partial_{τ} X = G (X) + H (τ),

(15)

where $v = \partial_{τ} x$ and

X = (\begin{matrix} x \\ v \end{matrix}); G (X) = (\begin{matrix} v \\ - F (x) - η v \end{matrix}); H (τ) = (\begin{matrix} 0 \\ f \cos (τ) \end{matrix}) .

(16)

In equations (15) and (16), the prime notation is omitted.

Taking the gradient of G yields the Jacobian matrix [35]

J (X) \equiv \nabla_{X} G (X) = (\begin{matrix} 0 & 1 \\ - E (x) & - η \end{matrix}) .

(17)

where $E (x)$ is the tangent stiffness [36]:

E (x) \equiv \frac{dF (x)}{dx} = α \frac{13 - 7 {(1 + x)}^{6}}{{(1 + x)}^{14}} .

(18)

The condition for $E (x) = 0$ is known as phonon instability [37,38], which leads to the appearance of imaginary frequencies. To avoid it, the condition $E (x) > 0$ will be assumed to hold, implying

x < {(\frac{13}{7})}^{1 / 6} - 1 .

(19)

The eigenvalues of matrix $J (X)$ are as follows

e_{1, 2} = \frac{- η \pm \sqrt{η^{2} - 4 E (x)}}{2} .

(20)

2.1.2. Autonomous system

To analyze stiffness, the reduction of equation (15) to the autonomous system is needed. Introduce a new 3-vector

Y = (\begin{matrix} x \\ v \\ w \end{matrix}),

(21)

where $w = \cos (τ)$ , then equation (15) can be reduced to the following autonomous system

\partial_{τ} Y = K (Y) .

(22)

Herein,

K (Y) = (\begin{matrix} v \\ - F (x) - η v + fw \\ \pm \sqrt{1 - w^{2}} \end{matrix}) .

(23)

Taking the gradient of equation (23) with respect to Y yields the Jacobian matrix

J \equiv \nabla_{Y} K = (\begin{matrix} 0 & 1 & 0 \\ - E (x) & - η & 1 \\ 0 & 0 & \pm w {(1 - w^{2})}^{- 1 / 2} \end{matrix}) .

(24)

The eigenvalues of matrix (24) are as follows

e_{1, 2} = \frac{- η \pm \sqrt{η^{2} - 4 E (x)}}{2}; e_{3} = \pm w {(1 - w^{2})}^{- 1 / 2} .

(25)

Expressions (26) imply

e_{3} = {\begin{matrix} 0, w = 0 \\ \mp \infty, w = \pm 1 \end{matrix} .

(26)

Now, following [33,39], introduce the stiffness ratio:

r_{s} = max \frac{| Re (e_{i}) |}{| Re (e_{j}) |}, i \neq j .

(27)

Condition (26) implies

r_{s} \to \infty at {\begin{matrix} e_{3} = 0 & & e_{1}, e_{2} \neq 0 \\ e_{3} = \pm \infty & & \forall e_{1}, e_{2} \end{matrix} .

(28)

Thus, the considered autonomous dynamical system is stiff, making the use of explicit solvers of the Runge–Kutta type inapplicable. Stiffness of the considered dynamical systems was overlooked in most previous studies [9,19 –21,27 –32].

2.2. ABM method

To integrate the autonomous dynamical system (22), the explicit–implicit ABM variable step-size method is employed [34,40 –44].

2.2.1. Predictor (Adams–Bashforth explicit step)

The explicit step of the considered method relies on N previous steps:

Y_{n + 1} = Y_{n} + h (\sum_{k = 1}^{N} a_{k} (K (Y_{n + 1 - k}))),

(29)

where

Y_{k} = Y (t_{k}); h = t_{k + 1} - t_{k},

(30)

and intergrating coefficients $a_{k}$ depend on the number of steps N in the integrating scheme. The coefficients $a_{k}$ depend on the order of the Lagrange interpolation polynomial [34,45 –47]:

\begin{matrix} N = 1 : a_{1} = 1 \\ N = 2 : a_{1} = \frac{3}{2}; a_{2} = - \frac{1}{2} \\ N = 3 : a_{1} = \frac{23}{12}; a_{2} = - \frac{16}{12}; a_{3} = \frac{5}{12} . \end{matrix}

(31)

2.2.2. Corrector (Adams–Moulton implicit step)

The implicit step relies on the following equation [40]

Y_{n + 1} = Y_{n} + h (\sum_{k = 0}^{N} b_{k} (K (Y_{n + 1 - k}))) .

(32)

Now, the integrating coefficients become [40]

\begin{matrix} N = 1 : b_{0} = \frac{1}{2}; b_{1} = \frac{1}{2} \\ N = 2 : b_{0} = \frac{5}{12}; b_{1} = \frac{8}{2}; b_{2} = - \frac{1}{12} \\ N = 3 : b_{0} = \frac{9}{24}; b_{1} = \frac{19}{24}; b_{2} = - \frac{5}{24}; b_{3} = \frac{1}{24} \\ N = 4 : b_{0} = \frac{251}{720}; b_{1} = \frac{646}{720}; b_{2} = - \frac{264}{720}; b_{3} = \frac{106}{720}; b_{4} = - \frac{19}{720} . \end{matrix}

(33)

These coefficients, as in the previous case, depend upon the order of the Lagrange interpolation polynomial.

2.3. Chaos characterization

To investigate the onset and evolution of chaotic behavior, Poincaré sections, bifurcation diagrams, and the following fractal measures are utilized.

2.3.1. Poincaré sections and bifurcation diagrams

The Poincaré section, known also as Poincaré map, can be defined as a set of points in the phase space, corresponding to the fixed phase [48]

C (s) : ⋃_{t > Λ / 2} (x (mod (t, T); v (mod (t, T))),

(34)

where Λ is the simulation time; $mod (t, T)$ is the remainder after division of t by the period T of the harmonic driving force

T = 2 π ω^{- 1} .

(35)

Following [49], the first half part of the simulation interval $t \in (0; Λ / 2)$ is discarded to eliminate data points associated with unstable oscillations.

Taking the projection of the right-hand side in (34) onto the first or second coordinate and varying the bifurcation parameter yields the desired bifurcation diagrams [49,50]:

D_{x} (s) : \Pr_{x} C (s); D_{v} (s) : \Pr_{v} C (s),

(36)

where s is the bifurcation parameter. In the following analysis, the amplitude f of the driving force plays the role of the bifurcation parameter.

2.3.2. LLE

The LLE quantifies the average exponential rate at which nearby trajectories in a system’s phase space diverge [51]. A positive LLE indicates sensitivity to initial conditions—a defining property of chaotic dynamics—whereas a negative or zero exponent corresponds to stable or quasi-periodic behavior, respectively. The LLE is therefore widely regarded as a primary diagnostic for chaos detection. Conceptually, if two trajectories start arbitrarily close to each other, their separation $δ (t)$ evolves approximately as

δ (t) : δ (0) e^{λ t},

(37)

where λ is the Lyapunov exponent. When λ > 0, the separation grows exponentially, indicating that small uncertainties in the initial state amplify exponentially with time—a hallmark of chaos. Conversely, $λ < 0$ corresponds to stable or convergent motion, while $λ : 0$ is typically associated with quasi-periodic behavior [52]. The largest exponent in the Lyapunov spectrum, the LLE, therefore serves as a concise numerical indicator distinguishing chaotic dynamics from regular motion.

In applied settings, the LLE can be estimated either analytically from the linearized system equations or numerically from time-series data through phase-space reconstruction. Among several approaches, the algorithms [53 –55] offer a computationally efficient method based on the average divergence of neighboring points along the same trajectory, reconstructed in an embedding space using appropriate delay and embedding dimension parameters; see also [56 –59]. Unlike earlier method [51] that rely on tracking the divergence of independent trajectories, the Rosenstein algorithm [54] uses time-shifted neighbors within a single observed signal, making it particularly suitable for experimental or short, noisy datasets. Despite its efficiency and robustness, accurate estimation still depends on the careful selection of embedding parameters and sufficient data length.

Although these algorithms provide reliable results, they often require long, noise-free data sets and careful tuning of embedding parameters such as the delay time and embedding dimension. Moreover, the computational expense of LLE estimation increases rapidly with the number of data points, which can limit its applicability in extensive parameter sweeps or real-time analyses.

2.3.3. CD

The CD provides a quantitative measure of the geometric complexity or fractal dimensionality of an attractor reconstructed from a system’s time series [60,61]. It is derived from the broader concept of fractal dimension, which characterizes how a set fills the surrounding phase space as the scale of observation changes. In the context of nonlinear dynamics, the CD reflects how the number of pairs of points within a certain distance ε scales with distance, thus revealing the intrinsic structure of the attractor.

Consider a reconstructed trajectory consisting of N points $(x_{i})_{i = 1 . . . N}$ in an embedding space, where

x_{i} = (x_{i}, x_{i + τ}, \dots, x_{i + (m - 1) τ}),

(38)

herein $x_{i}$ is a point in the time-series plot; m is the embedding dimension, and τ is the time delay. Now, the correlation sum is defined as

C (ε) = \frac{2}{N (N - 1)} \sum_{i < j}^{N} H (ε - | x_{i} - x_{j} |),

(39)

where H is the Heaviside function, and $| x_{i} - x_{j} |$ is the Euclidean norm. Following [60,61] at sufficiently small ε, the following asymptotic relation takes place

C (ε) : ε^{d}, ε \to 0,

(40)

from where, the CD d can be defined as

d = lim_{ε \to 0} \frac{\ln C (ε)}{\ln ε} .

(41)

Despite its conceptual simplicity, the method is sensitive to several factors—most notably, data length, noise, and the choice of embedding parameters. In particular, both the embedding dimension m and the time delay τ, strongly affect the computed CD, making the algorithm potentially unstable, if wrong values of m and τ are chosen [62]. Moreover, computing pairwise distances between all points scales as $O (N^{2})$ , making the algorithm computationally demanding for large datasets.

2.3.4. SE

Originally introduced by Shannon [63] in the context of information theory, entropy quantifies the average amount of information produced by a stochastic or deterministic process. When applied to time-series analysis, SE measures the diversity of system states and the degree of disorder present in the underlying dynamics [64 –66].

Following [63], a system entropy can be defined as

H = - \sum_{i = 1}^{N} p_{i} \log_{2} p_{i},

(42)

where $p_{i}$ is the probability of a discrete state. One of the most attractive realizations of the probability evaluation relates to the introduction of the time-delay embedding (38) in either time-series or phase space, this m-dimensional space is then divided into N hypercubes, and the probability $p_{i}$ is defined as

p_{i} = k_{i} {(\sum_{i = 1}^{N} k_{i})}^{- 1},

(43)

where $k_{i}$ is the number of points within the corresponding hypercube.

One of the main advantages of SE lies in its simplicity and versatility. It can be applied directly to finite and noisy datasets without the need for explicit phase-space reconstruction or trajectory tracking, unlike geometric methods such as the Lyapunov exponent or CD. Similarly to the CD method, the SE method is heavily dependent on the choice of embedding parameters m and τ along with the number of hypercube divisions N, making the method potentially unstable, if wrong values for m and τ are adopted [65].

2.3.5. Some other fractal measures

Several other methods for quantification of fractal dimension include (1) the Minkowski–Bouligand method, known also as the box-counting method [66 –70]; (2) the information dimension method [71,72]; (3) the entropy-based evaluation using the Recurrence Quantification Analysis (RQA) [73,74]; and (4) the 0–1 test for chaos via the Gottwald–Melbourne method [75 –77], which transforms the original time-series into an auxiliary random-walk like process and evaluates whether the resulting trajectory remains bounded (indicating regular motion) or unbound (indicating chaos).

3. Numerical analysis

3.1. The model

Following [78 –82], consider LJ oscillator (13) with parameters

α \in [0.01; 25]; η \in [0.01; 0.1]; f \in [0.01; 25],

(44)

resembling those used in molecular dynamic simulations of cross-linked rubber-like materials [78,79,83,84]. The initial conditions are taken as

x (0) = 0; \partial_{τ} x (0) = 0 .

(45)

According to recommendations [33,40,48], computations were performed using the following parameters: transient periods—500; record periods—1000; relative tolerance—10⁻⁹; absolute tolerance—10⁻¹²; number of stages in the explicit part of the ABM solver—3; and number of stages in the implicit part—4.

3.2. Bifurcation diagrams

The bifurcation diagrams—defined as the vertical and horizontal projections of the Poincaré sections (36) versus driving force amplitude—are shown in Figure 1 for $α = 24$ and $η = 0.05$ .

These bifurcation diagrams clearly reveal the presence of a questionable forcing amplitude band below 5, within which chaotic or quasi-periodic behavior may emerge, whereas all other values in the studied range appear to be associated with periodic motion. A more detailed bifurcation diagrams near the questionable region are shown in Figure 2. Note that at smaller $α < ~ 19.2$ and any η, no questionable regions occur.

The bifurcation diagram in Figure 2(a) (displacement projection) reveals a rich sequence of dynamical transitions. Starting from a period-1 orbit at low forcing (∼4.3–4.4), then the system undergoes a transition region near ∼4.45, leading to chaos. Intermittent periodic windows (e.g., at ∼4.5 and near ∼4.9) are split by thin chaotic strips. A prominent crisis-induced intermittency is visible near different forcing amplitudes within this range. The dense black bands indicate fully developed chaos, while sparse regions correspond to regular motion. This intermittent structure is observed, apparently for the first time. Another observation is the apparent absence of Feigenbaum period-doubling and Feigenbaum-like cascades [85] near the onset of chaotic regimes.

The bifurcation diagram in Figure 2(b) (velocity projection) exhibits complementary behavior. Periodic windows are marked by sharp horizontal lines forming isolated clusters.

3.3. Fractal measures

Herein, the LLE via the Rosenstein algorithm, CD via the Grassberger–Procaccia algorithm, SE via the Ambika algorithm, and the 0–1 test for chaos via the Gottwald–Melbourne method are applied for the analysis of the considered LJ oscillator within a smaller questionable region, $f \in [4.3; 4.5]$ .

The following parameters were adopted for computing fractal measures: (a) LLE: embedding dimension $m = 4$ , time delay $τ = 7$ , following practical recommendations for small data sets [54,58]; (b) CD: embedding dimension $m = 5$ , time delay $τ = 7$ , as commonly adopted [60,61]; (c) SE: embedding dimension $m = 4$ , time delay $τ = 7$ , identical to the CD settings [64]; and (d) 0–1 test for chaos: 100 randomly chosen frequencies in the range $c \in (0; π)$ , following guidelines [75 –77].

The corresponding plots for displacement are shown in Figure 3.

The plots in Figure 3 reveal that all the considered methods give consistent results in chaos identification. These measures collectively capture the system’s transition from regular to chaotic dynamics.

The LLE plot (Figure 3(a)) reveals distinct peaks at specific forcing amplitudes, indicating regions where the system exhibits sensitivity to initial conditions—a defining feature of chaos. In contrast, intervals with near-zero LLE values suggest stable, periodic behavior. Within interval $f \in (4.38; 4.45)$ , the Lyapunov exponent, together with peaks well above zero—a hallmark of chaos—most likely indicates quasi-periodic behavior. The CD plot (Figure 3(b)) complements this by quantifying the fractal dimensionality of the attractor. Sharp increases in CD correspond to the emergence of complex, chaotic structures, while lower values reflect simpler, more regular dynamics. The SE plot (Figure 3(c)) measures the unpredictability of the system’s output. Elevated entropy values align with chaotic regimes, whereas dips in entropy mark transitions back to periodic or quasi-periodic motion. Finally, the results of the 0–1 test are presented in Figure 3(d) confirm the observations of predominantly regular dynamics over most of the examined range of forcing amplitudes, as indicated by consistently low values of the K-parameter. However, within a narrow interval of the control parameter $f \in (4.38; 4.45)$ , a transition region is observed with pronounced fluctuations, intermittently attaining values close to unity, which is characteristic of chaotic motion. Together, these plots demonstrate that the LJ oscillator undergoes multiple transitions between order and chaos within this narrow parameter window. The consistency across all three measures reinforces the reliability of these indicators in detecting subtle dynamical changes and highlights the rich complexity of the system’s behavior under varying external forcing.

It is worth noting that the sharp, multiple transitions from regular to chaotic motion observed here are atypical for most well-known nonlinear harmonic oscillators. In many classical systems, such as the Duffing or Van der Pol oscillators [48,49,86], the onset of chaos is typically characterized by a gradual progression through a Feigenbaum cascade—a sequence of period-doubling bifurcations that systematically lead to chaotic dynamics. In contrast, the LJ oscillator exhibits abrupt and closely spaced transitions, suggesting a more intricate and sensitive dependence on the driving force amplitude. This behavior highlights the unique complexity of the system and may point to underlying mechanisms beyond those captured by conventional bifurcation scenarios.

3.4. Phase portraits, Poincaré sections, and amplitude spectra

3.4.1. Quasi-periodic regimes

Consider the interval $f \in [4.3; 4.5]$ , analyzed in the previous section. At a lower ( $f = 4.3$ ) bound, the quasi-periodic regime appears (see Figure 4).

The motion regime in Figure 4 can be characterized as quasi-periodic, as evidenced by the qualitative features of both the phase-space (Figure 4(a)) and frequency-domain (Figure 4(b)) representations. In the Poincaré section, the trajectory forms a closed torus-like attractor, rather than collapsing onto a finite set of points or a strange attractor. Such a structure is a hallmark of quasi-periodic motion, in which two (or more) independent oscillatory modes coexist with incommensurate frequencies [48,60,61]. This interpretation is further supported by the amplitude spectrum (Figure 4(b)), which displays a discrete set of peaks distributed across a range of frequencies that are not integer multiples of a single fundamental frequency. The presence of these incommensurate frequency components indicates that the system response arises from the nonlinear interaction of multiple oscillatory modes rather than from harmonic distortions of a single periodic oscillation.

3.4.2. Chaotic regimes

According to plots in Figure 3, chaotic motion should occur at $f = 4.4$ . The following analysis confirms this (see Figure 5).

The system response at $f = 4.4$ corresponds to a chaotic regime, as indicated by the combined features of the Poincaré section (Figure 4(a)) and amplitude spectrum (Figure 4(b)). The Poincaré section reinforces the interpretation of chaotic behavior: instead of a discrete set of points (typical of periodic motion) or a closed continuous curve (indicative of quasi-periodic motion), the Poincaré section represents a typical strange attractor. This dispersion signifies the breakdown of the toroidal structure in phase space and the onset of deterministic chaos, where trajectories no longer lie on a smooth invariant manifold but instead evolve on a strange attractor. The amplitude spectrum (Figure 4(b)) exhibits both subharmonic and superharmonic continuous components of the amplitude spectrum, along with discrete singular peaks containing the driving frequency, too. This broadband nature reflects the irregular and aperiodic character of the signal, in contrast to the sharp spectral lines observed in periodic or quasi-periodic regimes.

3.5. LJ stiffness variation

The analysis of bifurcation diagrams constructed for varying stiffness α of LJ oscillator reveals that the emergence of entangled chaotic motion persists over a wide range $α \in [5; 25]$ of LJ stiffness values (Figure 6).

The bifurcation diagrams reveal that increasing stiffness causes an intermittent chaotic region to appear and migrate to higher values of the driving force (Figure 6). This behavior demonstrates that the emergence of the intermittent chaotic region remains stable under variation of LJ stiffness. Note also that similar results were obtained for the linear viscous coefficient variation $η \in [0.03; 0.1]$ .

4. Concluding remarks

The investigation of the LJ harmonic oscillator (equation (13)) demonstrates a complex dynamic behavior characterized by alternating regimes of quasi-periodicity and intermittent chaos. A key outcome of this study is the discovery that chaotic oscillations can arise even at low forcing amplitudes, while the system remains predominantly non-chaotic over a broad parameter range (Figure 2). The emergence of strange attractors with distinct, densely arranged patterns highlights the presence of deterministic chaos embedded within otherwise regular motion (Figure 5). This phenomenon appears to have been observed for the first time.

Comprehensive nonlinear analyses—including bifurcation diagrams, LLEs, CD, SE, and 0–1 test for chaos—have been employed to confirm and quantify the onset of chaos. The use of the variable-step ABM method in an explicit–implicit formulation ensured both numerical accuracy and stability in capturing subtle dynamical transitions.

Overall, the findings provide important insight into the chaotic nature of interatomic and intermolecular dynamics described by the LJ potential. The emergence of intermittent chaos in the forced LJ oscillator is particularly noteworthy, as it exemplifies a classical Pomeau–Manneville type intermittency [87] route to chaos in a realistic single-degree-of-freedom system with strongly anharmonic, distance-dependent restoring forces. The observed sensitivity to small perturbations and the formation of strange attractors suggest potential implications for understanding energy transfer mechanisms, nonlinear resonance phenomena, and complex vibrational dynamics in molecular and condensed matter systems.

In the context of atomic force microscopy (AFM), where the cantilever tip interacts with a sample surface via forces well approximated by an LJ potential, the observed intermittent chaos provides a plausible dynamical origin for several experimentally reported phenomena, including irregular tip-sample contact events, noisy cantilever response near resonance, and apparent broadband power spectra in certain parameter regimes. Recognizing this chaotic intermittency is therefore essential for correct interpretation of force spectroscopy data, for distinguishing true sample properties from probe artifacts, and for improving the reliability of nanoscale imaging and manipulation on soft matter and biological surfaces.

Another important practical application concerns much larger-scale engineering devices, such as high-damping rubber bearings (HDRBs) and laminated rubber isolators widely employed in seismic base isolation and vibration mitigation systems for buildings, bridges, and critical infrastructure. These bearings, often fabricated from hyperelastic rubber compounds with embedded steel plates, exhibit strongly nonlinear force-deformation characteristics including strain-stiffening at large shear deformations that can be effectively approximated using hyperelastic potentials of the LJ type [88 –90]. The observation of intermittent chaos in the forced LJ oscillator, featuring long near-periodic phases punctuated by abrupt chaotic bursts, may carry direct implications for understanding irregular dynamic responses in real rubber-based isolators under extreme loading.

Footnotes

ORCID iD

Sergey V. Kuznetsov

Funding

The author received no financial support for the research, authorship, and/or publication of this article.

Declaration of conflicting interests

The author declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Research involving human participants and/or animals

This article does not contain any studies with human participants and/or animals.

Informed consent

The author confirms knowledge of what this study involves.

References

Jones

. On the determination of molecular fields.—I. From the variation of the viscosity of a gas with temperature. Proc R Soc Lond A 1924; 106: 441–462.

Albe

Nord

Nordlund

. Dynamic charge-transfer bond-order potential for gallium nitride. Philos Mag 2009; 89(34–36): 3477–3497.

Stone

. The theory of intermolecular forces. Oxford University Press, 2016.

Mie

. Zur kinetischen Theorie der einatomigen Körper. Annalen Der Physik 1903; 316(8): 657–697.

Fischer

Wendland

. On the history of key empirical intermolecular potentials. Fluid Ph Equilib 2023; 573: 113876.

Lennard-Jones

. Cohesion. Proc Phys Soc 1931; 43(5): 461–482.

Stephan

Thol

Vrabec

, et al. Thermophysical properties of the Lennard-Jones fluid: database and data assessment. J Chem Inf Model 2019; 59: 4248.

Le, Roy

Van, Kranendonk

. Anisotropic intermolecular potentials from an analysis of spectra of H2- and D2-inert gas complexes. J Chem Phys 1974; 61(11): 4750–4769.

Epstein

Showalter

. Nonlinear chemical dynamics: oscillations, patterns, and chaos. J Phys Chem 1996; 100(31): 13132–13147.

10.

Herman

Ogilvie

. An effective Hamiltonian to treat adiabatic and nonadiabatic effects in the rotational and vibrational spectra of diatomic molecules. Adv Chem Phys 1998; 103: 187–215.

11.

Hajigeorgiou

Le Roy

. A “modified Lennard-Jones oscillator” model for diatom potential functions. J Chem Phys 2000; 112(9): 3949–3957.

12.

Tsyganov

. The shock forced oscillator model: a harmonic approximation. High Temp 2014; 52: 48–56.

13.

Vitória

RLL

Belich

. Harmonic oscillator in an environment with a point-like defect. Phys Scr 2019; 94: 125301.

14.

Vitória

RLL

da, Silva

KAT

. Effects of small oscillations on a diatomic molecule in an elastic medium. Eur Phys J Plus 2023; 138: 172.

15.

Bezerra

JES

Vitória

RLL

Moucherek

FMO

. Lennard-Jones oscillations in an elastic environment. Appl Math 2025; 5(4): 129.

16.

Marx

Hutter

. Ab initio molecular dynamics: basic theory and advanced methods. Cambridge University Press, 2009.

17.

Child

. Semiclassical mechanics with molecular applications. Clarendon, 1991.

18.

Ueda

. Randomly transitional phenomena in the system governed by Duffing’s equation. J Stat Phys 1979; 20(2): 181–196.

19.

Ueda

. Explosion of strange attractors exhibited by Duffing’s equation. Annals NY Acad Sci 1980; 357: 422–434.

20.

Ueda

. Random phenomena resulting from non-linearity in the system described by Duffing’s equation. Int J Non-Linear Mech 1985; 20: 481–491.

21.

Ueda

. Survey of regular and chaotic phenomena in the forced Duffing oscillator. Chaos Solitons Fractals 1991; 1(3): 199–231.

22.

Negron-Marrero

Lopez-Serrano

. Minimal energy configurations of finite molecular arrays. Symmetry 2019; 11: 158.

23.

Yao

. Dressed active particles in spherical crystals. Soft Matter 2016; 12: 7020–7027.

24.

Yao

. Collective dynamics and shattering of disturbed two-dimensional Lennard–Jones crystals. Eur Phys J E 2022; 45: 88.

25.

Acharjee

Dutta

Devi

, et al. Nonlinear oscillations, chaotic dynamics, and stability analysis of bilayer graphene-like structures. Chaos 2023; 33(1): 013136.

26.

Baek

Lee

M-H

, et al. Characterization of a forced damped harmonic oscillator in a magneto-optical trap of 133Cs atoms. J Opt Soc Am 2023; B40: 1735–1741.

27.

Calvo

. Chaos and dynamical coexistence in Lennard-Jones clusters. J Chem Phys 1998; 108(16): 6861–6867.

28.

Sagues

Epstein

. Nonlinear chemical dynamics. Dalton Trans 2023; 2003: 1201–1217.

29.

Ciesla

Pawlowicz

Longa

. Molecular dynamics simulation of the Lennard-Jones polymers in a good solvent. Acta Phys Polonica B 2007; 38(5): 1727.

30.

Sehgal

Maroudas

Ford

. Phase behavior of the 38-atom Lennard-Jones cluster. J Chem Phys 2014; 140(10): 104312.

31.

Showalter

Epstein

. From chemical systems to systems chemistry: patterns in space and time. Chaos 2015; 25(9): 097613.

32.

Epstein

Pojman

. An introduction to nonlinear chemical dynamics: oscillations, waves, patterns, and chaos. Oxford Academic, 2020.

33.

Hairer

Nørsett

Wanner

. Solving ordinary differential equations I: nonstiff problems. 2nd ed. Springer-Verlag, 1993.

34.

Zhao

Zhang

. Comparison of different order Adams-Bashforth methods in an atmospheric general circulation model. Acta Meteorol Sin 2011; 25: 754–764.

35.

Hartman

. Ordinary differential equations. Wiley, 1964.

36.

Xiao

Belytschko

. A bridging domain method for coupling continua with molecular dynamics. Comp Meth Appl Mech Eng 2004; 193(17–20): 1645–1669.

37.

Kelly

Macmillan

. Strong solids. Clarendon, 1986.

38.

Dang

Wang

, et al. Extreme mechanics of nanoscale diamond towards functional device applications. Extreme Mech Lett 2023; 58: 101931.

39.

Lambert

. Numerical methods for ordinary differential systems: the initial value problem. Wiley, 1973.

40.

Hairer

Wanner

. Solving ordinary differential equations II: stiff and differential-algebraic problems. 2nd ed. Springer-Verlag, 1996.

41.

Misirli

Gurefe

. Multiplicative Adams–Bashforth–Moulton methods. Numer Algor 2011; 57: 425–439.

42.

Brun

Djeran-Maigre

, et al. Explicit/implicit multi-time step co-simulation in unbounded medium with Rayleigh damping and application for wave barrier. Eur J Environ Civ Eng 2020; 24: 2400–2421.

43.

Brun

Djeran-Maigre

, et al. Benchmark for three-dimensional explicit asynchronous absorbing layers for ground wave propagation and wave barriers. Comp Geotech 2021; 131: 103808.

44.

Terentjeva

. Wave fields and domination regions for the interior Lamb problem. Mech Solids 2015; 50(5): 508–520.

45.

Kuznetsov

. Love waves in layered anisotropic media. J Appl Math Mech 2006; 70(1): 116–127.

46.

Kuznetsov

. SH-waves in laminated plates. Quart Appl Math 2006; 64(1): 153–165.

47.

Kuznetsov

. Appearing ZGV point in the first flexural branch of Lamb waves in multilayered plates. Composite Struct 2022; 290: 115532.

48.

Kovacic

Brennan

. The Duffing equation: nonlinear oscillators and their behavior. John Wiley and Sons, 2011.

49.

Guckenheimer

Holmes

. Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. Springer-Verlag, 1983.

50.

Grebogi

Ott

Yorke

. Crises, sudden changes in chaotic attractors, and transient chaos. Physica D Nonlinear Phenom 1983; 7(1–3): 181–200.

51.

Benettin

Galgani

Giorgilli

, et al. Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems: a method for computing all of them, part 1: theory. Meccanica 1979; 15(1): 9–20.

52.

Shimada

Nagashima

. A numerical approach to ergodic problems of dissipative dynamical systems. Progress Theor Phys 1979; 61(6): 1605–1616.

53.

Wolf

Swift

Swinney

, et al. Determining Lyapunov exponents from a time series. Physica D Nonlin Phenom 1985; 16(3): 285–317.

54.

Rosenstein

Collins

De, Luca

. A practical method for calculating largest Lyapunov exponents from small data sets. Physica D Nonlin Phenom 1993; 65(1–2): 117–134.

55.

Bruijn

Meijer

Beek

, et al. Maximum Lyapunov exponents as predictors of global gait stability: a modelling approach. Med Eng Phys 2012; 34(4): 428–436.

56.

Beaudette

Howarth

Graham

, et al. On the use of a Euclidean norm function for the estimation of local dynamic stability from 3D kinematics using time-delayed Lyapunov analyses. Med Eng Phys 2016; 38(10): 1139–1145.

57.

Mehdizadeh

. The largest Lyapunov exponent of gait in young and elderly individuals: a systematic review. Gait & Posture 2018; 60: 241–250.

58.

Mehdizadeh

. A robust method to estimate the largest Lyapunov exponent of noisy signals: a revision to the Rosenstein’s algorithm. J Biomech 2019; 85: 84–91.

59.

Liu

H-F

Dai

Z-H

W-F

, et al. Noise robust estimates of the largest Lyapunov exponent. Phys Lett A 2005; 341: 119–127.

60.

Grassberger

Procaccia

. Measuring the strangeness of strange attractors. Physica D Nonlin Phenom 1983; 9(1–2): 189–208.

61.

Grassberger

Procaccia

. Characterization of strange attractors. Phys Rev Lett 1983; 50(5): 346–349.

62.

Kennel

Brown

Abarbanel

HDI

. Determining embedding dimension for phase-space reconstruction using a geometrical construction. Phys Rev A 1992; 45(6): 3403–3411.

63.

Shannon

. A mathematical theory of communication. Bell System Tech J 1948; 27(3): 379–423.

64.

Webber

Zbilut

. Dynamical assessment of physiological systems and states using recurrence plot strategies. J Appl Physio 1994; 76(2): 965–973.

65.

Misra

Harikrishnan

Ambika

, et al. The correlation entropy of X-ray light curves from GRS 1915+105. Astrophys J 2004; 609(1): 313–320.

66.

Il’yasov

Kravtsov

Kuznetsov

, et al. Exterior 3D lamb problem: harmonic load distributed over a surface. Mech Solids 2016; 51: 39–45.

67.

Bouligand

. Les Definitions Modernes de la Dimension. Hermann, 1935.

68.

Mandelbrot

. The fractal geometry of nature. Freeman and Company, 1982.

69.

Falconer

. Fractal geometry: mathematical foundations and applications. 3rd ed. Wiley, 2014.

70.

Ghasemi

Yun

. Fractal structures arising from interfacial instabilities in bio-oil atomization. Sci Rep 2021; 11: 411.

71.

Rényi

. On the dimension and entropy of probability distributions. Acta Math Acad Sci Hungaricae 1959; 10(1–2): 193–215.

72.

Verdu

. Rényi information dimension: fundamental limits of almost lossless analog compression. IEEE Trans Inform Theory 2010; 56(8): 3721–3748.

73.

Eckmann

Kamphorst

Ruelle

. Recurrence plots of dynamical systems. Europhys Lett 1987; 4(9): 973–977.

74.

Marwan

Romano

Thiel

, et al. Recurrence plots for the analysis of complex systems. Phys Rep 2007; 438(5–6): 237–329.

75.

Gottwald

Melbourne

. A new test for chaos in deterministic systems. Proc Roy Soc A 2004; 460: 603–611.

76.

Gottwald

Melbourne

. Testing for chaos in deterministic systems with noise. Phys D 2005; 212: 100–110.

77.

Gottwald

Melbourne

. On the implementation of the 0–1 test for chaos. SIAM J Appl Dyn Syst 2009; 8: 129–145.

78.

Davidson

Goulbourne

. A nonaffine network model for elastomers undergoing finite deformations. J Mech Phys Solids 2013; 61(8): 1784–1797.

79.

Davidson

Goulbourne

. Nonaffine chain and primitive path deformation in crosslinked polymers. Modelling Simul Mater Sci Eng 2016; 24: 065002.

80.

Dudchenko

. Vertical wave barriers for vibration reduction. Arch App Mech 2021; 91: 257–276.

81.

Kumar

Brassart

. On tube models of rubber elasticity: fitting performance in relation to sensitivity to the invariant I2. Mech Soft Mater 2023; 5(1): 54.

82.

Araujo

Kryven

Brassart

. Micromechanical modelling of rubbery networks: the role of chain pre-stretch. Int J Non-linear Mech 2024; 166: 104834.

83.

Van, der Hoef

. Free energy of the Lennard-Jones solid. J Chem Phys 2000; 113: 8142–8148.

84.

Toda

Morita

. Rubber elasticity of realizable ideal networks. AIP Advances 2018; 8(12): 125005.

85.

Feigenbaum

. Quantitative universality for a class of nonlinear transformations. J Statist Phys 1978; 19: 25–52.

86.

Kuznetsov

. Forced harmonic vibrations of duffing oscillator with cubic-quintic spring nonlinearity. Chaos Solitons Fractals 2025; 199: 116699.

87.

Pomeau

Manneville

. Intermittent transition to turbulence in dissipative dynamical systems. Commun Math Phys 1980; 74: 189–197.

88.

Bratov

Ilyashenko

Kuznetsov

, et al. Homogeneous horizontal and vertical seismic barriers: mathematical foundations and dimensional analysis. Mater Phys Mech 2020; 44(1): 61–65.

89.

Bratov

Kuznetsov

Morozov

. Seismic barriers filled with solid elastic and granular materials: comparative analysis. Math Mech Solids 2022; 27(9): 1761–1770.

90.

Oyerokun

Schweizer

. Microscopic theory of rubber elasticity. J Chem Phys 2004; 120(19): 9359–9370.