Non-Integrability of a Coupled Two-Mode Quartic Hamiltonian via Kovalevskaya Exponents

The reliable analysis of complex, multi-physics engineering systems frequently begins by reducing a high-dimensional model to a small number of interacting degrees of freedom. Coupled microbeam resonators, post-buckled composite panels, electromechanical transducers, and fluid--structure assemblies are all commonly represented, near an operating point, by two or three nonlinearly coupled modal oscillators that trade energy through geometric or field-mediated coupling. The predictive value of such a reduced-order model rests on whether the underlying dynamics admit additional conserved quantities: a second independent integral would allow the modal interaction to be decoupled or linearized exactly, whereas its absence signals that energy transfer between modes is an intrinsic feature that no change of coordinates can remove. Establishing which situation prevails, for a given choice of physical parameters, is therefore a question of direct engineering consequence, bearing on energy localization, complex resonant response, and the admissibility of simplified design models. The present work studies a representative two-mode Hamiltonian for which this question can be settled rigorously and analytically.

The dynamics of the early universe can be captured, within the framework of quantum cosmology, by finite-dimensional Hamiltonian systems obtained from a suitable mini-superspace reduction of the Einstein--Hilbert action. Among these, the Friedmann--Robertson--Walker (FRW) model coupled to a scalar field occupies a central position, both because of its physical richness and because the resulting system exhibits a wealth of dynamical phenomena, including chaos, resonance, and multi-family periodic orbits. The same system, as detailed below, provides the analytical proxy referred to above: its potential coincides in form with that of a symmetric two-mode nonlinear structure. The spacetime metric of the FRW model is given by:

where, $a(\eta)$ is the scale factor, $\eta$ is conformal time, $d\Omega_{2}^{2}$ is the round metric on the unit two-sphere, and $K\in\{0,\pm1\}$ is the spatial curvature index. The general action governing the coupled gravity—scalar-field system reads:

$I=\frac{c^{4}}{16\pi G}\int\left[R-2\Lambda-\frac{1}{2}\left(\nabla_{\mu}\varphi\nabla^{\mu}\varphi+V(\varphi)\right)+\xi R|\varphi|^{2}-\sigma\right]\sqrt{-g}\,d^{4}x$

where, $R$ is the Ricci scalar, $\Lambda$ is the cosmological constant, $V(\varphi)$ is the scalar potential, $\xi$ is the non-minimal coupling constant, and $\sigma$ is the energy density of a perfect fluid background. The potential includes a quadratic mass term $m^{2}|\varphi|^{2}$, where $m$ denotes the field mass.

For the conformal-coupling value $\xi=1/6$ and the quartic fluid density $\sigma=\lambda|\varphi|^{4}/24$, the mini-superspace reduction yields the two-degree-of-freedom potential:

where, $q_{1}$ and $q_{2}$ are real canonical coordinates, and $\lambda,\Lambda,m^{2}\in\mathbf{R}$, $K\in\{0,\pm1\}$.

Beyond its cosmological provenance, the potential defined in Eq.(2) coincides in structure with the potential energy of a symmetric two-mode nonlinear mechanical system, and it is this correspondence that this study exploits throughout. Read in engineering terms, $q_{1}$ and $q_{2}$ are modal coordinates measuring the amplitudes of two interacting vibration modes. The quadratic group $k(-q_{1}^{2}+q_{2}^{2})$ assigns the first mode an effective negative linear stiffness and the second a positive one, so that $k$ acts as an axial-load or pre-stress parameter that drives the first mode through a buckling threshold as it changes sign. The quartic coefficients $\Lambda$ and $\lambda$ are the cubic, Duffing-type stiffnesses of the two modes, the signature of geometric nonlinearity arising from large-amplitude deflection and mid-plane stretching in beams, plates, and panels; they harden the response when positive and soften it when negative. The cross term $m^{2}q_{1}^{2}q_{2}^{2}/2$ is the lowest-order amplitude--amplitude coupling between the modes, and its coefficient $m^{2}$ governs the strength of nonlinear modal interaction and hence the capacity of the system to transfer energy between modes, the mechanism underlying internal resonance and energy localization. The invariance of Eq. (2) under independent sign reversals of $q_{1}$ and $q_{2}$ mirrors the geometric symmetry of such structures. Figure 1 depicts a representative lumped realization together with the resulting configuration space. This correspondence motivates treating the system described by Eq. (2) as a proxy through which the integrability, stability, and energy-transfer properties of an entire family of engineering systems may be examined by exact analysis.

The potential defined in Eq.(2) and its special cases have appeared repeatedly in the cosmological dynamics literature. Calzetta and El Hasi [1] investigated the following special form and demonstrated the existence of chaotic behavior:

$V=\frac{1}{2}(y^{2}-x^{2})+\frac{b}{2}x^{2}y^{2}$

Subsequently, Llibre and Makhlouf [2] generalized the Calzetta--El Hasi model to the form as follows:

$V=\frac{1}{2}(y^{2}-x^{2})+\frac{b}{2}x^{2}y^{2}+\frac{a}{4}x^{4}+\frac{c}{4}y^{4}$

This coincides with Eq.(2) (up to a change of sign in the quadratic terms) and may be interpreted as a model of galactic dynamics. The $\mathbf{Z}_2 \times \mathbf{Z}_2$ symmetry of this potential (invariance under independent sign reversals of $x$ and $y$) was exploited by Pucacco et al. [3] to obtain quantitative predictions for the periodic orbit structure via detuned normal form theory. The identical quartic structure, comprising two grounded cubic stiffnesses and a single quartic cross-coupling, recurs in the engineering analysis of nonlinearly coupled oscillators, post-buckled plates, and two-mode microelectromechanical systems resonators, which is the basis for the proxy interpretation adopted in the study. The same model arises in the study of the cosmological and thermodynamic arrows of time within the quantum cosmology framework [4], [5]. The Hamiltonian system associated with the potential defined in Eq.(2) possesses two degrees of freedom and is integrable in the Liouville—Arnold sense [6] if and only if there exist two functionally independent, single-valued, analytic, mutually involutive first integrals. Since no general algorithm for establishing Liouville integrability is available, one typically resorts either to constructing the second integral explicitly or to proving its non-existence.

Several powerful criteria for non-integrability have been developed. The Painlevé analysis [7-8] seeks necessary conditions for homomorphic solutions and has been applied successfully to Hamiltonian systems of this type. Ziglin’s method [9] certifies non-integrability by examining the monodromy group of the variational equations along a known particular solution; Shibayama [10] refined this into an explicit non-existence criterion for systems with homogeneous potentials of integer degree. The differential Galois approach of Morales-Ruiz et al. [11] provides the sharpest available criterion, establishing non-integrability through the non-virtual-abelianness of the Galois group of the higher variational equations.

Among the available techniques, the Kovalevskaya exponent method occupies a distinguished position because it simultaneously furnishes necessary conditions for integrability and provides a direct algebraic route to the singularity structure of solutions. The method originated in the classical work of Kovalevskaya on the integrability of the rigid body and was systematized for modern Hamiltonian mechanics by Yoshida [12], who established that a necessary condition for the existence of an algebraic first integral is that all Kovalevskaya exponents be rational. If at least one exponent is irrational or complex, the system cannot possess an additional analytic first integral and is therefore non-integrable.

Building on this foundation, Yoshida [13] demonstrated that $\rho=-1$ is always a Kovalevskaya exponent and that non-rationality of a single remaining exponent suffices to preclude integrability. The conditions relating non-rationality of Kovalevskaya exponents to the non-existence of a second integral were further refined by Furta [14] and Goriely [15], who extended the theory to systems of ordinary differential equations beyond the Hamiltonian setting. Llibre [16] has investigated the connection between Kovalevskaya exponents and polynomial first integrals for the Lotka--Volterra system, and a comprehensive historical survey of these developments is given by Goriely [17]. For Hamiltonian systems in particular, Szumiński [18] established a pairing theorem for Kovalevskaya exponents: if the Hamiltonian is similarity invariant of weighted degree $h$, then the Kovalevskaya exponents come in pairs summing to $(h-1)$. The relation between Kovalevskaya exponents and Painlevé resonances was clarified by Yoshida et al. [19], who proved that, provided all particular solutions of the similarity-invariant system are non-trivial, the two sets of exponents coincide.

From the standpoint of engineering modeling frameworks, the appeal of the Kovalevskaya exponent method is that it operates as an inexpensive analytical pre-screen for the admissibility of reduced-order models. A second independent analytic integral is precisely the object whose existence would permit the two-mode interaction to be decoupled into independent normal coordinates; its certified absence therefore rules out any global linearizing or decoupling transformation, with direct consequences for design and control. This diagnostic reading is developed in Section~2 and revisited in the conclusions. The present work pursues two complementary objectives. First, the Kovalevskaya exponent method is applied to prove the non-integrability of the FRW scalar-field Hamiltonian in Eq.~\eqref{eq3} in generic parameter regimes, and the exceptional cases in which all Kovalevskaya exponents are rational are identified. Second, Lyapunov’s theorem is applied to holomorphic integrals [20], as developed in the context of rigid-body mechanics by Abdel-Aziz [21] and El-Sabaa et al. [22], to construct families of periodic solutions near the equilibrium points of the system and to classify those equilibria according to their dynamical type. An application of the rigid body problem with external and gyroscopic torques is studied [23] for engineering, aerospace, and mathematical applications.

The main contributions of this work are as follows. First, a two-mode quartic Hamiltonian drawn from FRW cosmology is recast as a proxy for coupled nonlinear engineering systems, and an explicit dictionary between its parameters and structural quantities is provided (Figure 1). Second, non-integrability is established in the generic parameter regime through the Kovalevskaya exponent method, and the three exceptional regimes in which all exponents are rational are delimited. Third, all nine equilibrium configurations are classified as stable operating points or buckling-type instabilities, and the number of periodic vibration families about each is related to the parameters. Fourth, the methodological implication for modeling frameworks is stated: the lack of a second integral confirms that decoupled reduced-order models are inadmissible in the corresponding regimes. The remainder of the study is organized as follows. Section 2 derives the Kovalevskaya matrix, computes the three classes of exponents, and establishes the non-integrability result together with the exceptional parameter regimes with rational Kovalevskaya exponents. Section 3 constructs the periodic solutions via Lyapunov’s method, classifies all nine equilibrium points, and discusses the number of periodic families associated with each stability type.

The Hamiltonian associated with the potential defined in Eq.~\eqref{eq2} is:

where, $\lambda,\Lambda,m^{2}\in\mathbf{R}$ and $k\in\{0,\pm1\}$. Note that the kinetic energy carries an indefinite signature, reflecting the pseudo-Riemannian character of the FRW mini-superspace metric. Hamilton’s equations $\dot{q}_{i}=\partial H/\partial p_{i}$ and $\dot{p}_{i}=-\partial H/\partial q_{i}$ yield:

To extract the dominant singularity structure, the quasi-homogeneous scaling is applied:

where the weight exponent $\bar{g}=2/(n-2)$ is determined by the degree $n$ of the leading potential. Since the highest-degree monomials in $V$ are of degree $n=4$, this study obtains $\bar{g}=1$. In the limit $\alpha\to\infty$, the lower-degree terms (proportional to $k$) are suppressed, and the system in Eq. (4) reduces to the similarity-invariant system in Eq. (6) with the corresponding homogeneous Hamiltonian in Eq. (7):

By virtue of its quasi-homogeneity, the system in Eq. (6) admits power-law particular solutions of the form:

where the weight exponents $(\bar{g}_{1},\bar{g}_{2},\bar{g}_{3},\bar{g}_{4})=(1,1,2,2)$ are consistent with the scaling in Eq. (5). Substituting Eq. (8) into Eq. (6) yields the algebraic balance conditions:

where, $\dot{x}_{i}=F_{i}(\boldsymbol{x})$ with $\boldsymbol{x}=(q_{1},q_{2},p_{1},p_{2})$.

To determine the Kovalevskaya exponents, the particular solution is perturbed by setting $x_{i}=(c_{i}+\eta_{i})t^{-\bar{g}_{i}}$ and linearizing. The resulting variational equation is:

The ansatz $\boldsymbol{\eta}\propto t^{\rho}$ transforms Eq. (10) into the eigenvalue problem $K\boldsymbol{\eta}=\rho\boldsymbol{\eta}$, where the Kovalevskaya matrix is:

A direct computation of the Jacobian of Eq. (6) evaluated at the particular solution in Eq. (8), followed by addition of the diagonal correction $\operatorname{diag}(\bar{g}_{1},\bar{g}_{2},\bar{g}_{3},\bar{g}_{4})=\operatorname{diag}(1,1,2,2)$, gives:

The eigenvalues of $K$, referred to as the Kovalevskaya exponents, encode the integrability properties of the system in Eq. (6). Two exponents are universal: Yoshida [12], [13] proved that $\rho=-1$ is always a Kovalevskaya exponent (arising from time-translation symmetry), and that $\rho=2\bar{g}+2=4$ is a Kovalevskaya exponent whenever $\operatorname{grad}H_{\bar{g}}(\boldsymbol{c})=0$ (arising from the weighted scaling symmetry of $H_{\bar{g}}$). The two remaining exponents depend on the particular solution class. Setting $c_{1}=0$ or $c_{2}=0$ in the balance conditions in Eq. (9) yields two degenerate classes; the generic class has both $c_{1}\neq0$ and $c_{2}\neq0$. These three classes produce the following sets of Kovalevskaya exponents:

where, $\Omega_{1}$ and $\Omega_{2}$ correspond to the branches $c_{1}=0$ and $c_{2}=0$, respectively, in which only $\lambda$ and $m^{2}$ as well as only $\Lambda$ and $m^{2}$ appear, respectively; $\Omega_{3}$ is the generic branch, in which all three parameters are present. Yoshida’s criterion [12] asserts that a necessary condition for the existence of an algebraic first integral of the system in Eq. (6) is that all Kovalevskaya exponents be rational numbers. Inspection of the discriminants in Eq. (13) reveals that the Kovalevskaya exponents of $\Omega_{1}$ are rational if and only if $\lambda^{2}-8m^{2}\lambda$ is a perfect square; analogously for $\Omega_{2}$ and $\Omega_{3}$. For generic values of the parameters $m^{2}$, $\Lambda$, and $\lambda$, at least one Kovalevskaya exponent is irrational or complex, and the system in Eq. (6) is therefore non-integrable.

The list of exceptional parameter regimes in which all three classes $\Omega_{1}$, $\Omega_{2}$, and $\Omega_{3}$ simultaneously yield rational exponents is:

These cases, in which the Kovalevskaya exponent method does not preclude the existence of a second integral, warrant independent investigation; however, a full integrability analysis for each of them lies beyond the scope of the present work.

The Kovalevskaya exponent test does more than classify singularities; it operates as a rigorous diagnostic for the admissibility of reduced-order models of the analogous engineering system. A second independent analytic integral is precisely what would allow the two-mode interaction to be decoupled, or globally linearized, into independent normal coordinates. Its demonstrable absence in the generic regime, therefore, certifies that no smooth change of variables reduces the system to a pair of uncoupled oscillators so that any reduced-order model tacitly assuming such an invariant can misrepresent the global response. The three exceptional regimes in Eq. (14) delimit the only parameter combinations for which this obstruction is lifted and a near-integrable, decoupled description may be sought. For design and control, the generic regime marks the parameter sets in which modal energy transfer cannot be removed by coordinate choice and must instead be confronted directly.

As for the full system governed by the potential in Eq. (2), the equations of motion read in terms of the configuration coordinates:

where, $V_{i}=\partial V/\partial q_{i}$, and the system admits the energy first integral:

The equilibrium points $(q_{10},q_{20})$ can be obtained by solving the equations $V_{1}=V_{2}=0$, i.e.,

The system in Eq. (17) admits at most nine solutions:

(i) $Q_{1}=(0,0)$, the origin, which is always an equilibrium.

(ii) $Q_{2,3}=(0,\pm\sqrt{-k/\lambda})$, which exists when $k/\lambda<0$.

(iii) $Q_{4,5}=(\pm\sqrt{k/\Lambda},0)$, which exists when $k/\Lambda>0$.

(iv) $Q_{6,7,8,9}=\left(\pm\sqrt{\dfrac{k(\lambda+m^{2})}{\Lambda\lambda-m^{4}}},\pm\sqrt{\dfrac{-k(\Lambda+m^{2})}{\Lambda\lambda-m^{4}}}\right)$, which exists when $k(\lambda+m^{2})(\Lambda\lambda-m^{4})>0$ and $k(\Lambda+m^{2})(\Lambda\lambda-m^{4})<0$.

To determine the dynamical type of each equilibrium, this study introduces small perturbations $q_{i}=q_{i0}+\xi_{i}$ and linearizes Eq. (15). After retaining only first-order terms and writing $\alpha=V_{11}(q_{10},q_{20})$, $\beta=V_{22}(q_{10},q_{20})$, and $\gamma=V_{12}(q_{10},q_{20})$, the linearized equations take the form:

where,

$\begin{aligned} & V_{11}=-k+m^2 q_{20}^2+3 \Lambda q_{10}^2 \equiv \alpha, \\ & V_{22}=k+m^2 q_{10}^2+3 \lambda q_{20}^2 \equiv \beta, \\ & V_{12}=2 m^2 q_{10} q_{20} \equiv \gamma\end{aligned}$

The stability character of each equilibrium is determined by the sign of the Hessian determinant $\alpha\beta-\gamma^{2}$ and the sum $\alpha+\beta$:

i) If $\alpha\beta-\gamma^{2}<0$, the equilibrium is a saddle point and the frequency in Eq. (23) has exactly one real positive root, yielding a single family of periodic orbits.

ii) If $\alpha\beta-\gamma^{2}>0$ and $\alpha+\beta>2\sqrt{\alpha\beta-\gamma^{2}}$, the frequency equation has two distinct real positive roots, yielding two families of periodic orbits about a stable (minimum) equilibrium.

The complete classification of all nine equilibria is summarized in Table 1. Interpreted in the engineering setting of Figure 1, the equilibrium classification of Table 1 is a map of the system’s operating points. A potential minimum corresponds to a statically stable configuration, a stable operating point about which the structure executes bounded oscillations; the two distinct real frequencies furnished by Eq. (23) are then the two vibration modes accessible in its neighborhood. A saddle corresponds to a buckling-type instability: the configuration is statically unstable, an escape direction exists, and only a single bounded oscillation family survives. The parameters therefore act as design levers. The sign and magnitude of the load parameter $k$ relative to the cubic stiffnesses $\Lambda$ and $\lambda$ determine which operating points exist and whether the off-axis configurations $Q_{6-9}$ are realized, whereas the coupling coefficient $m^{2}$ enters the off-diagonal Hessian term $\gamma=2m^{2}q_{10}q_{20}$ and thus sets the splitting of the two vibration frequencies and the degree of modal interaction at a stable operating point. Parameter choices that place the system at a minimum with well-separated frequencies favor predictable, weakly coupled response, whereas choices that drive the system toward a saddle signal the onset of buckling and the loss of one vibration family.

Equilibrium PointMaximumMinimumSaddle $Q_{1}$—$k=\pm1$$Q_{2,3}$$k=1$, $\lambda<0$, $m^{2}>0$, $1>-m^{2}/\lambda$$k=-1$, $\lambda>0$, $m^{2}>0$$k=1$, $\lambda<0$, $m^{2}>0$, $1<-m^{2}/\lambda$$Q_{4,5}$$k=-1$, $\Lambda<0$, $m^{2}>0$, $1>-m^{2}/\Lambda$$k=1$, $\Lambda>0$, $m^{2}>0$$k=-1$, $\Lambda<0$, $m^{2}>0$, $1<-m^{2}/\Lambda$$Q_{6,7,8,9}$$k=1$, $\Lambda<0$, $(m^{2}+\Lambda)<0$, $(m^{2}+\lambda)>0$, $(m^{4}-\Lambda\lambda)<0$$k=-1$, $\Lambda>0$, $(m^{2}+\Lambda)>0$, $(m^{2}+\lambda)<0$, $(m^{4}-\Lambda\lambda)<0$$k=\pm1$, $k^{2}(m^{2}+\Lambda)(m^{2}+\lambda)(m^{4}-\Lambda\lambda)<0$

Families of periodic solutions near each equilibrium are constructed by applying Lyapunov’s method [21]. Under the time rescaling $u=\nu t$, the following can be derived from Eq. (18):

This study seeks solutions in the form of asymptotic series in a small amplitude parameter $\varepsilon$:

where the frequency admits the expansion:

At first order in $\varepsilon$, after substituting Eq. (20) into Eq. (19) and representing the unknowns $x^{(i)}$ and $y^{(i)}$ as Fourier series in $u$ (with terms up to the fifth harmonic), matching trigonometric coefficients yields the leading-order periodic solution:

where, $c$ and $f$ are free parameters. The leading-order frequency is determined by the characteristic equation of the system in Eq. (19):

The double sign in Eq. (23) reflects the two possible oscillation frequencies of the coupled system: both are real and distinct precisely when the equilibrium is stable (minimum), giving two independent families of periodic orbits parameterized by $(c,f)$; when the equilibrium is a saddle, only the $+$ sign yields a real frequency, giving one such family. In the engineering reading, these families are the small-amplitude vibration modes about the operating point, and their number is the count of independent oscillatory responses the structure can sustain locally.

This study established the non-integrability of the two-degree-of-freedom Hamiltonian system governing a conformally coupled scalar field in the FRW universe by demonstrating that the Kovalevskaya exponents associated with each class of particular solutions of the similarity-invariant system are generically complex. Three exceptional parameter regimes in Eq. (14), in which all exponents are rational and the method is inconclusive, were identified; the integrability status of these cases remains an open question. In addition, families of periodic solutions near all nine equilibrium points of the full system were constructed by applying Lyapunov’s method. The stability classification (Table 1) shows that saddle equilibria support exactly one family of periodic orbits, whereas stable (minimum) equilibria support two independent families. These results provide a local characterization of the non-degenerate equilibria and nearby periodic solutions of the system and may serve as a basis for further investigation of global orbit structure and bifurcation phenomena. From an engineering standpoint, these findings carry direct implications for analysis and design. The parameter regimes in which the Kovalevskaya exponents are complex are precisely those in which the two-mode interaction admits no second invariant. Consequently, the response in these regimes is expected to be irreducibly nonlinear, with possible dynamical consequences including persistent modal energy transfer and pronounced sensitivity to initial conditions, although these are not directly quantified here. For analogous coupled structural, electromechanical, or fluid--structure systems operating in such regimes, a reduced-order model premised on decoupled modes is inadmissible, so that robust design margins and full multi-physics numerical simulation become necessary rather than optional. The exceptional regimes in Eq. (14), together with the stable operating points identified in Section 3, delineate regions where simplified, near-integrable descriptions may potentially be applicable. However, the Kovalevskaya exponent analysis is inconclusive in these regimes and does not establish integrability, near-integrability, decouplability, or guaranteed safety for modal reduction. Therefore, these regimes warrant further detailed integrability and model-reduction analysis before any definitive claims can be made regarding their suitability for modal reduction. A natural continuation is to evaluate the present diagnostic for the positive-definite kinetic structure characteristic of mechanical realizations and to couple it with multi-physics simulation of the identified non-integrable regimes.