Non-integrability of a model of elastic dumbbell satellite

We study the integrability of a model of elastic satellite whose centre of mass moves in a circular Keplerian orbit around a gravity centre. The satellite is modelled by two point masses connected by an extensible massless spring that obeys Hooke’s law. It is assumed that the distance between point masses is much smaller than the radius of the orbit, so the orbital motion of the satellite is not perturbed by its rotational motion. The gravity potential of the satellite is expanded into a series with respect to its size up to quadratic terms which describe the gravity gradient torque acting on the satellite. Two cases are considered with Hooke’s centre localised in the centre of mass of the dumbbell and at an arbitrary point along a line connecting both masses. It is shown that the first case appears to be integrable and super-integrable for selected values of the parameter of the system. In the second case, model depends effectively only on one parameter and is non-integrable. In the proof, differential Galois integrability obstructions are used. For the considered sysem, these obstructions are deduced thanks to the recently developed symplectic Kovacic’s algorithm in dimension 4. According to our knowledge, this is the first application of this tool to a physical model.


Introduction
We consider a model of an elastic dumbbell satellite following Sidorenko and Celletti [10]. It consists of two point masses connected by a massless spring and moves in a central gravity field. We assume that the dumbbell is short; thus, its orbital motion decouples from its rotational motion. The dumbbell satellite models have attracted the attention of scientists because they are suitable for an investigation of the general properties of a rigid or elastic deformable body motion in a gravity field.
The aim of this article is the integrability analysis of this model. Two sets of coordinates: Cartesian and spherical are used to show different aspects of the dynamics. The complete integrability result for physically important ranges of parameters variability is obtained using the differential Galois theory and the direct method of searching for first integrals.
The plan of the paper is as follows. In Sect. 2, we derive equations of motion for the considered model. We use two sets of coordinates which parametrise the configuration space of the problem. Using the first one, we can easily notice that the first model with Hooke's centre located at the mass centre of the dumbbell is integrable. In Sect. 3, we show that for particular values of parameters this model is super-integrable. The second set of coordinates is used in the proof of our main theorem which states that except explicitly distinguished integrable cases the system is not integrable. So, in this paper we give necessary and sufficient conditions for integrability of the considered system. In Sect. 4, we distinguish equilibria and certain invariant sets of the system. We demonstrate also its chaotic behaviour using numerical methods. Section 5 contains a proof of our main theorem. It is based on the symplectic Kovacic algorithm recently formulated in Combot and Sanabria [1]. This algorithm can be used to investigate systems of variational equations in dimension 4. As this tool is relatively new, we collect basic theoretical facts in "Appendix".

Equations of motion
Following Sidorenko and Celletti [10], we assume that two point masses m 1 and m 2 are connected by an extensible massless spring that obeys Hooke's law. The centre of mass of this elastic dumbbell moves in an orbit around a gravity centre located at the origin of the inertial frame. We assume that the dimensions of the dumbbell are much smaller than the dimension of the orbit.
The radius vectors of masses m 1 and m 2 are The gravitational potential energy of point masses can be expanded into the Taylor series where κ = G M is the gravitational parameter of the centre, see, e.g. Sidorenko and Celletti [10]. Thus, the gravitational potential of the dumbbell reads where e r = r/|r| is the unit vector in the direction to the centre of mass. We assume that the elastic dumbbell is permanently straight, and the potential energy of its elastic deformations is given by where d 0 is the length of the undeformed spring. If d 0 = 0, then Hooke's centre is localised at the centre of mass of the dumbbell. We will discuss this case later. In a typical study of the attitude motion of a satellite, it is usually assumed that its orbital motion is not perturbed by its rotation. Up to rotation of the Cartesian coordinates, we can assume that the mass centre of the dumbbell moves in the (x, y)-plane of the inertial frame along a circular Keplerian orbit with radius a := |r|. Thus, the Lagrange function of the system is It depends explicitly on time because components of e r are time dependent. To remove this dependence, we pass to the orbital frame with the origin at the centre of mass of the dumbbell and with axes {s, t, n}, where n and t are unit vectors normal and tangential to the orbit, and s = t × n = r/r = e r , see Fig. 1. We set where q i are component of vector d with respect to the orbital frame, and ω 2 = κ/a 3 is the orbital angular velocity. Then, According to our assumptions, s = e r and ω = ω[0, 0, 1] T . In terms of (q,q) variables, the Lagrangian reads Assuming that d 0 = 0, we can rescale variables q i → d 0 q i , and time t → ωt. Then, the above Lagrangian and γ := c/(μω 2 ) is a new parameter. Let us note that γ ≥ 0 that will be used in further considerations. The standard Legendre transformation of the above Lagrangian gives the following Hamiltonian function (2.14) The corresponding Hamilton equations read (2.15) Now, we introduce new variables (ξ, ϕ, θ) useful in further analysis (2.16) These variables are well defined for θ = 0, π. The angles ϕ and θ are defined as it is shown in Fig. 2. In these variables, Lagrangian L given in (2.13) takes the form and making its Legendre transformation, we obtain the following Hamilton function (2.18) The corresponding Hamilton equations take the forṁ At this point, we explain why we derived equations of motion in two sets of coordinates in the configuration space of the system. In the Cartesian coordinates, the state of the system is described by (q,q) or (q, p) which are 'almost' global, that is (q, p) ∈ (R 3 \ {0}) × R 3 . We have to remove q = 0 from the configuration space because if cd 0 = 0, then the elastic potential is not differentiable at q = 0. Moreover, if γ = 0, then the Hamiltonian (2.14) is a homogeneous quadratic polynomial function of (q, p), so equations of motion (2.15) are linear and thus, they are integrable. However, this property is not obvious in other coordinates, see Hamiltonian (2.18) and Eq. (2.19).
On the other hand, we will use coordinates (ξ, ϕ, θ) in our proof of the main theorem of this paper, namely Theorem 2.1 If γ > 0, then the system given by Hamiltonian (2.18) is not integrable in the Liouville sense with first integrals which are meromorphic in p ϕ , p θ , p ξ , ξ, cos θ, sin θ, cos ϕ, sin ϕ.
Remark that the variable change between Cartesian coordinates and coordinates (2.16) transforms mero-morphic functions in p, q, q 2 1 + q 2 2 + q 2 3 to meromorphic functions in p ϕ , p θ , p ξ , ξ, cos θ, sin θ, cos ϕ, sin ϕ. Thus, Theorem 2.1 forbids the existence of additional first integrals of Hamiltonian (2.14) meromorphic in p, q, q 2 1 + q 2 2 + q 2 3 . The proof of this theorem is quite long and the proper choice of coordinates allows to avoid several analytical difficulties and complications.
In the end, let us underline that we derive the Hamiltonians mentioned in the above theorem under assumption that d 0 = 0. The case d 0 = 0 is itself interesting, and it is considered in the next section.

Integrable and super-integrable cases
In order to have possibility to investigate cases with d 0 = 0, we rescale only time variable t → ωt in the Lagrange function (2.12). Then, L = μω 2 L, where now The corresponding Hamilton function is Let us notice that for d 0 = 0 or γ = 0 this Hamiltonian is a homogeneous polynomial of degree two with respect to phase variables (q, p) and thus, its equations of motion are linear in variables.
In the remaining part of this section, we will consider case d 0 = 0 with simplified Hamiltonian This system is integrable with commuting first integrals  and 4b the corresponding projections of these trajectories on plane (q 1 , q 2 ). The inclination θ oscillates between −π and 0, see Figs. 3c and 4c. The time changes of the azimuth angle ϕ are more complicated as the vector q rotate and the direction of these rotations is changing more or less periodically, see Figs. 3d and 4d. The deformation parameter of the dumbbell ξ oscillates with high frequency which illustrate Figs. 3e and 4e .
If we additionally assume that γ = 0, then one can find one more first integral. Thus, in this case the system has three additional first integrals They are functionally independent together with H and satisfy the following commuting relations Thus, if γ = 0, then the system is super-integrable, but H, I 1 , I 2 , I 3 , I 4 are algebraically dependent as the following relation holds true It is natural to ask if for other values of γ the system is super-integrable. To answer this question, we observe that for a generic value of γ the system is integrable so its invariant tori are three-dimensional manifolds in the phase space. With each such a torus, we have related three periods which are independent over Z. If one more additional first integral appears, then its common level with the invariant torus is typically a torus of dimension two. So, if the system is super-integrable, then the threedimensional torus is foliated by two-dimensional tori. This is why the three periods cannot be Z independent. Thus, if the system is super-integrable, then a resonance between basic frequencies appears.
To identify characteristic frequencies of this system, let us check the characteristic polynomial of this matrix After substitution λ = iω, characteristic equation takes the form Its solutions are of the form ±ω k , k = 1, 2, 3, where Thus, if the matrix is diagonalisable, using a canonical change of variables the Hamiltonian (3.3) can be transformed to the form where now (I i , ϕ i ) are canonical action-angle variables and Clearly, I 1 , I 2 and I 3 are first integrals of the system. If there is a resonance between frequencies of the form is a first integral of the system. The number |n| := |n 1 | + |n 2 | + |n 3 | is called the order of the resonance. Instead of ψ it is more useful to consider the first integral I 4 = sin ψ (or I 4 = cos ψ). Since sin ψ and cos ψ are polynomials of degree |m| in canonical variables, the additional first integral for cases corresponding to resonance conditions (3.8) can be chosen a homogeneous and polynomial function in variables (q, p) of degree equal to the order of resonance |m|.

Below, we list triples of integers
Let us notice that for γ = 0 we have ω 1 −ω 3 = 0 and ω 2 = 0, thus |m| = 2, and this is the only resonance of order 2 with γ ≥ 0. In this case, matrix A is not diagonalisable.
In a generic Liouville integrable case, connected compact common levels of three independent first integrals are three dimension tori. In a super-integrable case the dimension of invariant tori is smaller. Let us consider a common level of four independent first integrals where h, and α i are real constants. It is a twodimensional surface in R 6 . We can visualise it in the configuration space R 3 and eliminate momenta from polynomial equations defining level M. As result, we obtain a polynomial P(q, h, α 1 , α 2 , α 3 ). Its zero level defines an algebraic surface which is the image of the invariant torus in the configuration space. Polynomial P is a product of irreducible polynomials P = P 1 · · · P k , k > 1 and P i = 0 defines a connected component of the surface. In Fig. 5, we present components of these algebraic surfaces selected by the choice of the initial condition q The trajectory for this initial condition is shown as black line on these surfaces. Surfaces shown in Fig. 5c, f are cylinders as the polynomials defining them do not depend of q 3 . Thus, for these cases invariant manifolds M are not compact.
Time evolution of inclination θ , azimuth angle ϕ and the deformation parameter of the dumbbell ξ cor-132 T. Combot et al.  6 Time evolution of θ, ϕ and ξ for selected super-integrable cases. Initial conditions: responding to these super-integrable cases is shown in Fig. 6.
If relation (3.8) holds for two non-colinear m 1 , m 2 , m 3 , the system is doubly resonant, and then it would admit two additional first integrals. As ω 1 = √ γ + 1 > 0, we can divide by it relation (3.8), and noting we obtain that both ρ 1 , ρ 2 should be rational.
Let us recall that for γ = 0 matrix A is not diagonalisable and this is why the system is super-integrable in this case with just one additional first integral.
Proof Taking both Eqs. (3.10) and eliminating γ by a resultant formula, we obtain the relation We are looking for rational solutions of this curve. Noting s = ρ 2 1 , if ρ 1 , ρ 2 is a rational solution of (3.11), then s, ρ 2 is a rational solution of This curve is an elliptic curve, and its Weierstrass form is v 2 = u 3 + 191u + 99198. Magma computational algebra system manages to compute its Mordell-Weil group, and it admits only 6 rational points, leading to the following rational solutions (s, ρ 2 ) For the last one, s is not a rational squared, thus does not lead to a rational solution of (3.11). Solution (0, 1) gives γ = 0, and the others do not lead to a nonnegative γ .
Removing the radicals in (3.8), we find that γ should satisfy the following quartic polynomial Then, substituting m 1 , m 2 , m 3 by integers gives resonant γ 's as solutions of this equation. Relation (3.8) defines a straight line with rational slope in ρ 1 , ρ 2 coordinates for each triplet of (m 1 , m 2 , m 3 ), and the ones leading to nonnegative γ 's are represented in Fig. 7 for |m| ≤ 4. Compared with a generic Hamiltonian system with three degrees of freedom not all resonances for our system are possible because the quartic equation could have no nonnegative roots (compare with Karabanov and Morozov [4]).

Basic properties of the system
We start this section with numerical examples showing behaviour in time of the dumbbell with d 0 = 0 obtained by integration of equations of motion in angular variables (2.19). Presented simulations were made for γ = 500. All simulations were performed using software Mathematica with working precision at least 13 so that precision of 13 digits has been maintained during internal computations. We chose four initial conditions given below for which the dynamics looks different Spatial motion of the dumbbell vector q(t) in the orbital frame expressed in angular variables in (2.16) for the above initial conditions is presented in Fig. 8. The corresponding time evolutions of the inclination angle θ , the azimuth angle ϕ and the deformation parameter ξ are shown in Fig. 9. One can notice that for initial conditions IC 1 oscillations of ξ have small amplitude, so the end of vector q moves approximately on a sphere; its inclination oscillates and rotates around the q 3 -axis. For remaining considered initial conditions, the length of the dumbbell ξ changes considerably and the motion of the dumbbell vector q is a superposition of oscillations of θ and ξ and rotations of ϕ.
These numerical experiments just illustrate the complexity and the variety of dynamics of the system. However, for further analytical and numerical investigations we need to identify its simplest invariant sets.
Hamilton's equations (2.15) have the following equilibria: saddle-centre-centres L 1,2 , saddle-saddlecentres S 1,2 and for elastic dumbbell satellite rigid Eigenvalues of linearization of vector field (2.15) at the respective equilibria are following • for L 1,2 : • and for O 1,2 In variables (ϕ, θ, ξ, p ϕ , p θ , p ξ ), equilibria are given by The methods used to prove that the system is not integrable require that we know a non-equilibrium solution of the system. There is no general method how to find such a solution. However, the considered system admits an invariant two-dimensional manifold on which it reduces to a one degree of freedom Hamiltonian system describing oscillations of vertically oriented satellite.
The Hamiltonian system (2.15) admits also the invariant manifold on which it is a Hamiltonian system with two degrees of freedom with Hamiltonian where q = (q 1 , q 2 ). In angular coordinates (2.16), this manifold is defined by On this manifold, variable θ becomes the polar angle in the plane, Hamiltonian of the restricted system has the form and the corresponding Hamilton equations reaḋ 3 cos(2θ)) .

(4.4)
This two degrees of freedom system on M 4 does not seem to be integrable. This is suggested in the sequence of the Poincaré cross-sections shown in Fig. 10. Following Sidorenko and Celletti [10], we take γ = 500 that corresponds to considered in this paper value of parameter β = 1 γ = 0.002. Cross-sectional plane was chosen ξ = 0, and points are generated when p ξ > 0. Angle θ is taken modulo π . These cross-sections show that for low energies the dumbbell oscillates periodically with small amplitude around the direction to the gravitational centre, and this periodic motion is stable, see Fig. 10a. Increasing energy we can achieve that besides quasi-periodic oscillations rotations are possible. Moreover, visible chaos appears, but still the periodic solution with small amplitude is stable. However, when the energy is big enough, then it vanishes and a new hyperbolic periodic solution with small amplitude appears.
These numerical examples suggest moreover that the original system with three degrees of freedom is not integrable. In fact, if the system restricted to M 4 is not integrable, then the original system is not integrable. We tried to prove directly that the restricted system is not integrable using differential Galois tools, however, in vain because we were not able to find a non-equilibrium particular solution.

Proof of Theorem 2.1
We consider the complexification of the considered system in order to apply differential Galois methods to the integrability analysis. Applied notions and results of this theory are shortly described in "Appendix".

It has first integral
and its solutions are Let [ , , , P ϕ , P θ , P ξ ] T denote variations of variables [ϕ, θ, ξ, p ϕ , p θ , p ξ ] T . Then, the variational equations along this particular solution will take the form where ξ = ξ(t) satisfies (5.1). We notice that variational equations separate into two blocks: normal variational equations for variables [ , , P ϕ , P θ ] T and tangential equations for variables [ , P ξ ] T .
We will now analyse normal variational equations. This system of four equations can be rewritten as the fourth-order equations, e.g. for variable . We choose a particular solution corresponding to h = 0 just for simplification of expressions. After the change of independent variable t → z = ξ(t), we obtain the following equation

The companion matrix of this equation is
see Eq. (A.5) for the definition of the companion matrix. We check whether fourth-order differential operator defined by Eq. (5.2) is symplectic. To this aim, we look for a skew-symmetric matrix which satisfies the matrix equation see Combot and Sanabria [1]. Common factor w 0 of all entries of this matrix was distinguished just for simplifying the result form. Equation (5.4) has a solution with the following entries Because of the presence of the square root √ z(γ z + z + 2) matrix W defines a projective symplectic structure. This structure is not-degenerated because det W = 16z 6 (z + 1) 10 This implies that differential Galois group G of differential Eq. (5.2) is a subgroup of projectively symplectic matrices PSP(4, C). In our further considerations, we use positivity of parameter γ , γ > 0. Our proof is based on the necessary conditions formulated in Lemma A.5 contained in "Appendix". Thus, at first, we show that Eq. (5.2) does not have any hyperexponential solution. Function f (z) is called hyperexponential if its logarithmic derivative f (z)/ f (z) is a rational function. If γ ∈ R + \ {0, 3 5 , 1}, Eq. (5.2) has five regular singularities: γ +1 and z 5 = ∞. The sets of exponents at the respective singularities are following where c i are roots of the following polynomial of degree 4 They are given by , .

Equation (5.2) is Fuchsian; thus, it admits a hyperexponential solution if and only if it has the form
singularities, e i are exponents at z i , and there exists an exponent at infinity e ∞ such that the sum i e i +e ∞ is a non-positive integer, see, e.g. Singer and Ulmer [12]. Taking into account exponents at finite singularities one of roots of (5.5) must be a half-integer. When we substitute Z = n 2 to Eq. (5.5), then we obtain quadratic equation for γ In order to have at least one real root, the discriminant of this equation should be non-negative that gives condition Q(n) = 16(−15n 4 + 120n 3 − 272n 2 +128n + 256) ≥ 0.
As lim n→±∞ Q(z) = −∞, function Q(n) takes non-negative values only for a few integer n contained in the interval 2 If γ is different from the above-mentioned values, then the differential Eq. (5.2) does not admit a hyperexponential solution. On the other hand, for these specific values of γ the equation does not depend on parameters and using a computer algebra system, e.g. Maple, we can check that it does not have any hyperexponential solution.
To check the second assumption of Lemma A.5, we can transform the system (5.4) to an equation of sixth order and check how many hyperexponential solutions it has. For computation simplifications, it is better to calculate the second exterior power of Eq. (5.2) and check how many hyperexponential solutions it admits. The second exterior power of Eq. (5.2) is a sixth-order differential equation with eight regular singularities: γ +1 , and z 4 , z 5 , z 6 , z 7 which are roots of the fourth-order equation where γ i for i = 1, . . . , 11 are roots of the following equation The respective sets of exponents at singularities are the following: where d i are roots of equation For the second exterior power to have exponential solutions, the sum 8 i=1 e i , e i ∈ E i should be a nonpositive integer. We have one choice of exponents (e 1 , e 2 , e 3 , e 4 , e 5 , e 6 , e 7 , e 8 ) = (−2, −1/2, −1/2, 0, 0, 0, 0, 3) that sum to 0. Really, the second exterior power has a solution for every γ of the form In order to find other solutions, one can built them only using a root of Eq. (5.9). Looking on exponents at finite points, we deduce that such a root must be a halfinteger. Substitution Z = n 2 to (5.9) gives the following second-order equation for γ (n − 10)(n − 6) 2 (n − 2)γ 2 + 2 ((n − 12)n((n − 12)n The discriminant of this polynomial is 1024((n − 12)n((n − 12)n + 61) + 964) ≥ 0; thus, it has real roots γ 1 and γ 2 . Both of them are negative if simultaneously two inequalities due to Vieta's formulas hold where γ 1 , . . . , γ 11 satisfy Eq. (5.8), the differential Galois group contains the whole SP(4, C) that is not solvable and thus in particular is not Abelian. These selected values of γ correspond to confluence of singularities of normal variational Eq. (5.2) or of its second exterior power, or possibility of existence of an exponential solution of normal variational equation or of its the second exterior power. For these values γ , the system no longer depends on parameters and can then be treated with symbolic software as, e.g. Maple.
We obtain that for all these selected cases normal variational equation does not factorise thus its differential Galois group is not reducible. The second exterior power also does not factorise and has just one hyperexponential solution even for these special values of γ . Thus, differential Galois group of (5.2) is not irreducible. Therefore, the system is not integrable for all values γ > 0.
Roughly speaking, differential Galois group of variational equations is a linear algebraic group which preserves polynomial relations between solutions of these linear equations. As a linear algebraic group in general, it consists of several connected components in the Zariski topology, and among them the one containing identity is called the identity component. For a precise definition of the differential Galois group and differential Galois theory, see, e.g. Kaplansky [3] or Magid [7].
The presence of n commuting first integrals of a Hamiltonian system with n-degrees of freedom imposes a very special property of the differential group of variational equation that gives the necessary integrability condition formulated in Morales Ruiz [8].
Theorem A.1 Assume that a Hamiltonian system is meromorphically integrable in the Liouville sense in a neighbourhood of an analytic phase curve . Then, the identity component of the differential Galois group of the variational equations along is Abelian.
Applications of this theorem are hard because of difficulties related to • finding a non-equilibrium particular solution, • determination of differential Galois group of variational equations. Difficulty grows very quickly with dimension of the variational equations.
In applications, it happens that variational equations split into tangential and normal subsystems and one can restrict to analysis of differential Galois group of normal variational equations that have smaller dimension than full variational equations, for detailed explanation, see, e.g. in Morales Ruiz [8] or shortly in Maciejewski and Przybylska [6]. In the case, when the dimension of normal variational equations is two and coefficients of these equations are rational functions of an independent variable (sometimes after appropriate algebraic change of independent variable) there is a very powerful tool for determination of the differential Galois group. This tool called the Kovacic algorithm formulated in Kovacic [5] is a decision procedure in a finite number of steps if one can write a general solution of differential equation in a closed form, more precisely in a field of Liouvillian functions. This algorithm is constructive, it means that it enables to determine explicitly the form of solutions.
As the byproduct, it determines the differential Galois group of Eq. (A.1). In fact, it is built on the base of complete classification of algebraic subgroups of the group SL(2, C). Unfortunately, there is no an equivalent of the Kovacic algorithm for linear differential equations with rational coefficients of higher orders although many partial results especially for equations of the third and fourth orders are known, see, e.g. Singer and Ulmer [11,12], van Hoeij et al. [14], Ulmer [13] and references therein.
But in the case of Hamilton equations that can be written aṡ where I n is n-dimensional identity matrix, their variational equations in a neighbourhood of a particular solution ϕ(t) are also Hamiltoniaṅ 3) Here, variables Y are variations of z. Solutions of these variational equations have very special property. Namely, if y 1 and y 2 are arbitrary two solutions, then quantity ( y 1 , J y 2 ) is constant. But this implies that for arbitrary element G of differential Galois group holds const = ( y 1 , J y 2 ) = (G y 1 , JG y 2 ) = ( y 1 , G T JG y 2 ) and this implies J = G T JG. But this equality means that the differential Galois group of VEs is a subgroup of Sp(2n, C). Group Sp(2, C) is isomorphic to SL(2, C) but in higher dimensions Sp(2m, C) ⊂ SL(2m, C), for m > 1 and Sp(2m, C) is much smaller than SL(2m, C), thus classification of its subgroups should be much simpler. Using this classification, one can try to construct the equivalent of the Kovacic algorithm for symplectic differential operators and for dimension four this was formulated by Combot and Sanabria [1].
To describe it shortly, we introduce appropriate terminology. Let L be differential operator with coefficients in C(z) L(y) = y (n) + a n−1 y (n−1) + · · · + a 1 y + a 0 y = 0, a i ∈ C(z), (A.4) and A is its corresponding companion matrix (A.5) The 2n-th-order operator is • symplectic if there exists an invertible skew-symmetric matrix W with elements in C(z) which is a solution of following equation • projectively symplectic if there exists an invertible skew-symmetric matrix W with elements in C(z) which is a solution of equation for a certain λ ∈ C(z).
An operator L of order 2n is symplectic (respectively, projectively symplectic) when its Galois group is isomorphic to a subgroup of symplectic matrices Sp(2m, C) (respectively, of projectively symplectic matrices PSp(2m, C)) defined as If L is projectively symplectic, then up to a multiplication of a hyperexponential function, one can ensure that the operator is symplectic. Function f (z) is called hyperexponential if its logarithmic derivative f (z)/ f (z) is a rational function. This lemma can be proved by a direct check. It follows that if symplectic operator L has right factor L 1 , then its adjoint L * has a right factor L 1 of the same order. Thus, L has right factor L 1 , then it has also a left factor of the same degree. With a systemẋ = Ax, we can associate its external second power. It is a system of the forṁ for an antisymmetric matrix W , see "Appendix" in [2]. Thus, Eq. (A.6) is an equation for the external square of a dual to the system x = Ax. The classification theorem formulated in Combot and Sanabria [1] is the following. Lemma A.2 A Lie subgroup of Sp(4, C) is up to conjugacy generated by elements of the form:

upper block triangular matrices with diagonal
blocks of size at most 2 × 2, 2. 2×2 diagonal matrices and anti-diagonal matrices ⎡

full group Sp 4 (C).
This classification is constructed on the base of known classification of Lie subgroups of wider unimodular group SL(4, C) ⊃ Sp(4, C). Subgroups of the projective symplectic group are central extensions of these and so contain multiples of the identity matrix with non-unit determinant. But as these commute with all matrices, the possible structures of subgroups in items 1, 2 are unchanged. The next two lemmas characterise reducible case.

Lemma A.3
Let us consider the following block diagonal systeṁ where B, C and D are 2 × 2 matrices with rational coefficients. Then, Eq. (A.6) has the following particular solution where r is a rational function.
Lemma A.4 Let us consider differential operator where a, b, c, d ∈ C(z) and let A be its companion matrix. Then, Eq. (A.6) has the particular solution and det W = 0.
Both lemmas can be proved by a direct check.
In our proof of the non-integrability, we will use the following criterion. Lemma A.5 Assume that equation L(y) = y (4) + a 3 (z)y + a 2 (z)y is projectively symplectic and A is its corresponding companion matrix. If Eq. (A.13) does not admit a hyperexponential solution and Eq. (A.6) has exactly one hyperexponential solution, then the differential Galois group of (A.13) contains Sp(4, C).
Proof If (A.13) does not admit a hyperexponential solution, then it does not admit a right factor of order one. Thus, if it is reducible, it admits a right factor of order two. Hence, using Lemma A.4, Eq. (A.6) has a hyperexponential solution such that the corresponding antisymmetric matrix W 1 is singular. Equation (A.13) is projectively simplectic, and thus, Eq. (A.6) admits a hyperexponential solution such that the corresponding antisymmetric matrix W 2 is not singular. Thus, W 1 , W 2 cannot be the same, but this is impossible as by assumption (A.6) has exactly one hyperexponential solution. It implies that Eq. (A.13) is irreducible. If it admits Liouvillian solutions, then it admits at least two projective Poisson structures. But this is impossible as by assumption there exists exactly one hyperexponential solution of (A.6).
Item 1 in Lemma A.2 corresponds to reducible cases. Differential Galois group of equation L(y) = y (4) + a 3 (z)y + a 2 (z)y +a 1 (z)y + a 0 (z)y = 0, = d dz (A.14) is a reducible subgroup of Sp(4, C) if and only if the equation itself can be factorised. If a symplectic operator has a left factor, then it has a right factor of the same order. Thus for an order four symplectic operator only the following factorizations into factors with appropriate orders are possible: (i) 1,1,1,1 or (ii) 1,2,1 or (iii) 2,2.
If this equation has a factor of order one, then it is either a right factor of order one or a left factor of order one. Since a left factor of order one leads to a right factor of order one of the adjoint differential operator L * (y), then one can restrict to testing of a right factor of order one of L(y) or of L * (y). If an equation has a right factor of order one, then it has an hyperexponential solution. If L(y) = 0 is of Fuchsian type, then any hyperexponential solution must be of the form P(z) i (z − z i ) e i , where P(z) ∈ C[z], z i ∈ C are singularities except the infinity, e i are exponents at z i . Then, using lemma 3.1 in Singer and Ulmer [12] about the necessary condition for such a solution one can formulate the following proposition Proposition A.1 If a Fuchsian symplectic equation of order four has a factor of order one, then L(y) = 0 has a solution of the form P(z) i (z − z i ) e i , where P(z) ∈ C[z], z i ∈ C are singularities, e i are exponents at z i , and there exists an exponent at infinity e ∞ such that the sum i e i + e ∞ is a non-positive integer.
In the case when L(y) factorises into two factors L = L 1 L 2 of order two to each of them one can use the Kovacic algorithm. If one of them is not solvable, then the differential Galois group is not solvable.
Item 2 in Lemma A.2 corresponds to irreducible solvable cases. Then, operator L admits a factorisation in two operators of order two with coefficients in a quadratic extension of C(z). This case can be identified by the presence of two linearly independent rank two Poisson structures (requirement as for symplectic structure but without invertibility requirement) with coefficients in a quadratic extension of C(z), see Combot and Sanabria [1]. But the presence of two structures is related to the presence of two hyperexponential solutions of the exterior power of L. Thus, to identify the irreducible solvable cases we use the following proposition formulated in Combot and Sanabria [1] Proposition A. 2 The necessary condition for irreducible solvable cases is that the second exterior power of differential equation has at least two hyperexponential solutions.
The second exterior power of a differential operator L with solutions {y 1 , . . . , y n } is the equation having for its space of solutions the vector space generated by the 2×2 Wronskians of solutions y i y j y i y j , i < j of operator L, see, e.g. Schwarz [9]. To identify these solutions for Fuchsian operators, L one can use Proposition A.1. Item 3 in Lemma A.2 corresponds to non-solvable cases.