C T Bop

Y Kalugina

F Lique

Rotational excitation of NS + du H 2 -impact Rotational excitation of NS + by H 2 revisited: a new global potential energy surface and rate coefficients

Rotational excitation of NS+ by H-2 revisited: A new global potential energy surface and rate coefficients

C. T. Bop, Y. Kalugina, F.

I. INTRODUCTION

Recent observations of sulphur-bearing molecules in photodissociation regions (PDR) such as the Horsehead nebula 1 have revived investigations on the interstellar sulphur chemistry, which is still poorly understood. Indeed, in most of the astrophysical environements where S-bearing species are observed, the sum of observed S-species does not account for the sulfur cosmic abundance. The S-reservoir(s) question is still an unresolved important question.

Despite holding a small percentage of the S-budget, NS + plays a key role in the of sulphur chemistry. For instance, most of charge exchange and ion-neutral reactions of S-bearing molecules lead to NS + . 1 In addition, this cation is claimed to be a valuable PDR tracer since it is more abundant in the PDR than in the core of the Horsehead nebula. 1 This enhanced abundance towards the outer layers of the Horsehead nebula is a signature of the dependence of the NS + abundance on the external UV field. Moreover, being present at different evolutionary stages of dark molecular clouds, NS + is considered as good barometer and thermometer in these regions. 2 The large variation of the NS/NS + abundance ratio observed between cold molecular clouds (30 -50) and hot cores (> 500) suggests that the chemistry underlying the formation and destruction of these species varies from one medium to another. 2 In molecular clouds, NS + likely forms from reactions between atomic nitrogen and SO + /SH + . 2 In external layers of PDRs, i.e. UV-irradiated regions, S + + NH, SH + + N and H 2 S + + N are expected to be the favorite formation routes of NS + . 1 In addition, charge exchange reactions between NS and C + /S + are less critical in PDRs than in cold molecular clouds.

Therefore, constraining the chemistry of NS and NS + in cold molecular clouds and PDRs is crucially needed. For B1b molecular cloud, the chemical model of Cercnicharo et al. 2 reproduced the observed NS/NS + abundance ratio within a factor of ∼ 2. However, this finding must be treated with caution since the authors highlighted that they used approximate rate coefficients to interpret the observational spectra of NS + and NS.

Accurate abundance modeling requires prior determination of collisional rate coefficients. Without these data, local thermodynamic equilibrium (LTE) conditions which are rarely verified in space, would be considered. This approach has led to underestimate the abundance of NS by a factor of 1.5 -2.5. 2 Since collisional rate coefficients of NS and NS + were not available in the leterature, Cernicharo et al. 2 used the rate coefficients of CS 3 for both species. They neglected the fine structure of NS and used a scaling factor of 5 to take into account the ionic nature of NS + .

Since its interstellar detection, NS + has become the target of rotational energy transfer studies induced by collisions. Trabelsi et al. 4 reported rate coefficients for the 11 low-lying rotational levels of this cation due to helium-impact for temperatures ranging up to 300 K. The calculations were based on a three dimensional potential energy surface (which was thereafter averaged using the NS + ground vibrational wave functions) computed at the explicitly correlated coupled cluster method with single, double and non-iterative triple excitation in conjunction with the augmented-correlation consistentpolarized valence triple zeta basis (hereafter denoted as CCSD(T)-F12/aug-cc-pVTZ). Few months later, Cabrera et al. 5 used a CCSD(T)/aug-cc-pV5Z 3D PES (a similar accuracy to the PES mentioned above) to revisit the NS + -He col-ACCEPTED MANUSCRIPT / CLEAN COPY lisional system. They extended the rotational basis up to the 27 low-lying energy levels and resolved the hyperfine structure since such emission lines were also reported. 2 Both authors compared their results with the collisional rate coefficients of CS-He 3 used in the NS + abundance modeling. Despite the similarity of their data, Trabelsi et al. 4 validated the use of CS as template for NS + while Cabrera et al. 5 highlighted the limits of such an approximation.

Due to the fact that He is seen as a bad template for para-H 2 (p-H 2 ) in the case of charged species, Bop 6 reassessed the hyperfine rate coefficients corresponding to the 23 low-lying rotational levels of NS + for temperatures up to 100 K, but using p-H 2 ( j 2 = 0) as projectile. These calculations were based on a 4D PES, computed at the CCSD(T)-F12/aug-cc-pVTZ level of theory, which was thereafter averaged over the orientations of H 2 (denoted hereafter as model treatment). The comparison with the previous rate coefficients of NS + led to deviations of 2-4 depending on the transition and on the temperature.

The strong anisotropy of the PES of NO + -H 2 , 7 an isoelectronic analogue of NS + -H 2 , with respect to the orientations of H 2 let us reconsider the use of an average potential for the scattering of NS + with H 2 . In practice, we revisited the scattering of NS + due to p-H 2 impact using a full dimensional PES. Since this molecular cation was also observed towards warm/hot environments such as the Horsehead nebula, 1 where the ortho-to-para-H 2 ratio is not negligible, rate coefficients by collision with both ortho-H 2 (o-H 2 ) and p-H 2 are needed.

This paper is structured as follows: the computational details are briefly given in Section II, the results are shown and discussed in Section III and concluding remarks are displayed in Section IV.

II. COMPUTATIONAL DETAILS A. Potential energy surface and analytical fit

The potential underlying the interaction between NS + ( 1 Σ + ) and H 2 ( 1 Σ + g ), in their ground electronic state, was calculated using the Jacobi coordinate system. The PES was calculated under the assumption of rigid rotor-rigid rotor interaction using the bond lengths averaged with the ground vibrational wave functions, <r H 2 > 0 = 1.449 a 0 and <r NS + > 0 = 2.7189 a 0 . 2,8 In practice, 4 degrees of freedom are needed to describe the geometry of the collision complex (see Fig. 1): the intermolecular distance (R), the polar angles of NS + and H 2 relative to the R vector (θ NS + and θ H 2 , respectively) and the dihedral angle between the half-planes that contain NS + and H 2 (φ ).

In this work, we used as Bop 6 the state-of-the art CCSD(T)-F12/aug-cc-pVTZ level of theory 9,10 implemented in the MOLPRO molecular package version 2010 11 for all ab initio calculations. We improved the PES Bop 6 by (i) fully describing the H 2 rotations and (ii) explicitly treating the long-range interaction potential using a multipolar expansion. In fact, the previous PES of the NS + -H 2 complex was averaged over 3 H 2 rotations : (θ H 2 , φ ) = (0, 0), (90, 0) and (90, 90). The basis set superposition error correction was taken into account using the counterpoise method. 12 The interaction energy V (R,

θ NS + , θ H 2 , φ ) was expanded in bispherical harmonics A L 1 L 2 L (θ NS + , θ H 2 , φ ): A L 1 L 2 L (θ NS + , θ H 2 , φ ) = 2L 1 + 1 4π L 1 L 2 L 0 0 0 ×P L 1 0 (θ NS + )P L 2 0 (θ H 2 ) + 2 min(L 1 ,L 2 ) ∑ M=1 L 1 L 2 L M -M 0 ×P L 1 M (θ NS + )P L 2 M (θ H 2 )cos(Mφ ) (1)
where

L = |L 1 -L 2 |, ..., L 1 + L 2 ,
V (R, θ NS + , θ H 2 , φ ) = ∑ L 1 L 2 L v L 1 L 2 L (R)A L 1 L 2 L (θ NS + , θ H 2 , φ ).
(2) The expansion coefficients v L 1 L 2 L (R) at each point R (31 Rgrid points from 4.5 to 29 a 0 ) were obtained using 11-point (5-point) Gauss-Legendre quadrature for θ NS + (θ H 2 ) and 9 Gauss-Chebyshev quadrature for φ . Angles, θ NS + , θ H 2 and φ were set to [11.98, 27.49, 43.10, 58.73, 74.36, 90.00, 105.64, 121.27, 136.90, 152.51, 168.02], [13.12, 30.11, 47.20, 64.32, 81.44] and [10, 30, 50, 70, 90, 110, 130, 150, 170], respectively. The use of these quadrature allows us to expand the potential in the coupled basis with L max 1 = 10 and L max 2 = 6 which gives us a total of 142 radial expansion coefficients v L 1 L 2 L (R).

ACCEPTED MANUSCRIPT / CLEAN COPY

The long-range part of the potential can be represented using analytical multipolar expansion through the order of R -6 . 13

V (R, θ NS + , θ H 2 , φ ) = 1 3 T αβ q A Θ B αβ + 1 105 T αβ γδ q A Φ B αβ γδ - 1 2 T α T β (q A ) 2 α B αβ + T α T β γ q A α B αβ µ A γ - 1 2 T αγ T β δ α B αβ µ A γ µ A δ - 1 3 T αβ γ Θ B αβ µ A γ + 1 9 T αβ γδ Θ B αβ Θ A γδ - 1 45 T αβ γδ ε Θ B αβ Ω A γδ ε - C 0 6 6α A α B T αβ T γδ α A αγ α B β δ , (3)
In practice, 100 Gauss-Legendre quadrature nodes are considered to determine the isotropic term, C 0 6 = 40.3027 E h a 6 0 , using the (equation-of-motion) EOM-CCSD method 14 implemented in MOLPRO. The average polarizability α and the tensor T αβ γ...ν are defined as 1/3(α xx + α yy + α zz ) and ∇ α ∇ β ∇ γ ...∇ ν R -1 , respectively. Superscrips A and B denote molecules NS + and H 2 , respectively. We implied in Eq. 3 the dependence of multipole moments and polarizabilities of molecule A on angle θ NS + and of molecule B on angles θ H 2 and φ . The values of multipole moments and polarizabilities used in analytical calculations are presented in Table I. These values were calculated at the CCSD(T)/aug-cc-pVTZ level of theory using finite-field method. 15 We noticed that for very large intermolecular separations, CCSD(T)-F12a method started to deviate from the analytical long-range approximation. Therefore, we performed further calculations using standard CCSD(T) along with the augcc-pVTZ basis set. Indeed, at intermediate distances R = 20 -22 a 0 , the CCSD(T) and CCSD(T)-F12a computed energies agree very well whereas discrepancies occur at R > 29 a 0 . Using the long-range analytical potential as reference at R = 100 a 0 , we obtained slight deviations of ∼ 10 -5 % for standard CCSD(T) whereas CCSD(T)-F12a leads to relative errors of up to 30%. In fact, the agreement with CCSD(T) is not surprising since the multipole moments and polarizabilities of monomers (Eq. 3) are computed at the same level of theory. Therefore, to obtain a smooth and accurate interaction potential, we used CCSD(T)-F12a for R ≤ 20 a 0 , CCSD(T) for 20 < R ≤ 25 a 0 and the analytical multipolar expansion for R > 25 a 0 .

The root-mean-square error (rmse) of the analytical fit (Eq. 2) is about [10 -4 -10 -3 ] cm -1 and [10 -2 -10 -1 ] cm -1 in the attractive and repulsive regions, respectivly. The expansion coefficients obtained from the fit, i.e. for 4.5 ≤ R ≤ 29 a 0 , were interpolated using cubic spline routine. For R > 29 a 0 , the potential is represented solely by analytical long-range multipolar expansion.

B. Cross sections and rate coefficients

The scattering processes we are studying can be summarised as follows:

NS + ( j 1 ) + o-; p-H 2 ( j 2 ) → NS + ( j 1 ) + o-; p-H 2 ( j 2 )(4)
where j 1 and j 2 are the NS + and H 2 rotational quantum numbers, respectively. The hyperfine structure due to the nuclear spin of nitrogen is not considered in this work but it can be derived from our data using the Infinite Order Sudden (IOS) approximation. 16,17 The cross sections (σ ) were calculated for the 15 low-lying rotational levels of NS + ( j 1 = 0 -14) for total energies ranging up to at least 500 cm -1 . The energy range was spanned using a fine step size of 0.1 cm -1 to correctly describe the resonances. We used the close-coupling quantum mechanical approach 18 implemented in the MOLSCAT code. [START_REF] Hutson | Collaborative computational project[END_REF] The log derivative-airy propagator was used to solve the coupled equations. [START_REF] Alexander | [END_REF] Preliminary test calculations were performed to set the integration parameters. For NS + , we considered a quite wide rotational basis which includes at least 10 channels above the considered rotational energy levels. Typically, the size of the basis was varied from 23 to 30 rotational levels ( j 1 = 22 -29).

Regarding the H 2 rotational basis, only j 2 = 0 and j 2 = 1 were used for pand o-H 2 colliders, respectively. For p-H 2 , convergence was checked by comparing the sum of partial cross sections, σ ( j 2 = 0) and σ ( j 2 = 0 -2), over the 11 first total angular momenta (J = 0 -10). For total energies of 100 cm -1 , 250 cm -1 and 500 cm -1 , the deviation of σ ( j 2 = 0) with respect to σ ( j 2 = 0 -2) is less than 30%. In the case of o-H 2 , we calculated partial cross sections σ ( j 2 = 1) and σ ( j 2 = 1 -3) up to J = 5 to reduce the computational cost due to the inclusion of j 2 = 3. As shown for HCCNC and HNCCC, 21 the more total angular momentum we use, the better the convergence would be. Since the first rotational energy level of H 2 is shifted to ∼ 120 cm -1 above the ground level, we used total energies of 220 cm -1 , 370 cm -1 and 620 cm -1 . The disagreement between σ ( j 2 = 1) and σ ( j 2 = 1 -3) is less than 20% except at low energy. For a better appreciation of the convergence, we show in Table II the variation of cross sections as a function of the increment of the H 2 rotational manifold.

The STEPS-parameter was varied from 60 to 10 to keep the integration step small enough. The rotational energy levels were constructed using the spectroscopic constants of NS + [B 0 = 0.835 cm -1 and D 0 = 1.17 × 10 -7 cm -1 ] and H 2 [B 0 = 59.322 cm -1 and D 0 = 0.047 cm -1 ]. 2,8 By thermally averaging the so-calculated cross sections over the Maxwell-Boltzmann velocity distribution, we derived collisional rate coefficients (k) for temperatures up to 50 K,

k j 1 , j 2 → j 1 , j 2 (T ) = ( 8 π µβ ) 1/2 β 2 ∞ 0 E k σ j 1 , j 2 → j 1 , j 2 (E k )e -β E k dE k (5)
where E k is the kinetic energy, µ = 1.931 au is the reduced mass of the complex and β the inverse of the Boltzmann constant multiplied by the temperature. TABLE II. The variation of cross sections ( Å2 ) as a function of the increment of the H 2 rotational manifold. For p-H 2 and o-H 2 , we summed partial cross sections over the 11 (J = 0 -10) and 6 (J = 0 -5) first total angular momenta, respectively.

ACCEPTED MANUSCRIPT / CLEAN COPY

p-H 2 o-H 2 E = 100 cm -1 E = 250 cm -1 E = 500 cm -1 E = 220 cm -1 E = 370 cm -1 E = 620 cm -1 j 1 j 1 j 2 = 0 j 2 = 0 -2 j 2 = 0 j 2 = 0 -2 j 2 = 0 j 2 = 0 -2 j 2 = 1 j 2 = 1 -3 j 2 = 1 j 2 = 1 -3 j 2 = 1 j 2 = 1 -3 1 0 4.01×10 -1 5.20×10 -1 1.03×10 -1 1.38×10 -1 6.72×10 -2 6.02×10 -2 4.16×10 -1 5.85×10 -1 9.76×10 -2 9.20×10 -2 4.14×10 -2 4.87×10 -2 3 2
2.54 2.93 5.08×10 -1 6.05×10 -1 1.83×10 -1 1.85×10 -1 9.45×10 -1 9.82×10 -1 1.74×10 -1 1.69×10 -1 6.22×10 -2 6.44×10 -2 7 5

5.73 5.48

8.14×10 -1 7.55×10 -1 2.09×10 -1 1.67×10 -1 1.08 1.19 2.23×10 -1 2.26×10 -1 5.52×10 -2 5.54×10 -2 9 7 1.11×10 1 1.12×10 1 9.55×10 -1 9.52×10 -1 2.27×10 -1 2.34×10 -1 1.16 1.29 1.83×10 -1 1.97×10 -1 4.48×10 -2 4.75×10 -2

III. RESULTS

We display in Fig. 2 bi-dimensional cuts of the 4D PES of the NS + -H 2 complex through the global minimum (see below). The potential exhibits a strong anisotropy except in panel (c) where it depends weakly on φ . The global minimum corresponds to a planar structure where the axes of the H 2 and NS + molecules form an acute angle. Typically, we obtained

θ NS + = 124.3 • , θ H 2 = 81.8 • , φ = 180.0 • and R = 5.77 a 0 with ∆E = -849.88 cm -1 .
The potential well depth is thus similar to that of Bop 6 which is located at -848.24 cm -1 below the dissociation limit ACCEPTED MANUSCRIPT / CLEAN COPY of the complex. With respect to the interaction potentials where helium was used as projectile, 4,5 the new global minimum is about 5 times deeper. We plot in Fig. 3 the variation of the potential as a function of θ H 2 at various radial distances. As one can see, the PES remains sensitive to the H 2 orientations at long-range. For example at R = 45 a 0 , the potential varies by nearly 2 cm -1 . In fact the long-range is led by charge-quadrupole interactions which slowly decay as R -3 and sensitive to the orientation of the quadrupole, i.e. H 2 . In contrast, we notice very slight variations of the PES relative to θ NS + for R > 20 a 0 . This comes from the fact that, charge-quadrupole interactions do not depend on θ NS + at long-range and the other higher-order interactions (dipole-quadrupole and quadrupole-quadrupole) remain weak relative to charge-quadrupole.

We show in Fig. 4 the variation of the rotational cross sections of NS + in collision with both oand p-H 2 as a function of the kinetic energy. The amplitude quickly drops with increasing kinetic energy-logarithm, exhibiting the so-called Langevin behaviour which is known for ion-molecule collisional cross sections. In fact, for ion-molecule collisional systems, the cross sections are inversely proportional to the square root of the collision energy. Moreover, both projectiles lead to very close cross sections with slight deviations below 10 cm -1 . This behavior is a typical feature of molecular ions. 22 Compared to the previous calculations of Bop, 6 the shape of the cross sections does not change significantly. For a better assessment of the improvements we carried on the potential, we show in Fig. 5 a systematic comparison of the NS + rotational cross sections induced by p-H 2 computed in this work and those reported by Bop 6 at 100 cm -1 , 250 cm -1 and 500 cm -1 . The cross sections derived from the full 4D PES agree mostly within a factor of 2 with the data mentioned above. A comparison using a full energy range would show even more discrepancies than the few points seen in the top panel of Fig. 5. Indeed, the use of two different PESs will lead to a shift of the locations of the resonances, i.e. at low kinetic energy, where the cross sections are very sensitive to the PES. However, it is worth noting that the differences obtained are surprising since reasonable agreement are usually found between full and reduced dimensional approaches for collisions involving ions. 23,24 To figure out the origin of these discrepancies, we averaged the new potential over the H 2 rotation (i.e. (θ H 2 , φ ) = (0, 0), (90, 0) and (90, 90) as previously done by Bop. 6 The bottom panel of Fig. 5 shows a perfect agreement clearing out any suspicion on the different coordinate systems used in the construction of the PESs. From this comparison, we show the need for PESs describing all the orientations of H 2 in the case of strongly anisotropic interactions.

To settle upon the ambiguity highlighted in the introduction, i.e. whether CS is a good template for NS + or not, we compare in Fig. 6 our new collisional data obtained for the p-H 2 collider to the CS-He rate coefficients. 25 Despite being scaled up by a factor of 5, the CS-He rate coefficients 25 underestimates the NS +p-H 2 ones by nearly an order of magnitude for most transitions. The discrepancy is even greater at 10 K, i.e. the typical temperature of cold media where NS + was observed. 2 These strong deviations are not surprising since scaling up rate coefficients of the neutral CS molecule cannot efficiently incorporate the ionic nature of the NS + cation, i.e. the Langevin behavior seen in the cross sections. Besides the different shapes, Fig. 6 shows that a factor of at least ACCEPTED MANUSCRIPT / CLEAN COPY 10 would have been needed to describe reasonably the dominant rate coefficients (k > 10 -11 cm 3 s -1 ). The need of such a high multiplier is directly related to the interaction potentials. In fact, the CS-He potential 25 has a well depth of ∼ 22 cm -1 which is very small compared to that of NS + -H 2 (∼ 850 cm -1 ).

IV. CONCLUSION

The first 4D PES of the NS + -H 2 van der Waals complex was calculated using the CCSD(T)-F12/aug-cc-pVTZ level of theory. The long-range interaction was explicitly described using analytical multipolar expansion. Rotational cross sections of NS + induced by collision with ortho-H 2 ( j 2 = 1) and para-H 2 ( j 2 = 0) were derived, for total energies up to at least 500 cm -1 , using the close-coupling approach. We considered transitions between the 15 low-lying rotational energy levels and we retrieved downward rate coefficients, for temperatures up to 50 K, by thermally averaging the cross sections. The collision data obtained using ortho-H 2 ( j 2 = 1) and para-H 2 ( j 2 = 0) as projectiles are very similar suggesting the possibility of using only p-H2 collisional data in radiative transfer calculations. The explicit description of the anisotropy of the PES with respect to the orientations of H 2 led to deviations of a factor of 2 in the cross sections. Concerning the use of CS as template for NS + , we showed that a scaling factor cannot incorporate the ionic nature.

The data presented in this work, supplemented by the NS-He rate coefficients briefly described in the recent paper of Hily-Blant et al., 26 will be of great importance to better constrain the evolution of the NS + /NS abundance ratio in cold dense clouds.

6

 6

FIG. 1 .

 1 FIG. 1. Coordinate system of NS + -H 2 complex.

FIG. 2 .

 2 FIG. 2. Contour plots of 2D cuts of the NS + -H 2 4D PES. The panels (a), (b), (c) and (d) showV (R, θ NS + , θ H 2 = 81.8 • , Φ = 180 • ), V (R, θ H 2 , θ NS + = 124.3 • , Φ = 180 • ), V (θ NS + , Φ, R = 5.77 a 0 , θ H 2 = 81.8 • ) and V (θ NS + , θ H 2 , R = 5.77 a 0 , Φ = 180 • ), respectively.

1 )FIG. 3 .

 13 FIG.3. The dependence of the potential energy with respect to θ H 2 at various radial distances. The angles θ NS + and φ are fixed at 124.3 • and 180 • , respectively.

FIG. 4 .

 4 FIG.4. Variation of the rotational cross sections of NS + induced by collision with oand p-H 2 ( j 2 = 0) as a of the kinetic energy for selected transitions.

FIG. 6 .

 6 FIG.6. Systematic comparison of rotational rate coefficients (in cm 3 s -1 ) for selected temperatures. The x-axis shows the CS-He collisional rate coefficients scaled by a factor of 5 and the y-axis represents the NS + -p-H 2 data. The dashed lines delimit the region where the agreement is better than a factor of 2.

TABLE I .

 I Electric properties of the NS + and H 2 monomers computed at the CCSD(T)/aug-cc-pVTZ level of theory. All values are in a. u.

	Param.	Definition	NS +	H 2
	µ z	dipole moment	-0.7805	0.
	Θ zz	quadrupole moment	1.0029	0.4866
	Ω zzz	octupole moment	-0.7002	0.0
	Φ zzzz	hexadecapole moment	-7.79	0.36
	Φ xxxx		-2.92	0.14
	α xx	dipole polarizability	13.8399	4.7608
	α zz		23.4455	6.7291

ACKNOWLEDGEMENTS

The authors acknowledge the European Research Council (ERC) for funding the COLLEXISM project, the Programme National "Physique et Chimie du Milieu Interstellaire" (PCMI) of Centre National de la Recherche Scientifique (CNRS)/Institut National des Sciences de l'Univers (INSU) with Institut de Chimie (INC)/Institut de Physique (INP) cofunded by Commissariat a l'Energie Atomique (CEA) and Centre National d'Etudes Spatiales (CNES). F.L. acknowledges the Institut Universitaire de France. Results of ab initio calculations and the fit of the PES have been obtained under support of the RSF grant No. 17-12-01395.

CONFLICTS OF INTEREST

There are no conflicts of interest to declare.

SUPPLEMENTARY MATERIAL

A Fortran routine for generating the potential is available as supplementary material.

DATA AVAILABILITY

The data underlying this manuscript are available from the corresponding author upon reasonable request.

ACCEPTED MANUSCRIPT / CLEAN COPY


Sasho Ivanov Popov
email: sashopopov2003@yahoo.com

Witold Respondek

Jean-Marie Strelcyn
email: jean-marie.strelcyn@univ-rouen.fr

On rational integrability of Euler equations on Lie algebra so(4, C), revisited

Keywords: Euler equations, Rational and polynomial first integrals, Darboux polynomials 2 6

We consider the Euler equations on the Lie algebra so(4, C) with a diagonal quadratic Hamiltonian. It is known that this system always admits three functionally independent polynomial first integrals. We prove that if the system has a rational first integral functionally independent of the known three ones (so called fourth integral), then it has a polynomial fourth first integral. This is a consequence of a more general fact that for these systems the existence of a Darboux polynomial with non vanishing cofactor implies the existence of a polynomial fourth integral.

Introduction

For a given system of (polynomial) ordinary differential equations depending on parameters, a question is how to identify those values of the parameters for which the equations have (rational or polynomial) first integrals? Except for some simple cases, this problem is very hard and there are no general methods to solve it.

In this paper we present a partial result concerning the integrability problem for the so-called Euler equations on Lie algebras [1-4, 13, 14], for which the problem is also largely open.

Let us recall some relevant definition. Let (L, [•, •]) be a finite dimensional (real or complex) Lie algebra and L * its dual. For f, g ∈ C ∞ (L * ) their Lie-Poisson bracket {f, g} is defined by {f, g}(x) = x, [df (x), dg(x)] ,

where x ∈ L * , df (x), dg(x) ∈ (L * ) * = L * * ≈ L, and for x ∈ L * and y ∈ L, we denote x, y = x(y).

Let {e * 1 , . . . , e * n } be the dual basis to a fixed basis {e 1 , . . . , e n } of L. An element x ∈ L * can be written as x = n i=1 x i e * i , where the coefficients x i (real or complex) depend smoothly on x thus defining the functions ∈ C ∞ (L * ), for 1 ≤ i ≤ n.

For a given function H ∈ C ∞ (L * ), the system of differential equations

dx i dt = {x i , H}, 1 ≤ i ≤ n, (1.1)
is called Euler equations on the Lie algebra L with the Hamiltonian H.

Recall that a function F ∈ C ∞ (L * ) is a Casimir function of the Lie algebra L if {f, F } = 0, for every f ∈ C ∞ (L * ). It is clear (e.g., [START_REF] Olver | Applications of Lie groups to differential equations[END_REF]) that a function F defined on L * is a first integral of system (1.1) if and only if {F, H} = 0. In particular, the Hamiltonian H and any Casimir function of the Lie algebra L are first integrals of system (1.1).

Only for Hamiltonians H that are functionally independent of the Casimir functions, the right sides of system (1.1) does not vanish identically. That is why we will always suppose that the Hamiltonian H under considerations is functionally independent of the Casimir functions.

From now on we will concentrate only on the complex six dimensional Lie algebra so(4, C), the Lie algebra of the complex Lie group SO(4,C), and study one of the simplest Euler equations on it, namely the Euler equations corresponding to the so called diagonal quadratic Hamiltonian.

In an appropriate basis of the Lie algebra so(4, C) (see [START_REF] Adler | The algebraic integrability of geodesic flow on SO(4)[END_REF]), the Euler equations corresponding to the diagonal quadratic Hamiltonian 1 take the following elegant form:

dx 1 dt = (λ 3 -λ 2 )x 2 x 3 + (λ 6 -λ 5 )x 5 x 6 , dx 2 dt = (λ 1 -λ 3 )x 1 x 3 + (λ 4 -λ 6 )x 4 x 6 , dx 3 dt = (λ 2 -λ 1 )x 1 x 2 + (λ 5 -λ 4 )x 4 x 5 , dx 4 dt = (λ 3 -λ 5 )x 3 x 5 + (λ 6 -λ 2 )x 2 x 6 , dx 5 dt = (λ 4 -λ 3 )x 3 x 4 + (λ 1 -λ 6 )x 1 x 6 , dx 6 dt = (λ 2 -λ 4 )x 2 x 4 + (λ 5 -λ 1 )x 1 x 5 , (1.2)
where λ := (λ 1 , . . . , λ 6 ) ∈ C 6 . Exactly the same construction holds for the Lie algebra so(4, R), where λ := (λ 1 , . . . , λ 6 ) ∈ R 6 and the form of equations (1.2) remains unchanged. The Lie algebra so(4, C) admits two functionally independent polynomial Casimir functions. Thus any system of Euler equations on it always admits three functionally independent first integrals. Indeed, the above system possesses three first integrals:

H 1 = x 1 x 4 + x 2 x 5 + x 3 x 6 , H 2 = 6 i=1 x 2 i , H 3 = 6 i=1 λ i x 2 i , (1.3)
where the first integrals H 1 and H 2 are Casimir functions of the Lie algebra so(4, C) and H 3 is the Hamiltonian. These three first integrals are functionally independent unless all the λ i , 1 ≤ i ≤ 6, are equal, in which case the right hand sides of system (1.2) vanish.

For the Lie algebra so(4, C), on the constant level manifolds of two functionally independent Casimir functions, any Euler system, at least locally, can be reduced to the standard Hamiltonian equations with two degrees of freedom (see Secs. 6.1-6.2 and Theorem 6.22 from [START_REF] Olver | Applications of Lie groups to differential equations[END_REF]). Therefore whatever the chosen notion of integrability, the system (1.2) in order to be integrable needs a supplementary first integral H 4 , functionally independent of H 1 , H 2 and H 3 , called shortly a fourth integral. The only known cases when a fourth integral exists are the Manakov case, defined by the condition

M = λ 1 λ 4 (λ 2 + λ 5 -λ 3 -λ 6 ) + λ 2 λ 5 (λ 3 + λ 6 -λ 1 -λ 4 ) + λ 3 λ 6 (λ 1 + λ 4 -λ 2 -λ 5 ) = 0,
and the product case, defined by the conditions

λ 1 = λ 4 , λ 2 = λ 5 , λ 3 = λ 6 .
In both cases the fourth integral can be found among the polynomials of degree 2 at most (see [START_REF] Adler | The algebraic integrability of geodesic flow on SO(4)[END_REF][START_REF] Maciejewski | The Euler equations on Lie algebra so(4): an elementary approach to integrability condition[END_REF][START_REF] Popov | On rational integrability of Euler equations on Lie algebra so(4, C)[END_REF]).

In our study, we will concentrate only on the existence of a fourth rational integral. In fact, it is well known, its absence implies the absence of an algebraic fourth integral [START_REF] Kummer | On a result of Bruns[END_REF][START_REF] Whittaker | A treatise on the analytical dynamics of particles and rigid bodies[END_REF][START_REF] Yoshida | Necessary condition for the existence of algebraic first integrals. I. Kowalevski's exponents[END_REF] as well as the absence of a meromorphic fourth integral defined on some neighbourhood of 0 of C 6 [START_REF] Ziglin | Bifurcation of solutions and the nonexistence of first integrals in Hamiltonian mechanics[END_REF].

In [START_REF] Popov | On rational integrability of Euler equations on Lie algebra so(4, C)[END_REF], two of us proved the following theorem:

Theorem 1.1. If for some λ ∈ C 6 , the Euler equations (1.2) admit a rational fourth integral, then they admit a polynomial fourth integral.

Let us note that from the validity of Theorem 1.1 in complex setting, its validity in real one follows immediately. The proof of Theorem 1.1 presented in [START_REF] Popov | On rational integrability of Euler equations on Lie algebra so(4, C)[END_REF], based on Holomorphic Rectification Theorem (see Theorem 1.18 in [START_REF] Ilyashenko | Lectures on analytical differential equations[END_REF]) is elementary but quite long and quite involved. The aim of this note is to present a simpler proof based on the more powerful Holomorphic Frobenius Integrability Theorem (see Theorem 2.9 in [START_REF] Ilyashenko | Lectures on analytical differential equations[END_REF]). The present paper is self-contained and independent of [START_REF] Popov | On rational integrability of Euler equations on Lie algebra so(4, C)[END_REF]. To put it in such a form some overlaps with [START_REF] Popov | On rational integrability of Euler equations on Lie algebra so(4, C)[END_REF] were unavoidable.

The paper is organized as follows. In Sec. 2 we collect various facts needed for the proof. In Sec. 3, Theorem 1.1 is obtained with the help of more general Theorem 3.1. Actually, the proof of Theorem 1.1 is based on the study of so called Darboux polynomials (see Sec. 2.1) for the Euler equations (1.2) and the rich symmetry properties of these equations and Theorem 1.1 is a direct consequence of Theorem 3.1 concerning Darboux polynomials. Let us stress that all proofs are completely elementary.

Finally let us note that in [START_REF] Maciejewski | Darboux polynomials and first integrals of natural polynomial Hamiltonian systems[END_REF] exact counterparts of Theorems 1.1 and 3.1 are proved for so called natural polynomial hamiltonian systems of an arbitrary degree of freedom.

Preliminaries

Darboux polynomials

Consider a polynomial system of ordinary differential equations defined in

C n dx j dt = V j (x 1 , . . . , x n ), 1 ≤ j ≤ n. (2.1)
For a holomorphic function F defined on some open subset of C n , let us define

d(F ) = n i=1 ∂F ∂x i V i .
The operator d is called the derivation associated with system of differential equations (2.1).

A polynomial P ∈ C[x 1 , . . . , x n ] \ C is called a Darboux polynomial of system (2.1) if for some polynomial S ∈ C[x 1 , . . . , x n ] we have

d(P ) = SP.
(2.

2)

The polynomial S is called a cofactor of the Darboux polynomial P . When S = 0, P is called a proper Darboux polynomial. When S = 0, P is nothing but a first integral of system (2.1).

Here we mention some properties of the Darboux polynomials that we will need throughout:

(D1) Let P 1 and P 2 be non-zero relatively prime polynomials that are not first integrals of system (2.1). Then the rational function P 1 /P 2 is a first integral of system (2.1) if and only if P 1 and P 2 are its proper Darboux polynomials with the same cofactor. (D2) All factors of a Darboux polynomial of system (2.1) are also its Darboux polynomials. (D3) If P 1 and P 2 are two Darboux polynomials of system (2.1) with cofactors S 1 and S 2 , respectively, then P 1 P 2 is also its Darboux polynomial with cofactor S 1 + S 2 . (D4) Let us suppose that the right-hand sides of system (2.1) are homogeneous polynomials of the same degree. Let P be a Darboux polynomial of system (2.1). Then its cofactor S is homogeneous and all homogeneous components of P are also Darboux polynomials of system (2.1).

See [START_REF] Moulin-Ollagnier | On the non-existence of constants of derivations: the proof of a theorem of Jouanolou and its development[END_REF] for more details.

Permutational symmetries

The Euler equations (1.2) possess an invariance property, called permutational symmetry. In general, permutational symmetries can be described as follows. Let

x = (x 1 , . . . , x n ) ∈ C n , λ = (λ 1 , . . . , λ n ) ∈ C n , and let V (x, λ) = (V 1 (x, λ), . . . , V n (x, λ)) depends holomorphically on (x, λ) ∈ C 2n .
Let us consider the following system of differential equations

dx dt = V (x, λ). (2.3)
Let σ be an element of the symmetric group S n , i.e., the group of all permutations of {1, . . . , n}. For a = (a 1 , . . . , a n ) ∈ C n , we will denote σ(a) = (a σ(1) , . . . , a σ(n) ).

A permutation σ ∈ S n will be called a permutational symmetry of system (2.3) if for all (x, λ) ∈ C 2n , we have

V k (σ(x), σ(λ)) = εV σ(k) (x, λ), 1 ≤ k ≤ n,
where ε = ±1 is a constant independent of k. Clearly, all permutational symmetries of system (2.3) form a group.

Theorem 2.1. Let σ be a permutational symmetry of system (2.3). (a) If F = F (x) is a first integral of system (2.3), then the function

F = F • σ -1 is a first integral of the system dx dt = V (x, σ(λ)). (2.4) (b) If P = P (x) is a Darboux polynomial of system (2.3) (see (2.2)), then d( P ) = S P ,
where

P = P • σ -1 , S = S • σ -1
, and d is the derivation associated with system (2.4).

For the proof of (a) see Sec. II of [START_REF] Maciejewski | The Euler equations on Lie algebra so(4): an elementary approach to integrability condition[END_REF]. The proof of (b) follows the same line.

The group of permutational symmetries of the Euler equations (1.2) consists of 24 elements. Among others, it contains the following five permutations:

τ 2 (1, 2, 3, 4, 5, 6) = (2, 1, 3, 5, 4, 6), τ 3 (1, 2, 3, 4, 5, 6) = (3, 2, 1, 6, 5, 4), τ 4 (1, 2, 3, 4, 5, 6) = (4, 2, 6, 1, 5, 3), τ 5 (1, 2, 3, 4, 5, 6) = (5, 4, 3, 2, 1, 6), τ 6 (1, 2, 3, 4, 5, 6) = (6, 2, 4, 3, 5, 1). (2.5)
For more details see Sec. II of [START_REF] Maciejewski | The Euler equations on Lie algebra so(4): an elementary approach to integrability condition[END_REF] where, in its notations,

τ 2 = σ 1 , τ 3 = σ 3 , τ 4 = σ 7 , τ 5 = σ 8 • σ 1 and τ 6 = σ 7 • σ 3 .
Let P be a proper Darboux polynomial of system (1.2), that is d(P ) = SP , where d is the corresponding derivation and S ∈ C[x 1 , . . . , x 6 ]\{0}. Since the right-hand sides of system (1.2) are homogeneous of the same degree, it follows from (D4) that the cofactor is a homogeneous linear form, i.e.,

S = 6 i=1 α i x 1 ,
where α i ∈ C, 1 ≤ i ≤ 6, are some complex constants and at least one of them is non-zero, say α i 0 = 0.

According to (2.5) τ i 0 (i 0 ) = 1. Now, Theorem 2.1b implies that without any loss of generality, we can always assume that α 1 = 0. This fact will be used in the proof of Theorem 1.1.

From now on, d will always denote the derivation associated with the Euler equations (1.2).

Another invariance property

Beside permutational symmetries, the Euler equations (1.2) possess also another invariance property related to the change of signs of the couples of variables (x 1 , x 4 ), (x 2 , x 5 ), and (x 3 , x 6 ) respectively. More precisely, let us denote:

τ 14 (x 1 , x 2 , x 3 , x 4 , x 5 , x 6 ) = (-x 1 , x 2 , x 3 , -x 4 , x 5 , x 6 ), τ 25 (x 1 , x 2 , x 3 , x 4 , x 5 , x 6 ) = (x 1 , -x 2 , x 3 , x 4 , -x 5 , x 6 ), τ 36 (x 1 , x 2 , x 3 , x 4 , x 5 , x 6 ) = (x 1 , x 2 , -x 3 , x 4 , x 5 , -x 6 ). (2.6)
It is easy to see that for (ij) = ( 14), (ij) = (25), and (ij) = (36),

τ ij -1 • d • τ ij = -d,
which means that under any of the transformations of (2.6), the right side of equations (1.2) changes the sign.

For the polynomial

T ∈ C[x 1 , . . . , x 6 ], let us denote T (ij) := T • τ ij . Thus if T is a first integral of the system (1.
2), then T (14) , T (25) , and T (36) also are first integrals of this system.

Moreover, if P is its Darboux polynomial, that is, d(P ) = SP , then d(P (ij) ) = -S (ij) P (ij) . In particular, if (2.9)

d (P )(x) = (α 1 x 1 + α 2 x 2 + α 3 x 3 + α 4 x 4 + α 5 x 5 + α 6 x 6 )P (x), (2.7) then d (P (14) )(x) = (α 1 x 1 -α 2 x 2 -α 3 x 3 +α 4 x 4 -α 5 x 5 -α 6 x 6 )P (14) (x) (2.

Explicit form of some linear differential operators

On C 6 , equipped with coordinates (x 1 , . . . , x 6 ), we consider the linear differential operator X ij , for 1 ≤ i < j ≤ 6, defined by the formula

X ij (G) = det ∂(H 1 , H 2 , H 3 , G) ∂(x 1 , . . . , xi , . . . , xj , . . . , x 6 )
where G is a holomorphic function on C 6 , the functions H 1 , H 2 , H 3 are given by (1.3), and xr means the absence of x r .

Two such operators X 25 and X 36 will play a crucial role in the proof of Theorem 1.1 and we need an explicit formula of them. To simplify the notations, we write: λ ij = λ i -λ j for i = j, 1 ≤ i, j ≤ 6. A direct calculation gives the formulae:

X 25 = λ 63 x 1 x 3 x 6 + λ 34 x 2 3 x 4 + λ 46 x 4 x 2 6 ∂ ∂x 1 + λ 16 x 2 1 x 6 + λ 41 x 1 x 3 x 4 + λ 64 x 2 4 x 6 ∂ ∂x 3 + λ 13 x 1 x 2 3 + λ 61 x 1 x 2 6 + λ 36 x 3 x 4 x 6 ∂ ∂x 4 + λ 31 x 2 1 x 3 + λ 14 x 1 x 4 x 6 + λ 43 x 3 x 2 4 ∂ ∂x 6 , X 36 = λ 52 x 1 x 2 x 5 + λ 24 x 2 2 x 4 + λ 45 x 4 x 2 5 ∂ ∂x 1 + λ 15 x 2 1 x 5 + λ 41 x 1 x 2 x 4 + λ 54 x 2 4 x 5 ∂ ∂x 2 + λ 12 x 1 x 2 2 + λ 51 x 1 x 2 5 + λ 25 x 2 x 4 x 5 ∂ ∂x 4 + λ 21 x 2 1 x 2 + λ 14 x 1 x 4 x 5 + λ 42 x 2 x 2 4 ∂ ∂x 5 .
It is easy to see that outside of some very special subcases of the Manakov case, these two differential operators are not identically zero. Note that

X 25 (H r ) = X 36 (H r ) = 0, 1 ≤ r ≤ 3.

Linear partial differential equations

Let us consider the following system of k, where

1 ≤ k ≤ n -1, linear partial differential equations n i=1 a ij (x) ∂F ∂x i = 0, 1 ≤ j ≤ k, (2.10)
where a ij , 1 ≤ i ≤ n, are holomorphic functions defined on some open subset U ⊂ C n .

Theorem 2.2. Let x 0 ∈ U be such that rankA(x 0 ) = k, where the matrix A(x) = (a ij (x)) 1≤i≤n,1≤j≤k . Let us suppose that F and F 1 , . . . , F n-k are holomorphic on U solutions of system (2.10) such that F 1 , . . . , F n-k are functionally independent at x 0 , that means that the vectors (grad F r )(x 0 ), 1 ≤ r ≤ n -k, are linearly independent. Then there exists a neighbourhood V of the point (F 1 (x 0 ), . . . , F n-k (x 0 )) and a holomorphic function Ω defined on V, such that for any x in some neighbourhood of x 0 we have

F (x) = Ω(F 1 (x), . . . , F n-k (x)). (2.11)
This result is a direct consequence of Holomorphic Frobenius Integrability Theorem (see Sec. 2.11.20 in [START_REF] Narasimhan | Analysis on real and complex manifolds[END_REF], Theorem 2.9 in [START_REF] Ilyashenko | Lectures on analytical differential equations[END_REF] and Appendix 8 in [START_REF] Kobayashi | Foundations of differential geometry[END_REF]).

Throughout, U denotes a subset of C 6 defined by the condition that for all 1 ≤ i < j ≤ 6, and any point z ∈ U, the vectors (grad H 1 )(z), (grad H 2 )(z), (grad H 3 )(z), (grad x i )(z), (grad x j )(z) are linearly independent. Unless all λ i , 1 ≤ i ≤ 6, are equal, U is always an open dense subset of C 6 . Below, when saying that identity (2.11) is locally fulfilled, we will understand that this is so on a neighbourhood of some point from U.

Proof of Theorem 1.1.

Let us suppose that the irreducible rational fraction P 1 /P 2 , where P 1 , P 2 ∈ C[x 1 , . . . , x 6 ], is a first integral of system (1.2) and that P 1 (and thus also P 2 ) is not its first integral. Then (D1) from Sec. 2.1 implies that P 1 and P 2 are proper Darboux polynomials of system (1.2). Since the right-hand sides of system (1.2) are homogeneous of the same degree, it follows from (D2) and (D4) that system (1.2) admits also an irreducible homogeneous proper Darboux polynomial P and its cofactor is a homogeneous linear form, i.e.,

S = 6 i=1 α i x 1 ,
where α i , 1 ≤ i ≤ 6, are some constants. Since S = 0, at least one of its coefficients is not zero. As explained in Sec. 2.2, without any loss of generality we can assume that α 1 = 0.

Theorem 1.1 is now a direct consequence of Theorem 3.1. If for some λ ∈ C 6 , the Euler equations (1.2) have a proper Darboux polynomial, then they have a polynomial fourth integral.

Proof. Let P be a proper Darboux polynomial of the Euler equations (1.2). Without any lost of generality we can suppose that P is an irreducible and homogeneous polynomial. The proof is naturally divided into three almost independent parts. Part 1. Construction of a polynomial first integral. From (2.7) and (2.8) it immediately follows that R = P P ( 14) is a Darboux polynomial of system (1.2) with cofactor 2(α 1 x 1 + α 4 x 4 ), i.e.,

d(R)(x) = 2(α 1 x 1 + α 4 x 4 )R(x), (3.1)
Thus from (2.9) one deduces that the polynomial U = R (25) satisfies

d(U )(x) = -2(α 1 x 1 + α 4 x 4 )U (x),
and finally (see (D3) in Sec. 2.1) that

d(V ) = 0, X 25 (R) = det ∂(H 1 , H 2 , H 3 , R) ∂(x 1 , x 3 , x 4 , x 6 ) = 0. (3.5)
Let us note that not only U = R • τ 25 but also U = R • τ 36 . This is so because R is a homogeneous polynomial of even degree and contains monomials that have only an even sum of powers of x 1 and x 4 . Thus the monomials of R containing an even sum of powers of x 2 and x 5 contain also an even sum of powers of x 3 and x 6 and respectively, the monomials of R containing an odd sum of powers of x 2 and x 5 contain an odd sum of powers of x 3 and x 6 .

As U = R • τ 36 , exactly in the same way as (3.5), we prove that f 36 = 0 or, equivalently, that

X 36 (R) = det ∂(H 1 , H 2 , H 3 , R) ∂(x 1 , x 2 , x 4 , x 5 ) = 0. (3.6)
Equations (3.5) and (3.6) represent a system of linear homogeneous partial differential equations for R. This system has four solutions: H 1 , H 2 , H 3 and R. Consequently, H 1 , H 2 , H 3 , and R are also solutions of the linear homogeneous partial differential equation

Y (Φ) = 0, (3.7)
where Y = [X 25 , X 36 ] is the commutator (Lie bracket) of the vector fields X 25 and X 36 .

The three vector fields X 25 , X 36 and Y are necessarily linearly dependent on C 6 . Indeed, let us suppose that they are linearly independent. System (3.5)-(3.7) has three functionally independent solutions H 1 , H 2 and H 3 . Then, according to Theorem 2.2, the fourth solution R of this system is represented locally as a function of H 1 , H 2 , and H 3 , i.e., R is a first integral of system (1.2). But this is a contradiction because R is a proper Darboux polynomial of (1.2).

We compute the vector field Y (using Maple) and obtain

Y = λ 25 x 2 x 5 (λ 63 x 1 x 3 x 6 + λ 34 x 2 3 x 4 + λ 46 x 4 x 2 6 ) + λ 63 x 3 x 6 (λ 52 x 1 x 2 x 5 + λ 24 x 2 2 x 4 + λ 45 x 4 x 2 5 ) + (λ 42 x 2 2 + λ 54 x 2 5 )(λ 13 x 1 x 2 3 + λ 61 x 1 x 2 6 + λ 36 x 3 x 4 x 6 ) + (λ 34 x 2 3 + λ 46 x 2 6 )(λ 12 x 1 x 2 2 + λ 51 x 1 x 2 5 + λ 25 x 2 x 4 x 5 ) ∂ ∂x 1 + (λ 14 x 1 x 2 + 2λ 45 x 4 x 5 )(λ 13 x 1 x 2 3 + λ 61 x 1 x 2 6 + λ 36 x 3 x 4 x 6 ) + (2λ 51 x 1 x 5 + λ 14 x 2 x 4 )(λ 63 x 1 x 3 x 6 + λ 34 x 2 3 x 4 + λ 46 x 4 x 2 6 ) ∂ ∂x 2 + (λ 41 x 1 x 3 + 2λ 64 x 4 x 6 )(λ 12 x 1 x 2 2 + λ 51 x 1 x 2 5 + λ 25 x 2 x 4 x 5 ) + (2λ 16 x 1 x 6 + λ 41 x 3 x 4 )(λ 52 x 1 x 2 x 5 + λ 24 x 2 2 x 4 + λ 45 x 4 x 2 5 ) ∂ ∂x 3 + λ 52 x 2 x 5 (λ 13 x 1 x 2 3 + λ 61 x 1 x 2 6 + λ 36 x 3 x 4 x 6 ) + λ 36 x 3 x 6 (λ 12 x 1 x 2 2 + λ 51 x 1 x 2 5 + λ 25 x 2 x 4 x 5 ) + (λ 13 x 2 3 + λ 61 x 2 6 )(λ 52 x 1 x 2 x 5 + λ 24 x 2 2 x 4 + λ 45 x 4 x 2 5 ) + (λ 21 x 2 2 + λ 15 x 2 5 )(λ 63 x 1 x 3 x 6 + λ 34 x 2 3 x 4 + λ 46 x 4 x 2 6 ) ∂ ∂x 4 + (2λ 12 x 1 x 2 + λ 41 x 4 x 5 )(λ 63 x 1 x 3 x 6 + λ 34 x 2 3 x 4 + λ 46 x 4 x 2 6 ) + (λ 41 x 1 x 5 + 2λ 24 x 2 x 4 )(λ 13 x 1 x 2 3 + λ 61 x 1 x 2 6 + λ 36 x 3 x 4 x 6 ) ∂ ∂x 5 + (2λ 31 x 1 x 3 + λ 14 x 4 x 6 )(λ 52 x 1 x 2 x 5 + λ 24 x 2 2 x 4 + λ 45 x 4 x 2 5 ) + (λ 14 x 1 x 6 + 2λ 43 x 3 x 4 )(λ 12 x 1 x 2 2 + λ 51 x 1 x 2 5 + λ 25 x 2 x 4 x 5 ) ∂ ∂x 6 .
Let N be the 3 × 6 matrix composed by the coefficients of vector fields X 25 , X 36 and Y . The linear dependence of those vector fields means that the condition rankN ≤ 2

should be fulfilled. We consider two cases:

1. λ 1 = λ 4 , 2. λ 1 = λ 4 .
Case 1. Let N 123 be the minor of N that consists of first, second and third column. In particular, condition (3.8) implies

δ = det N 123 ≡ 0.
The Maple computation of δ gives

δ = λ 14 δ 1 δ 2 δ 3 ,
where

δ 1 = λ 31 x 2 1 x 2 3 + λ 16 x 2 1 x 2 6 + 2λ 63 x 1 x 3 x 4 x 6 + λ 34 x 2 3 x 2 4 + λ 46 x 2 4 x 2 6 , δ 2 = λ 56 x 1 x 5 x 6 + λ 64 x 2 x 4 x 6 + λ 45 x 3 x 4 x 5 , δ 3 = λ 12 x 2 1 x 2 2 + λ 51 x 2 1 x 2 5 + 2λ 25 x 1 x 2 x 4 x 5 + λ 42 x 2 2 x 2 4 + λ 54 x 2 4 x 2 5 .
As we consider now the case λ 1 = λ 4 , then the condition δ ≡ 0 implies that either δ 1 ≡ 0 or δ 2 ≡ 0 or δ 3 ≡ 0. From δ 1 ≡ 0 we obtain

λ 1 = λ 3 = λ 4 = λ 6
that is a subcase of the Manakov case. From δ 2 ≡ 0 we obtain

λ 4 = λ 5 = λ 6
that is a subcase of the Manakov case. From δ 3 ≡ 0 we obtain

λ 1 = λ 2 = λ 4 = λ 5
that is a subcase of the Manakov case.

In this way we conclude that in Case 1 condition (3.8) is eventually fulfilled in some subcases of the Manakov case only.

Case 2. In this case, condition (3.8) is fulfilled. Equations (3.5) and (3.6) have four functionally independent solutions. We know three of them explicitly: they are H 1 , H 2 and H 3 . As it is easily seen, equations (3.5) and (3.6) have also the solution

R = λ 42 x 2 2 + λ 45 x 2 5 .
We prove that out of the Manakov case, R is functionally independent of and therefore either λ 6 = λ 4 or λ 2 = λ 4 = λ 5 . As above, both possibilities lead to subcases of the Manakov case. So far we have proved that outside of the Manakov case the functions H 1 , H 2 , H 3 and R are four functionally independent solutions of equations (3.5) and (3.6). By Theorem 2.2 we know that locally any other holomorphic solution of those two equations is a function of H 1 , H 2 , H 3 and R. Thus our Darboux polynomial R, being a solution of equations (3.5) and (3.6), is represented locally in the following way

H
R = W (H 1 , H 2 , H 3 , R),
where W is some holomorphic function.

The polynomial R satisfies equation (3.1) and therefore

d(R) = ∂W ∂R d(R) = 2(α 1 x 1 + α 4 x 4 )W,
which shows that the functions H 1 , H 2 , H 3 , R, d(R) and α 1 x 1 + α 4 x 4 are functionally dependent. Thus the determinant ∆ of their Jacobi matrix is identically zero, i.e.

∆ = 16 x 4 x 5 x 6 x 1 x 2 x 3 x 1 x 2 x 3 x 4 x 5 x 6 λ 4 x 1 λ 2 x 2 λ 3 x 3 λ 4 x 4 λ 5 x 5 λ 6 x 6 0 λ 42 x 2 0 0 λ 45 x 5 0 ∂d(R) ∂x 1 ∂d(R) ∂x 2 ∂d(R) ∂x 3 ∂d(R) ∂x 4 ∂d(R) ∂x 5 ∂d(R) ∂x 6 α 1 0 0 α 4 0 0 ≡ 0.
It is easy to compute ∆ (using Maple) but this is not necessary. As ∆ equals zero identically, we can compute it only for some fixed particular values of certain variables x i , 1 ≤ i ≤ 6. For x 1 = 0, x 5 = 0, and x 6 = 0 we obtain easily (and without any computer algebra tools) that

∆ x 1 =x 5 =x 6 =0 = 16α 1 x 2 x 3 x 2 4 λ 34 λ 42 (λ 64 λ 42 x 2 2 + λ 34 λ 45 x 2
3 ) ≡ 0. (3.9)

As α 1 = 0, identity (3.9) is fulfilled in three cases only:

A. λ 2 = λ 4 , B. λ 3 = λ 4 , C. λ 4 = λ 5 = λ 6 .
Case A. We compute ∆ for x 1 = 0, x 2 = 0, and x 3 = 0 and obtain Case C. This case is directly a subcase of the Manakov case. Thus the assumption that H 1 , H 2 , H 3 and V are functionally dependent when R and U are relatively prime can eventually be true in some very special subcases of the Manakov case only.

∆ x 1 =x 2 =x 3 =0,λ 2 =λ 4 = 16α 1 λ 2 45 λ 43 λ 64 x 2 4 x 3 5 x 6 ≡ 0. ( 3
Remark. We have to note here that there are, indeed, some subcases of the Manakov case when our procedure does not lead to a fourth integral. For example, when λ 1 = λ 4 = λ 5 = λ 6 = 0 and λ 2 = -λ 3 (subcase of the case C), the polynomial P = x 2 + x 3 is a proper Darboux polynomial of the Euler equations (1.2). Applying our procedure on P, one obtains, however, a polynomial first integral that is functionally dependent on H 3 .

But we know that in the Manakov case there always exists a polynomial fourth integral. That is why we do not exclude the Manakov case from the condition of the theorem.

Part 3. R and U are not relatively prime polynomials.

We have for R and U R = P P [START_REF] Perelomov | Integrable systems of classical mechanics and Lie algebras[END_REF] and U = P (25) P (14)(25) .

Since the polynomial P is irreducible, the polynomials P [START_REF] Perelomov | Integrable systems of classical mechanics and Lie algebras[END_REF] , P (25) and P (14)(25) are also irreducible. Thus polynomials R and U are not relatively prime in the following 8 cases only:

1 According to (2.9) the cofactor of P (25) is -α 1 x 1 + α 2 x 2 -α 3 x 3 -α 4 x 4 + α 5 x 5 -α 6 x 6 .

Since P and P (25) are equal in the case under consideration, comparing the two cofactors we find α 1 = 0, α 3 = 0, α 4 = 0, α 6 = 0.

This contradicts, however our assumption that α 1 = 0. In the same way, cases 2, 3 and 4 also lead to α 1 = 0.

As an example of an application of the constructing procedure of the fourth integral described in the above proof, let us consider the product case when λ 1 = λ 2 and λ 1 = λ 3 . One can easily see that in this case the polynomial is a fourth integral of (1.2). In fact, in this example, already P P (25) is a fourth integral. The explicit form of the polynomial fourth integral when λ 2 = λ 1 and λ 2 = λ 3 or when λ 3 = λ 1 and λ 3 = λ 2 follows now from Theorem 2.1b applied to the permutational symmetries τ = τ 2 • τ 3 and τ 2 respectively.

Remark. When comparing our system (1.2) with its "twin brother", which are the Euler-Poisson equations of heavy rigid body motion (see [START_REF] Borisov | Poisson structures and Lie algebras in hamiltonian mechanics[END_REF][START_REF] Fomenko | Integrability and nonintegrability in geometry and mechanics[END_REF][START_REF] Perelomov | Integrable systems of classical mechanics and Lie algebras[END_REF][START_REF] Trofimov | Introduction to geometry of manifolds with symmetry[END_REF][START_REF] Whittaker | A treatise on the analytical dynamics of particles and rigid bodies[END_REF]), we conclude from [START_REF]Bifurcation of solutions and the nonexistence of first integrals in Hamiltonian mechanics[END_REF] (see also [START_REF] Maciejewski | Differential Galois approach to the nonintegrability of the heavy top problem[END_REF]) that for these equations the exact counterpart of Theorem 1.1 holds. Nevertheless, the exact counterpart of Theorem 3.1 for the Euler-Poisson equations fails. Indeed, in the non-integrable so-called Hess-Appelrot case, proper Darboux polynomial exists.

  8) and d (P (25) )(x) = (-α 1 x 1 + α 2 x 2 -α 3 x 3 -α 4 x 4 + α 5 x 5 -α 6 x 6 )P (25) (x).

. 10 )x 1 =x 2 =x 3 =0,λ 3 =λ 4 = 16α 1 λ 54 λ 2 46 λ 42 x 2 4 x 5 x 3 6

 10346 As α 1 = 0 identity (3.10) is fulfilled if and only if either λ 5 = λ 4 or λ 3 = λ 4 or λ 6 = λ 4 . Each of these possibilities together with the conditions λ 1 = λ 4 and λ 2 = λ 4 leads to a subcase of the Manakov case.Case B. We compute ∆ also for x 1 = 0, x 2 = 0 and x 3 = 0 and obtain∆ ≡ 0. (3.11)As α 1 = 0, identity (3.11) is fulfilled if and only if either λ 5 = λ 4 or λ 6 = λ 4 or λ 2 = λ 4 . Each of these possibilities together with the conditions λ 1 = λ 4 and λ 3 = λ 4 leads to a subcase of the Manakov case.

P = λ 21 c x 2 + x 3 + λ 21 c x 5 + x 6 ,- λ 21 λ 13 (x 2 + x 5 ) 2 + (x 3 + x 6 ) 2 2

 2561325232 where c = √ λ 13 λ 21 , is a proper Darboux polynomial of system (1.2) with cofactor c(x 1 + x 4 ). Here P = P[START_REF] Perelomov | Integrable systems of classical mechanics and Lie algebras[END_REF] and thus R = P P (14) = P 2 and U = (P 2 ) (25) = P 2 (25) . Finally, the polynomial V = RU = P P (25) 2 =

  [START_REF] Adler | The algebraic integrability of geodesic flow on SO(4)[END_REF] , H 2 and H 3 . Indeed, if R were functionally dependent on H 1 , H 2 , and H 3 , then R would be a first integral of the Euler equations (1.2). We compute d(R) and obtaind(R) = 2(λ 24 λ 34 x 1 x 2 x 3 + λ 54 λ 64 x 1 x 5 x 6 + λ 24 λ 64 x 2 x 4 x 6 + λ 54 λ 34 x 3 x 4 x 5 ).A necessary condition for the equation d(R) ≡ 0 to hold is that eitherλ 2 = λ 4 or λ 3 = λ 4 . Let us first suppose that λ 2 = λ 4 . Then d(R) = 2λ 54 x 5 (λ 64 x 1 x 6 + λ 34 x 3 x 4 )and therefore either λ 5 = λ 4 or λ 3 = λ 4 = λ 6 . Both possibilities together with the condition of Case 2, i.e., λ 1 = λ 4 , lead to subcases of the Manakov case.Let us suppose now that λ 3 = λ 4 . Then d(R) = 2λ 64 x 6 (λ 54 x 1 x 5 + λ 24 x 2 x 4 )

  . P = P (25) ; 2. P = -P (25) ; 3. P = P (14)(25) ; 4. P = -P (14)(25) ; 5. P (14) = P (25) , that is, equivalent to case 3; 6. P (14) = -P (25) , that is, equivalent to case 4; 7. P (14) = P (14)(25) , that is, equivalent to case 1; 8. P (14) = -P (14)(25) , that is, equivalent to case 2.Let us examine case 1. The cofactor of P is α 1 x 1 + α 2 x 2 + α 3 x 3 + α 4 x 4 + α 5 x 5 + α 6 x 6 .

Acknowledgements

We would thank Andrzej J. Maciejewski (University of Zielona Góra, Poland) for interesting discussions.

where V := RU = RR (25) = (P P ( 14) )(P P ( 14) ) (25) = P P (14) P (25) P (14)(25) .

This means that V is a polynomial first integral of the Euler equations (1.2).

The main difficulty is to decide whether V is a fourth integral. We will prove that this is always the case outside of some very special subcases of the Manakov case. This is proved in Part 2 if the polynomials R and U are relatively prime and in Part 3 if this is not the case. As in the Manakov case the polynomial fourth integral always exists, this will prove Theorem 3.1.

Part 2. R and U are relatively prime polynomials.

We have to decide when the first integrals

3)) and V are functionally independent. Let us suppose that they are functionally dependent.

Then for all α i , 1 ≤ i ≤ 6

We will prove that outside of very special subcases of the Manakov case this contradicts α 1 = 0.

We have supposed that the polynomials R and U are relatively prime and thus (3.2) shows that either R divides X ij (R), i.e.,

where f ij is a homogeneous polynomial of second degree, or X ij (R) = X ij (U ) = 0. For the first possibility, according to (3.2) and (3.3), we have that

In particular X 25 (R) = f 25 R and X 25 (U ) = -f 25 U . Applying to the first identity the change of variables τ 25 (see Sec. 2.3), we conclude that X 25 (U ) = (f 25 •τ 25 )U and finally that f 25 = -f 25 •τ 25 . But this is impossible because f 25 cannot depend on x 2 and x 5 . Indeed, the maximal powers of x 2 and of x 5 in X 25 (R) are never greater than their respective maximal powers in R. Thus f 25 = 0 and consequently X 25 (R) = X 25 (U ) = 0.

We have thus proved that R satisfies the equation