Ralf Winkler

Susanne Fuchs

Pascal Perrier
email: pascal.perrier@gipsa-lab.grenoble-inp.fr

Mark Tiede
email: tiede@haskins.yale.edu

Biomechanical Tongue Models: An Approach to Studying Inter-speaker Variability

Keywords: Biomechanical tongue models, speaker-specific articulation, variability

Speakers of a given language vary with respect to their acoustics, articulation, and motor commands. This variation is driven by a variety of influences, such as emotional states, communicative interaction, and individual properties of the vocal tract. In this work we focus on the latter. First, we build speaker-specific biomechanical tongue models. Second, we discuss the impact of the relative position of the bending in the vocal tract on the basis of extensive simulations with two different models. We focus on /i,a,u/ by defining target regions in the acoustic space, and discuss the corresponding speakerspecific articulatory and motor command variability observed.

Introduction

Similar to ideas recently developed in the domain of embodied cognition we suppose that the realization of abstract categories like phonemes is constrained by the physical properties of the speech production apparatus. In order to understand the realization of phonemes, we try to understand the functioning of speech production mechanisms and their corresponding patterns of speech motor control. To do so, over the last decade we have carried out a number of kinematic studies focusing on control strategies in the production of different phonemes, the interaction of the tongue with its enclosing vocal tract boundaries, and the impact of speech rate on the temporal organization of inter-articulatory coordination. Moreover, we have compared the experimental data with simulations of a physically realistic tongue model in order to test our predictions concerning the underlying speech motor control mechanisms. In the current step of our work, we have built speaker-specific biomechanical tongue models to explain speaker-specific articulatory behavior on the basis of differences in vocal tract shape. We consider physically realistic models to be primarily biomechanical, although geometrical models for instance can also involve a physically realistic description of the tongue and the surrounding vocal tract in the mid-sagittal plane. However, what is different in biomechanical models is that they include a description of the macro muscle fibers which are recruited to produce a certain speech movement. One might critically ask the questions: What kind of additional information can biomechanical models provide in comparison to geometrical models? Why is it worth using biomechanical models given their high computational cost? First of all, we think the choice of the model crucially depends on the task to be fulfilled and one should choose a model appropriate to the research question. Biomechanical tongue models have the advantage that muscular activity, muscular forces, stiffness and muscular co-contraction underlying a given speech movement can be studied. Hence, they provide a window into the motor command level and speech motor control. Geometrical models may also follow biomechanical principles if they are based on large amounts of experimental data, but the movement between different tongue positions usually consists of a linear interpolation and does not follow physical principles. To study speaker-specific behavior, as is the aim of this paper, it is possible to build either speakerspecific geometrical or biomechanical models [START_REF] Fuchs | Benoit Project: Speaker specific vowel articulation[END_REF]. The crucial test, however, is to map the parameter space of one speaker onto another. In geometrical models the physical realism of this mapping is unknown. We think that for such a task biomechanical models provide a more elegant solution since macro fibers are mapped onto the tongue surface/volume and all articulatory movements and degrees of freedom are a result of muscle activation, not on the articulatory parameter space of some reference speaker. This reason and the fact that biomechanical models support the investigation of the relations between variability of the motor commands, of the articulation and of the acoustics motivated us to build these models.

Modeling methodology

The following processing steps describe how speaker-specific models are constructed, and are grouped by the overall aims of modeling speaker-specific articulation and speaker-specific acoustics. This section ends with an overview of how all the components act together.

General tasks

Image data acquisition

Speakers were recorded by means of Magnetic Resonance Imaging (MRI) with a Philips Gyroscan T10-NT Powertrack 1000 scanner. The image matrix was composed of 256 x 256 pixels, each with a spatial resolution of 1 mm in the ydirection and 1.4 mm in the x-direction. Data were originally collected for 10 isolated vowels to study inter-speaker acoustic and articulatory variability [START_REF] Apostol | Étude et simulation des caractéristiques individuelles des locuteurs par modélisation du processus de production de la parole[END_REF]. For each vowel three 18 slice series of 3.6 mm thick parallel sections with a 0.4 mm distance between slices were recorded. The three slice series differ in their respective slice orientation (transverse, coronal, oblique), chosen to be approximately orthogonal to the tract airway over some portion of its length.

Airway segmentation

Within each MR image the airway was segmented from the surrounding tissues manually. We used the itk-SNAP (version 2.1.4) software [START_REF] Yushkevich | User-guided 3D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability[END_REF] for segmentation. Several rules were established to ensure reliable segmentation suitable for speaker-specific articulatory and acoustic modeling. The biomechanical tongue model already includes standard teeth and standard lips. For that reason segmentation of the air channel terminates at the incisors. For acoustic purposes, the front tube for the lip region used during area function calculation is invariant. Another convention was to always exclude the epiglottis from segmenting the tongue, because the tongue contour is handled separately in the articulatory model. In cases where the epiglottis was located in the middle of the airway, the tissue of the epiglottis was ignored and the contour of the tongue body and the pharyngeal walls determined the cross-sectional area. If the epiglottis touched the tongue body, and hence tissues were hard to separate, tissue of the epiglottis was ignored on the basis of a smooth contour of the tongue back. For acoustic modeling the area function in the region of the epiglottis was systematically decreased (see Sec. 2.3.3). A last convention regards the uvula, which is not part of the biomechanical tongue model. In order to allow for comparable distance function values in the corresponding region, the uvula was excluded during image segmentation as well.

Spatial reconstruction

Airway segmentation was accomplished by producing binary images of each MR image with one value for segmenting a pixel as air surrounded by tissue, and a second value otherwise. The resulting contours were sub-sampled with 100 equally spaced points. The single contours were combined to form a vocal tract by re-orientating each contour of the three series of 18 slices according to their original spatial orientation. The derived 3D wire-frame model of the relevant vocal tract airway constitutes the basis for speaker-specific modeling.

Biomechanical modeling

The biomechanical component of the model is based on the 2D tongue model [START_REF] Perrier | Influences of tongue biomechanics on speech movements during the production of velar stop consonants: A modeling study[END_REF]. This model has been shown capable of accounting for some important kinematic characteristics of speech articulation such as velocity profiles, trajectory shapes [START_REF] Perrier | Influences of tongue biomechanics on speech movements during the production of velar stop consonants: A modeling study[END_REF], or relations between trajectory curvature and speed [START_REF] Perrier | Speed-curvature relations in speech production challenge the one-third power law[END_REF]. It mainly consists of a deformable Finite Element Mesh (FEM) embedded in rigid vocal tract walls in the mid-sagittal plane. The geometry of the mesh has been specifically designed to facilitate anatomical implementation of the muscles within the tongue. The model is controlled according to the -model [START_REF] Feldman | Once More on the Equilibrium-Point Hypothesis (Lambda-Model) for Motor Control[END_REF], which generates muscle force as a function of the difference between a centrally specified threshold length and the actual muscle length. Two basic hypotheses underlie the design of the speaker-specific models: (1) the general anatomical arrangements accounted for by the mesh geometry is common to all human beings, (2) variations across speakers in muscle lengths and muscle orientations within the tongue are the results of global variations of the head morphology such as variations in larynx height, mandibule ramus length, head size, and mid-sagittal palate shape. Hence anatomical landmarks have been selected that enable characterization of each speaker using these morphological parameters: the mid-sagittal palatal contours and tongue contours at rest, the lower and upper limits of the tongue insertions on the mental spine (P1 and P2), and the styloid process (P3) were the landmarks chosen.

Creating the mid-sagittal vocal tract contour

The individual mid-sagittal palatal and tongue contours were extracted from MRI data for the subjects at tongue rest position and used to align the coronal and transverse contours of the vocal tract resulting from the airway segmentation (Section 2.1.2). After determining a global left and a global right extremum of the vocal tract, the mid-sagittal plane was established as the vertical plane passing through the center of these extrema. Then, adjustment of the mid-sagittal position was done by hand for each speaker interactively to ensure that the plane bisected all the contours.

Determination of anatomical landmarks

P1, P2 and P3 were measured from a high-definition midsagittal view of the speaker's head recorded with anatomical MRI. P1 and P2 were determined on the basis of grey level changes in the mental spine region. P3 was determined using the internal contour of the sphenoid bone in the mid-sagittal plane. It was located at 1/3 of the global length of this contour from the bottom.

Model creation and adaptation

Given these anatomical landmarks matching the original generic biomechanical model to speaker-specific anatomy was straightforward. First the upper contour of the tongue model was projected onto the mid-sagittal tongue contour measured for the subject. This contour corresponds to the upper limit of the new speaker-specific FEM tongue mesh. Second, the distribution of the nodes along this new upper contour was done proportionally to the distribution of the nodes in the original model, selected for anatomical appropriateness. Third, the lower and upper attachment points of the new tongue mesh on the mandible were assigned to points P1 and P2. Then, the distribution of the nodes within the mesh was obtained by deforming the original mesh linearly from the nodes on the upper contour to the insertion nodes of the mesh into the mental spine. Finally, the extremity of the external Styloglossus fiber was attached to point P3. This matching procedure fully determines the geometry of the new mesh and consequently muscle arrangement within the new speakerspecific tongue model. It preserves the original topology of the mesh while accounting for the speaker-specific morphology.

Acoustic modeling

In order to determine acoustics from articulation, the 2D distance function resulting from the biomechanical tongue model has to be converted to its corresponding area function. In our approach this is accomplished by determining the distance function based on a pre-defined grid and subsequent reconstruction of the area according to the coefficients associated to the respective grid line.

Grid orientation

Different kinds of grids have been suggested for reconstructing area functions from vocal tract distance functions, following at least two constraints which are important for determining acoustics: grid lines should ideally be perpendicular to the vocal tract midline and area values should not change abruptly. In our approach we use the grid proposed in [START_REF] Perrier | Influences of tongue biomechanics on speech movements during the production of velar stop consonants: A modeling study[END_REF] and adapt it to the vocal tract morphology of a particular speaker. The applied grid additionally allows the precise determination of coefficients necessary for the calculation of the area function, since the grid line orientation roughly follows the orientation of the slices during MRI acquisition.

Determining α coefficients

For the purpose of area function calculation we adapted a modified version of the α-β-model proposed in [START_REF] Perrier | Vocal tract area function estimation from midsagittal dimensions with CT scans and a vocal tract cast: Modeling the transition with two sets of coefficients[END_REF]. Essentially α-coefficients have different values depending on the vocal tract region (i.e. lips, palate, etc.), but α values also vary within one region for different speakers. For that reason we determine α-values for each speaker as follows: 3D wireframe representations of the three cardinal vowels as well as the neutral vowel are used to relate the distance values to the area values of a speaker. We further determine adequate thresholds corresponding to small and large distance values, leading to two speaker-specific αvalues per contour.

Determining area functions and acoustics

Area values are calculated based on the mid-sagittal distances by applying the corresponding α coefficient of large values, if a distance is larger than the pre-determined threshold for large sagittal distances. If a distance is less than the pre-determined threshold for small distances, α-coefficient of small values were applied. For distances in between, α-coefficients were linearly interpolated before calculation of an area value. The segmentation of the airway in the image volumes ends in the region of the incisors. To control for variation in the lip region the area values at the lips are modeled by a tube of length 10 mm and 19.5 mm diameter for open vowels and 20 mm and 8 mm respectively for closed vowels. Determining acoustics from area functions does not involve any speaker-specific adaptation. In our approach formants are simply computed by coupling an acoustic analog of the vocal tract with the reconstructed area functions.

Simulating tongue movements

In order to simulate a tongue movement, motor commands are interpolated starting at values resulting from tongue rest position and reaching the specified target motor command. Within this modeling framework we are capable of relating every tongue target position represented as set of six muscle activations to the resulting tongue shape and the corresponding first two formant values.

Studying speaker-specific variability

In the following section we present a first application of the individual biomechanical models. Two speaker-specific models were used to investigate the impact of the location of the bending vocal tract on the inter-speaker variability in the articulatory and motor command space.

Motivation

Acoustic models of speech production suppose that the length of the vocal tract, and the length and location of the constriction defining the resonance cavities are crucial parameters for the description of the spectral properties of vowels. The location of the vocal tract bend seems to play a negligible role with respect to the acoustics, assuming the length of the tract is kept constant [START_REF] Sondhi | Resonances of a bent vocal tract[END_REF]. However, the relative location of the bend within the vocal tract is crucial with respect to determining the vertical and horizontal dimensions of the tract and is likely to influence the degrees of freedom of articulatory motion in the respective directions. Moreover, the location of the bend probably also affects the biomechanical properties of the tongue, since the tongue is soft tissue and deforms to its surrounding structures. Hence, the location of the vocal tract bend may be one source of inter-speaker articulatory variation and may also be reflected in the motor command level. We hypothesize that the location of the vocal tract bend and the corresponding relation between horizontal and vertical portions of the vocal tract affect the degrees of freedom in the articulatory domain. For instance, a speaker with a relatively long vertical vocal tract dimension and a short horizontal dimension (i.e. a rather anterior location of the bending vocal tract) should have a relatively large degree of freedom moving vertically up and down, but should be constrained in the horizontal direction, moving the tongue front and back.

Method

In order to study inter-speaker variability two individual biomechanical tongue models were produced. The two distinctive speakers used for building the models were selected from a MRI database of 9 French speakers. One female speaker (AV) has a vocal tract with approximately equal lengths of the pharynx (vertical dimension) and the palate (horizontal dimension from the velum to the incisors). The male speaker (CS) has, in comparison to AV, a longer pharynx and a slightly shorter palate. Thus, the two speakers represent vocal tract morphologies typically found in female and male speakers, and they show some differences in the relative location of the vocal tract bend. They also show some differences in the shape and steepness of the palate. These differences in vocal tract shape have already been quantified by means of a Principal Component Analysis in [START_REF] Fuchs | Do Speakers' Vocal Tract Geometries Shape their Articulatory Vowel Space?[END_REF]. For each model we ran a series of 8500 simulations in order to account for the complete vowel space. To run the simulations, motor commands of six major tongue muscles were randomly varied with respect to the values at rest position. Assuming that auditory targets are the primary goals for vowel production, we kept the acoustic variability across speakers constant. Three formant ellipses were defined in the F1/F2 plane for the three corner vowels. Prior to the definition of the formant ellipses, the area functions were manipulated to match the acoustic vowel space of the two speakers. In a first step the length of the epilaryngeal tube of speaker CS was shortened to match that of speaker AV. Secondly, the original area values in the laryngeal region were fixed to a physiologically realistic value of 0.4 cm 2 . In the region of the epiglottis the area values were decreased by 50% for both speakers to avoid an artificially large cavity resulting from the excised epiglottis in the tongue model. Finally, the area functions were re-scaled to a vocal tract length of 17 cm, since we wanted to disentangle the impact of vocal tract length and the location of vocal tract bending. For each acoustic vowel target 60 simulations per speaker were randomly selected whose respective first two formant values are equally distributed over the corresponding ellipse area. To summarize, we define the same acoustic target regions for the three corner vowels and investigate the corresponding speaker-specific articulatory variability underlying these acoustic goal regions. Moreover, we obtain some indication of the variability of the motor commands.

Results

Figure 1 shows the articulatory results for model CS and model AV. Dispersion ellipses are plotted at the constriction location for /i/, for /u/ and for /a/. For /i/ the main diagonal of the dispersion ellipses follows the palatal contour in both models. The shape of the palate, which differs among the two models, is therefore crucial for the production of /i/. Model AV shows greater variability in the horizontal dimension (as would be predicted) for the most posterior dispersion ellipse, since this is at a location where the palate starts to be straight. The main difference in terms of motor commands between both models is the variability in Styloglossus activation. Model AV shows large variation of Styloglossus activity in comparison to model CS. For /u/ and /a/ the articulatory results confirm our hypothesis. Model AV always shows a larger variability in the horizontal direction than model CS, but it is more constrained in the vertical direction than model CS. In /u/ model AV shows a greater variability in Hyoglossus activation than model CS and for /a/ the variability of Hyoglossus and Styloglossus are greater for model CS than model AV. The linkage between articulatory positions and muscle activation is not straightforward, because tongue positioning is the result of the activation of several muscles simultaneously, which makes the linkage highly nonlinear. However, certain combinations of muscle activation lead to the characteristic constrictions for vowels and so far we consider these muscles to be important for the study of variability. For /a/ a constriction is formed in the pharyngeal region. This constriction is produced by a combination of Hyoglossus (down and back) and Styloglossus (up and back) activation. Both muscle activations display a larger variability in the CS model than in the AV model, but the AV model also shows smaller λ-values, and then larger activation levels. For /i/ the palatal constriction is primarily accomplished by a combination of Posterior Genioglossus (GGP, up and front) and Styloglossus activation. Model AV allows more variability in Styloglossus activation than model CS. Note however that this is not the case for GGP. This reflects the very precise requirement for a small constriction area located in the front palatal region, where the force of GGP (front direction) is directed. For /u/ the velar constriction results primarily from the activation of Styloglossus and GGP (similar to /i/, but with a different balance), in combination with the activation of Hyoglossus, which together move the constriction location backward, as in the production of /a/. Our data do not show any difference in λ-activations across subjects. There is also a trend in both subjects for greater variability for /u/ as compared to /a/ and /i/. This may reflect the fact that the requirements for tongue control accuracy are less for /u/ than for /i/ or /a/ due to the rounded lips.

Conclusion

Speaker-specific biomechanical models have been designed to study articulatory and motor control variability. This approach has been shown to give interesting insights into individual articulatory strategies, and can be used to relate them to individual morphological characteristics. Future work will focus on interpreting this data in relation to the achievement of perceptual goals.

Figure 1 :

 1 Figure 1: Dispersion ellipses showing articulatory results for model CS on top and model AV at the bottom. Dashed tongue contour corresponds to rest position.

Acknowledgements

This work was supported by the German Research Council to the SPEECHart project (Grant Nr. FU 791/1-1) and by the French-German University to the PILIOS project. It was initiated thanks to the Christian Benoit Award given to the second author in 2007.


P Le Quéré
email: plq@limsi.fr

ON THE COMPUTATION OF SOME EXTERNAL OR PARTIALLY ENCLOSED NATURAL CONVECTION FLOWS

We address the issue of the computation of natural convection in external or partially enclosed domains, such as natural convection flow in vertical open channels or along heated vertical plates which raises the question of the proper choice of artificial boundary conditions and of their numerical implementation. To this aim we perform a singular value decomposition of the discrete Stokes operator that arises from the discretization of the governing Navier-Stokes equations, using the classical staggered grid formulation. We show that some choices of boundary conditions lead to an increase of the kernel of this operator making the solution of the full nonlinear equations indefinite. This arises in particular when Neumann type boundary conditions are imposed on the velocity component normal to the main inlet and outlet boundaries. The existence of one or more velocity-pressure combinations solutions of the homogeneous Stokes operator can be used to derive a superposition algorithm in which a linear combination of these modes is added to a particular solution of the nonlinear equations in order to satisfy given constraints such as one or more pressure conditions for instance. Sample calculations are performed to demonstrate the effectiveness of this methodology, which has many applications in other configurations such as pipe flows with several outlets or free plane jets for instance.

INTRODUCTION

Natural convection flows abound in nature or in industrial configurations (e.g [START_REF] Gebhart | Buoyancy-induced flows Hemisphere Publishing Corporation 2[END_REF]). Such flows can occur either in confined spaces or in unbounded domains, very often also in partially bounded domains, that is in domains connected to the ambient surroundings through one or several apertures. Surprisingly enough external or partially enclosed flows have received much less interest in the computational literature than fully enclosed flows. Since both types of flows are governed by the same partial derivative equations, the explanation has to be found in the difficulty of deriving and imposing appropriate boundary conditions. Whereas for forced convection the imposed flow can be imposed at the entrance boundaries provided those are set sufficiently far upstream, this procedure can no longer be applied for natural convection as the flow magnitude and direction at the inflow boundary are determined by what happens close to the heated source located downstream the entrance. This difficulty has generally led people to consider boundary conditions that exert as less constraints as possible on the flow, in the aim of letting the flow establish itself in magnitude and direction.

In particular, the computation of flows such as developing external natural convection boundary layers along a heated plate, of channel flows, or flows due to isolated heated bodies necessitates considering artificial boundaries of the computational domain located at some distance away from the heated parts.

GOVERNING EQUATIONS

We assume that the governing equations are those of an incompressible flow of a newtonian fluid under the Boussinesq approximation. These equations are made dimensionless using H as unit length, Vref= (/Η) Ra 1/2 as reference velocity where  is the thermal diffusivity, (Η 2 /) Ra -1/2 as reference time and ΔT as unit temperature difference. In dimensionless form, these equations read:

T Ra z T w x T u t T T u Ra z P z w w x w u t w u Ra x P z u w x u u t u z w x u                    2 / 1 2 / 1 2 / 1 1 Pr Pr Pr 0                          
where u and w stand for, respectively, the dimensionless horizontal and vertical component velocities, P the pressure and T the dimensionless temperature. Ra is the Rayleigh number (=gβΔTH 3 /ν), and Pr is the Prandtl number (=ν/κ) Let us assume that these equations are discretized in time with a time step Δt using an implicit discretization of the linear diffusive terms and an explicit treatment of the nonlinear terms. Thus at each time step, the velocity and pressure satisfy the following linear system of equations known as the unsteady Stokes problem:

Sw z P w H Su x P u H z w x u w u               0
where H u and H w are the elliptic operators (Δ -λ I) with λ a constant that depends on Δt and Ra. Typically λ is equal to C Ra 1/2 / Δt with C a constant that depends on the specific time stepping scheme. At steady state, the velocity and pressure satisfy the same system with λ = 0. This elliptic system has to be supplemented with boundary conditions. Let us

ANALYSIS

Continuous equations: It is well known than the corresponding Stokes problem is singular. This is due to the fact that the 3-uplet (u,w,P)=(0,0,1) satisfies the homogeneous equations whatever the imposed boundary conditions. This corresponds to the fact that in an incompressible flow, the pressure is defined up an arbitrary constant. Now the basic question is the following: do the imposed boundary conditions allow for other linearly independent modes solution of the homogeneous equations? Or, in other words, do the imposed boundary conditions increase the size of the kernel of the Stokes operator? For the channel flow with the above specified boundary conditions, the answer is obvious as is known from any undergraduate textbook : a Poiseuille flow, with parabolic w-velocity component in x and linear pressure in z, satisfies the homogeneous steady Stokes problem, making the kernel of the steady Stokes problem of dimension at least 2. This means that if a steady state solution of the inhomogeneous equations is obtained in some way, it is not unique since another solution equal to the previous one plus any Poiseuille is solution of the inhomogeneous equations. This indeterminacy probably explains some of the discrepancies that have been observed in the literature. Now the following questions arise. Are there other "boundary modes" in the kernel? What happens for the type of geometry considered by Kettleborough [START_REF] Kettleborough | Transient laminar convection between heated vertical plates including entrance effects[END_REF] for instance? What happens for flows around heated bodies? As soon as the geometry becomes less trivial, the answer has to be obtained directly from the discretized equations. Looking at this question numerically may raise new difficulties owing to the classical question of the compatibility of the velocity and pressure spaces, which may add additional spurious modes in addition to those of the continuous problem. Discretized equations: We will restrict ourselves here to the classical staggered grid discretization, which is known to provide compatible velocitypressure approximation spaces, giving rise to the only constant pressure mode in the case of Dirichlet boundary conditions. The restriction to the staggered grid discretization may seem arbitrary, but it is not our goal to study this question in general but rather to provide a general framework which can be used in other types of spatial discretization, arising from finite element or spectral type approximations. Consider the unit square Ω = [0,1]x[0,1]. Let Ω be covered with NxM cells. With the staggered grid arrangement, the component u is defined by Nu = (N+1)xM degrees of freedom, the component w is defined by Nw = Nx(M+1) degrees of freedom, and the pressure by Np = NxM. The Stokes operator is thus of order (Nu+Nw+Np). The boundary conditions are imposed in the appropriate way. In order to explicitly build the operator, all fields are considered as 1D vectors resulting from the natural ordering. The explicit construction of the Stokes operator allows one to perform its Singular Value Decomposition. The SVD allows one to fully characterize the matrix, as it provides the dimension of the kernel of the operator, and orthonormal basis of its kernel and of its image as well as those of its transpose. More precisely any matrix S, here square, can be written as S = U Λ V t where Λ is the diagonal matrix of the singular values. The number of null singular values gives the dimension of the kernel of S, and the columns of V corresponding to the null singular values form an orthonormal basis of Ker(S) while the columns of U corresponding to the non zero singular values form an orthonormal basis of Im(S). Since S t = V Λ U t , Ker(S t ) and Im(S t ) are characterized as well, allowing a complete determination of the domain of S and of its range.

The SVD of the steady Stokes operator of the channel with its specified boundary condition was performed with Scilab for values of N and M small enough for the problem to remain tractable but large enough to get a generic answer. For the channel flow is suffices that N and M be larger than 3. The number of null singular values and hence the number of modes of the kernel is 2. Although the two modes that come out of the SVD are mixed, the boundary condition mode can be obtained if one assumes (and checks) that the intrinsic mode is the mode (0, 0, 1) and performing the orthonormalisation. It is verified that this mode is the discrete counterpart of the continuous Poiseuille flow, although even in this simple case, it has to be determined numerically due to the location of the velocity component at half grid cells. In the case of a heated vertical plate (Fig. 2), the dimension of the kernel is also equal to 2. The mode of the kernel is a 1D flow with slip at the outer boundary, corresponding to a linear pressure gradient.

Additional computations were performed for other types of geometries. Let us consider the geometry sketched in fig. 3 which corresponds to the case studied by Kettleborough [START_REF] Kettleborough | Transient laminar convection between heated vertical plates including entrance effects[END_REF] in order to set the entrance conditions at artificial boundaries away from the channel entrance. We impose the following boundary conditions 0 .

0 .        V n n V n
where n and τ are the unit vectors normal and tangential to the open boundaries. The use of a cell phase function allowed us to systematize the investigation of the various geometrical shapes shown below.

Figure 3 Sketch of the channel geometry considered in [START_REF] Kettleborough | Transient laminar convection between heated vertical plates including entrance effects[END_REF].

The artificial boundaries are the dotted lines.

Performing the SVD of the Stokes operator in this case reveals that its kernel is of dimension 3, which leaves 2 undetermined Stokes modes in addition to the constant pressure mode. Let us further consider the case for the flow around a heated body such as sketched in Fig. 4. If we impose Neuman conditions on the normal component on all four external boundaries, the number of null singular values is equal to 4 yielding 3 undetermined Stokes modes in addition to the constant pressure mode.

Figure 4 Configuration of a heated body

A SUPERPOSITION METHOD

General: Recognizing that with the specified boundary conditions the steady state solution of the Stokes problem is undefined up to a linear combination of the non-trivial Stokes modes solution of the homogeneous Stokes problem opens the way to the determination of a correct solution with the help of a superposition method. Assuming a particular solution has been obtained, one can superpose to this solution a linear combination of the non-trivial Stokes modes so as to satisfy additional constraints on the solution. The number of these constraints has to be equal to the dimension of the kernel. The additional constraints themselves depend on the modeling assumptions, as will be shown below. The present superposition method will thus allow one to compute the unique solution that satisfies the additional constraints which depend on the modeling assumption.

Channel flow: Assume a particular solution has been obtained. On can superpose to this solution a linear combination of the Poiseuille flow requiring that the pressure drop between inlet and outlet is equal to a specified value, which is the global condition that the channel flow has to satisfy. This value can be either zero, which says that the motion pressure outside the channel remains constant or to -½ G 2 if one takes into account the pressure drop needed to accelerate the fluid from infinity to the channel entrance [START_REF] Dalbert | Influence des conditions d'entrée sur l'écoulement de convection naturelle dans un thermosiphon vertical[END_REF]. One could then proceed as follows. A particular solution is obtained through some procedure, either through time integration or solving iteratively for the steady state equations. When steady state is obtained the mean pressure drop between inlet and outlet is computed and the multiplicative constant of the Poiseuille flow is determined so as to produce the required pressure drop. The Poiseuille flow times this constant is then added to the particular solution, yielding the "true" solution.

In fact it turns out that it is not so easy to compute directly the particular steady solution precisely because this solution is not completely defined. In particular it may continuously keep shifting in time and very often diverge, in particular when the chimney effect is triggered. It is better then to control the solution in a time dependent way by using the superposition principle at each time step. In that case the boundary condition mode is no longer the steady Poiseuille flow but its unsteady counterpart that has to be determined numerically. That mode corresponds to (0, w(x), z) where w is solution of 1)=0. In this version, the superposition method is done at each time step. This results in a stable algorithm. It should be noted that the unsteady Stokes mode has to be computed only once (assuming Δt is constant), and the additional cost of the superposition principle is completely negligible. Sample computations were performed to test the effectiveness of this unsteady superposition principle.

1 2 2   w x w    with w(0)=w(
We consider the specific case of an asymmetric channel (fig. 5) investigated experimentally in [START_REF] Webb | High Rayleigh number laminar natural convection in an asymmetrically heated ver Transfer[END_REF] and was considered as a benchmark exercise in the french convection community [START_REF] Desrayaud | Etude numérique comparative des écoulements thermoconvectifs dans un canal vertical chauffé asymétriquement[END_REF]. The numerical algorithm integrates the equations in unsteady form using a prediction-projection splitting procedure. The prediction step uses an implicit treatment of the viscous terms coupled with an explicit treatment of the convective terms. The projection step requires solving a Poisson type equation for the pressure correction which is done using a multigrid algorithm.

Figure 5

Asymmetric channel hated at uniform heat flux on half its vertical left wall.

We have computed the solution of the governing equations with the specified boundary conditions for one set of values of parameters corresponding to the benchmark exercise: A=5, Ra w =5 10 5 . Fig. 6 presents the w-velocity profile at various heights, showing that the flow enters with a parabolic profile and exits with a boundary layer profile. The plot also shows that above a certain height downward velocities are found which correspond to the existence of a recirculation zone with fluid entering from the top. The existence of such recirculation zones with fluid entering from top has been the subject of long debate in the heat transfer literature (see [START_REF] Kihm | Onset of flow reversal and penetration length of natural convective flow between isothermal v[END_REF]). Needless to say that the confident and accurate numerical simulation of these subtle configurations requires an appropriate treatment of the inlet and outlet conditions. Fig. 7 shows the vertical pressure distribution at mid-width. It starts at -0.5 G 2 , first shows a linear drop corresponding to fully developed flow, followed by a monotonous increase up to the end value.

               0 P m 0 z      0 z w z u                          0 x P 0 w u 1 x m u w z 0 x 2 A 2H/d 1 outlet inle A A/2 A/2 t                  2 m G 5 . 0 P 0 0 u 0 z w         Figure 6
Vertical velocity profiles at various heights We have shown that the imposition of Neumann type boundary conditions may lead to one or several nontrivial combinations of velocity-pressure fields which satisfy the homogeneous Stokes operator, in addition to the unavoidable constant pressure mode. This recognition leads to proposing an algorithm in which the solution is sought as a combination of particular solution of the inhomogeneous Stokes or unsteady Stokes problem plus a linear combination of the modes of the kernel so as to satisfy global conditions such as a mean pressure difference between outlet and inlet. The algorithm was tested on two classical configurations, a channel flow and a heated vertical plate, demonstrating its effectiveness and its efficiency. The algorithm can also more generally apply for pipe networks with one or several outlets, which is the unique way not to specify a priori the

  Figure1Sketch of the channel geometry

Figure 2

 2 Figure 2 Sketch of a heated vertical plate. The artificial boundaries are the dotted lines.

Figure 7

 7 Figure 7 Pressure distribution at mid-width. Convergence with mesh refinement is shown Let us mention that for a Ra w number of 10 7 , a global unsteadiness of the solution was observed as shown in Fig. 8 that presents the time variation of the Nusselt number. The figure shows that the asymptotic evolution of the Nusselt number displays repetitive large amplitude fluctuations of period approximately equal to 400 convective time units. The dynamics is made of a long decrease followed by a sudden jump to the maximum value. We have been able to relate this behavior to a slow increase in time of the

Figure 9

 9 Figure 9 Computational domain for heated plate calculation Fig. 10 presents the vertical velocity profiles at various heights along the plate, showing the development of the boundary layer and that all the flow comes from the lateral boundary since the vertical velocity at z=1/4 is zero. Figure 10 shows that the computed solution verifies the similarity conditions.

Figure 10

 10 Figure 10Vertical velocity profiles at various heights scaled with similarity variables

Figure 11

 11 Figure 11 Sketch of vented cavity

Figure 13

 13 Figure 13 Fully non linear solution for Ra=10 5 Pipe flow with multiple outlets: A classical flow configuration is that of pipe flow with multiple outlets which raises the issue of the way the inlet flow rate splits between the various branches. If one does not want to specify a priori the flow rates in the various branches, one is led to impose Neuman type conditions on the normal component at the various outlets. Consider for instance the flow geometry sketched in figure 14, consisting of one inlet and two successive branches with 4 outlets.

Figure 14

 14 Figure 14 Sketch of pipe flow with 4 outlets

Figure 16

 16 Figure 16Nonlinear solution for a Rey

ACKNOWLEDGMENTS

I would like to acknowledge many stimulating discussions with my friends and colleagues O.


Andrzej J Maciejewski

Maria Przybylska
email: maria.przybylska@proton.if.uz.zgora.pl

Non-integrability of the three body problem

Keywords: integrability, non-integrability criteria, monodromy group, differential Galois group, hypergeometric equation, Hamiltonian equations, three body problem

published or not. The documents may come

Introduction and results

We consider three point masses m 1 , m 2 and m 3 in a plane. We assume that they interact mutually according to a generalized low of the gravitation. Namely, the force of the attraction between two points is proportional to the product of their masses and inversely proportional to a certain power of the distance between them.

Let r 1 := (x 1 , x 2 ), r 2 := (x 3 , x 4 ) and r 3 := (x 5 , x 6 ) denote the inertial Cartesian coordinates of the masses and (y 1 , y 2 ), (y 3 , y 4 ), (y 5 , y 6 ) their respective linear momenta in this frame. Then the Hamiltonian of the problem has the form and

K = 1 2m 1 (y
r 12 := (x 1 -x 3 ) 2 + (x 2 -x 4 ) 2 ,
r 23 := (x 3x 5 ) 2 + (x 4x 6 ) 2 , r 31 := (x 5x 1 ) 2 + (x 6x 2 ) 2 .

We assume that 2n ∈ N. The system admits three additional classical first integrals: two components of the total momentum, and the total angular momentum Y 1 :=y 1 + y 3 + y 5 , Y 2 := y 2 + y 4 + y 6 , C :=x 1 y 2x 2 y 1 + x 3 y 4x 4 y 3 + x 5 y 6x 6 y 5 .

(1.3) However, those first integrals do not pairwise commute, as we have

{C, Y 1 } = Y 2 , and {Y 2 , C} = Y 1 . (1.4)
We perform a canonical reduction and we eliminate two degrees of freedom. To this end we make a linear canonical transformation x = S q, y = S -T p, S -T := S -1 T , (

where q := [q 1 , . . . , q 6 ] T , p := [p 1 , . . . , p 6 ] T ,

S :=         1 0 0 0 -1 0 0 1 0 0 0 -1 0 0 1 0 -1 0 0 0 0 1 0 -1 σ 1 0 σ 2 0 σ 3 0 0 σ 1 0 σ 2 0 σ 3         , ( 1.6)
and

σ i = m i m for i = 1, 2, 3; m = m 1 + m 2 + m 3 .
In other words (q 1 , q 2 ) is the vector between m 3 and m 1 , i.e., (q 1 , q 2 ) := r 1r 3 ; similarly, (q 3 , q 4 ) := r 2r 3 is the vector between m 3 and m 2 , but (q 5 , q 6 ) are coordinates of the mass centre

(q 5 , q 6 ) := 1 m 3 ∑ i=1 m i r i .
The transformed Hamiltonian (1.1) reads

H( q, p) := K(S q, S T p) = H r (q, p) + 1 2m (p 2 5 + p 2 6 ), (1.7) 2
where

q := [q 1 , . . . , q 4 ] T , p := [p 1 , . . . , p 4 ] T , H r = H r (q, p) = T r (p) + U r (q), (1.8) T r (p) := 1 2µ 1 p 2 1 + p 2 2 + 1 2µ 2 p 2 3 + p 2 4 + 1 m 3 (p 1 p 3 + p 2 p 4 ) , (1.9) U r (q) := - m 1 m 2 [(q 1 -q 3 ) 2 + (q 2 -q 4 ) 2 ] n - m 2 m 3 [q 2 3 + q 2 4 ] n - m 3 m 1 [q 2 1 + q 2 2 ] n , (1.10)
and

µ 1 := m 1 m 3 m 1 + m 3 , µ 2 := m 2 m 3 m 2 + m 3 .
(1.11)

Clearly q 5 and q 6 are cyclic coordinates. Moreover, Hamilton's equations with Hamiltonian (1.7) split into a direct product of Hamilton's equations with Hamiltonian H r and Hamiltonian H c := (p 2 5 + p 2 6 )/(2m). The reduced system with four degrees of freedom governed by Hamiltonian H r has one additional first integral

F := q 1 p 2 -q 2 p 1 + q 3 p 4 -q 4 p 3 , (1.12)
which is the total angular momentum of the system. The system given by (1.8) we call the partially reduced three body problem. One can eliminate one more degree of freedom using first integral (1.12). This reduction is described, e.g., in § 161 of [START_REF] Whittaker | A Treatise on the Analytical Dynamics of Particle and Rigid Bodies with an Introduction to the Problem of Three Bodies[END_REF]. The obtained system has three degrees of freedom and we call it the fully reduced three body problem.

We consider the Hamiltonian system generated by (1.8) in the complex phase space which is an open subset of C 8 with canonical coordinates q = (q 1 , . . . , q 4 ) and momenta p = (p 1 , . . . , p 4 ). Our main result is formulated in the following theorem. Theorem 1.1. Assume that 2n ∈ N \ {2} and masses m 1 , m 2 and m 3 are positive. Then the Hamiltonian system given by (1.8) is not integrable in the Liouville sense.

The non-integrability of the non-reduced classical gravitational 3 body problem corresponding to 2n = 1 was proved in the framework of the Ziglin theory in [START_REF] Ziglin | On involutive integrals of groups of linear symplectic transformations and natural mechanical systems with homogeneous potential[END_REF] and later by means of the so-called Morales-Ramis theory in [START_REF] Morales-Ruiz | On the meromorphic non-integrability of some N-body problems[END_REF][START_REF] Simon | On the meromorphic non-integrability of some problems in celestial mechanics[END_REF]. The first proof of the non-integrability of the fully reduced classical planar three body problem was done by Alexei Tsygvintsev in [START_REF] Tsygvintsev | La non-intégrabilité méromorphe du problème plan des trois corps[END_REF][START_REF] Tsygvintsev | The meromorphic non-integrability of the three-body problem[END_REF], and by an application of the differential Galois approach by Delphine Boucher and Jacques-Arthur Weil [START_REF] Boucher | theorem to test the non-complete integrability of the planar three-body problem[END_REF].

In [START_REF] Maciejewski | Partial integrability of Hamiltonian with homogeneous potentials[END_REF], basing on ideas of [START_REF] Morales-Ruiz | On the meromorphic non-integrability of some N-body problems[END_REF], we have found a surprisingly simple proof of the nonintegrability of the classical three body problem. In this paper we generalised this result to a class of potentials of the form (1.2) with an arbitrary positive integer 2n.

Our proof of Theorem 1.1 is based on an application of a general theorem concerning the integrability of homogeneous potentials. In the next section we reformulate this theorem according to our needs. Here it is worth to mention that this general result takes its origin from a brilliant work of H. Yoshida [START_REF] Yoshida | A criterion for the nonexistence of an additional integral in Hamiltonian systems with a homogeneous potential[END_REF]. The main statement of this theorem can derived from an analysis of the monodromy group of the variational equations, see [START_REF] Yoshida | Justification of Painlevé analysis for Hamiltonian systems by differential Galois theory[END_REF][START_REF] Ziglin | On involutive integrals of groups of linear symplectic transformations and natural mechanical systems with homogeneous potential[END_REF], as well as from an analysis of their differential Galois group [START_REF] Morales-Ruiz | A note on the non-integrability of some Hamiltonian systems with a homogeneous potential[END_REF]. Remark 1.2 In Theorem 1.1 we excluded the case of the Jacobi problem when the potential of interactions between bodies is homogeneous of degree -2. In this case, the above mentioned Morales-Ramis theorem does not give any obstruction to the integrability. The non-integrability of the Jacobi three body problem was investigated by Emmanuelle Julliard Tosel in [START_REF] Tosel | Meromorphic parametric non-integrability; the inverse square potential[END_REF]. She proved the non-integrability of the Jacobi problem when masses of two points are equal. It was also proved that the Jacobi problem (with arbitrary masses) is not integrable with rational first integrals which are meromorphic functions of the masses.

Theory

Let us consider a complex Hamiltonian system with n degrees of freedom given by a natural Hamiltonian function of the following form

H = 1 2 p T M p + V(q), (2.1)
where q = (q 1 , . . . , q n ) ∈ C n and p = (p 1 , . . . , p n ) ∈ C n are canonical coordinates and momenta, V(q) is a homogeneous function of degree k ∈ Z ⋆ , and M is a symmetric non-singular n × n matrix. The phase space of this system is C 2n which is considered as C-linear symplectic space with the canonical symplectic form

ω = n ∑ i=1 dq i ∧ dp i .
Hence, Hamilton's equations have the standard canonical form

d dt q = M p, d dt p = -V ′ (q), (2.2)
where V ′ (q) := grad V(q). Moreover, we assume also that the time t is a complex variable.

The basic assumption of our considerations is that there exists a non-zero vector

d ∈ C n such that MV ′ (d) = d. (2.3)
It is called a proper Darboux point of potential V. It defines a two dimensional plane in the phase spaces C 2n , given by

Π(d) := (q, p) ∈ C 2n | q = ϕd, p = ψM -1 d, (ϕ, ψ) ∈ C 2 . (2.4)
This plane is invariant with respect to the system (2.2). Equations (2.2) restricted to Π(d) have the form of one degree of freedom Hamilton's equations

d dt ϕ = ψ, d dt ψ = -ϕ k-1 , (2.5)
with the following phase curves

Γ k,ε := (ϕ, ψ) ∈ C 2 | 1 2 ψ 2 + 1 k ϕ k = ε ⊂ C 2 , ε ∈ C. (2.6)
In this way, a solution (ϕ, ψ) = (ϕ(t), ψ(t)) of (2.5) gives a solution (q(t), p(t)) := (ϕd, ψM -1 d) of equations (2.2) with the corresponding phase curve

Γ k,ε := (q, p) ∈ C 2n | (q, p) = (ϕd, ψM -1 d), (ϕ, ψ) ∈ Γ k,ε ⊂ Π(d).
(2.7)

In order to find necessary conditions for the Liouville integrability we consider the variational equations along an arbitrary phase curve Γ k,ε with ε = 0. These variational equations have the form

ẍ = -ϕ(t) k-2 MV ′′ (d)x, (2.8)
where V ′′ (d) is the Hessian of V calculated at d. A linear change of variables x = Aξ transforms the above system into the following one ξ = -ϕ(t) k-2 Jξ, (

where

J = A -1 MV ′′ (d)A = diag(J i 1 (λ 1 ), . . . , (J i p (λ p )),
C n ∋ q -→ F(q) := MV ′ (q) -q ∈ C n . (2.11) If det F ′ (d) = det(MV ′′ (d) -E n ) = 0, then d is an isolated zero of F. Thus, if d is not isolated, then the Jacobi matrix F ′ (d) is singular. Equivalently, λ = 1 is an eigenvalue of matrix MV ′′ (d) as we claimed.
Necessary conditions for the Liouville integrability of Hamiltonian systems of the form (2.2) which come from an analysis of the differential Galois group of the variational equations (2.9) can be formulated in the following form. Theorem 2.2. Assume that the Hamiltonian system defined by Hamiltonian (2.1) with a homogeneous potential V ∈ C(q) of degree k ∈ Z ⋆ satisfies the following conditions:

1. there exists a non-zero d ∈ C n such that MV ′ (d) = d, and 2. the system is integrable in the Liouville sense with first integrals which are meromorphic in a connected neighbourhood U of phase curve Γ k,ε with ε = 0, and independent on U \ Γ k,ε .

Then for each eigenvalue λ ∈ spectr(MV ′′ (d)), pair (k, λ) belongs to an item of the following table case k λ

1. ±2 arbitrary
2. k p + k 2 p(p -1) 3. k 1 2 k -1 k + p(p + 1)k 4. 3 - 1 24 + 1 6 (1 + 3p) 2 , - 1 24 + 3 32 (1 + 4p) 2 - 1 24 + 3 50 (1 + 5p) 2 , - 1 24 + 3 50 (2 + 5p) 2 5. 4 - 1 8 + 2 9
(1 + 3p) 2

(b) if k = -1, then λ = 1; (c) if k = 1, then λ = 0.
We denote by M k a subset of rational numbers λ specified by the table in the above theorem for a given k, e.g., for |k| > 5 we have

M k = p + k 2 p(p -1) | p ∈ Z ∪ 1 2 k -1 k + p(p + 1)k | p ∈ Z . (2.13)
Let us remark that the above theorem does not give any obstruction for the integrability if k = 2 or k = -2.

Under assumption M = E n , where E n is n-dimensional identity matrix, the first part of Theorem 2.2 coincides with Morales-Ramis theorem [START_REF] Morales-Ruiz | A note on the non-integrability of some Hamiltonian systems with a homogeneous potential[END_REF], and the second part coincides with Theorem 1.3 in [START_REF] Duval | Jordan obstruction to the integrability of Hamiltonian systems with homogeneous potentials[END_REF].

The generalisations of the above cited theorems to the case M = E n are straightforward. Since we assumed that matrix M is non-singular and symmetric, there exists a non-singular matrix B such that (2.15)

B T MB = E n . ( 2
The above shows that if det M = 0, then, in the context considered here, we can assume that M = E n . However, in practise, it is more convenient to avoid the described transformation, as even simple examples show that it introduces quite complicated expressions. It is also instructive to notice the following fact. If d is a proper Darboux point of the potential V, then c := B T d is a proper Darboux point of the transformed potential U(Q). Thus, applying Theorem 2.2 or Theorem 1.1 in the original variables, we use the eigenvalues of matrix MV ′′ (d) which is, in general non-symmetric. On the other hand, while working with the transformed system, we use symmetric matrix U ′′ (c). However, in general matrix U ′′ (c) is a complex matrix, so it is not necessarily diagonalisable. Of those two matrices are similar, as we have

′′ (c) = C T V ′′ (d)C = C -1 (CC T V ′′ (d))C = C -1 MV ′′ (d))C, (2.16)
because matrix M can be written as M = CC T , see (2.14).

In the formulation of Theorem 2.2 we assumed for the simplicity that the considered potential is a rational homogeneous function. However, this theorem remains valid for a class of homogeneous potentials which are meromorphic in a neighbourhood of the considered particular solution.

Proof of Theorem 1.1

The Hamiltonian (1.8) has the form (2.1) with matrix M = M r given by

M r :=            1 µ 1 0 1 m 3 0 0 1 µ 1 0 1 m 3 1 m 3 0 1 µ 2 0 0 1 m 3 0 1 µ 2            . (3.1)
In order to apply Theorem 2.2 we need the existence of a proper Darboux point. In fact such a point exists and it is related to the Euler configuration.

+ ρ), 0), a > 0, ρ > 0, (3.2)
where ρ a unique positive root of the following polynomial (3.3) and a = α

P(ρ) := m 1 (1 + ρ) 2n+1 1 -ρ 2n+2 + m 2 (1 + ρ) 2n+2 -ρ 2n+2 + m 3 ρ 2n+1 1 -(1 + ρ) 2n+2 ,
-1 2(n+1) is given by = ρ 2n+1 (ρ + 1) 2n+1 [m 3 + (m 1 + m 3 )ρ] 2nm [m 3 ρ 2m+1 + m 1 (1 + ρ) 2n+1 ] . (3.4)
Proof. For vector e of the form (3.2) equation ′ r (e) = M -1 r e reduces to the following equations

-m 3 + m 2 ρ + 2αnm m 3 - m 2 ρ 2n+1 = 0, -m 3 (ρ + 1) -m 1 ρ + 2αnm m 1 ρ 2n+1 + m 3 (ρ + 1) 2n+1 = 0. (3.5)
Eliminating α from the above equations we find that ρ is a root of polynomial (3.3), and that α is given by (3.4).

Clearly deg P(ρ) = 4n + 3. Moreover, P(0) > 0, and P(ρ) → -∞ as ρ → ∞. We write P(ρ) the following form

P(ρ) = 4n+3 ∑ i=0 P i ρ i .
From (3.3) it is easy to deduce that P i < 0 for i > 2n + 1, and P i > 0 for i ≤ 2n + 1. Hence we have exactly one change of sign of the coefficients, and this proves that P(ρ) has one positive root. Remark 3.2 It can be easily shown that if a and ρ are such that e given (3.2) is a proper Darboux point, then also e s := (a cos s, a sin s, a(1 + ρ) cos s, a(1 + ρ) sin s), (

is a proper Darboux point for an arbitrary s ∈ R. Hence, e = e 0 is not isolated. Now, our proof of Theorem 1.1 follows the following steps:

1. We show the matrix M r U ′′ r (e) has eigenvalues (1, -(2n + 1), λ, -(2n + 1)λ), where λ > 0.

2. If system is integrable, then λ 1 λ, and λ 2 = -(2n + 1)λ satisfy conditions of Theorem We show that this implies that λ = 1.

3. Finally, we prove that for positive masses m 1 , m 2 and m 3 the equality λ = 1 is impossible.

The next three propositions prove the statements formulated in the above steps. Proof. Potential U r is homogeneous of degree k = -2n. This is why k -1 = -(2n + 1) is an eigenvalue of M r U ′′ r (e). Moreover, by Remark 3.2, e is not isolated proper Darboux point. So, by Remark 2.1, λ = 1 is an eigenvalue of M r U ′′ r (e). Let p(z) := det(M r U ′′ r (e) -zE 4 ) with e given by (3.2) be the characteristic polynomial of matrix M r U ′′ r (e). The coefficients of this polynomial depend rationally on α and ρ. Using (3.4) we eliminate from them α. Next we notice that the coefficients of polynomial P(ρ) given by (3.3) are linear with respect to the masses. Thus, assuming that ρ is the positive root of P(ρ), from P(ρ) = 0, we obtain

m 2 = m 1 r 1 (ρ) + m 3 r 3 (ρ), where r 1 (ρ) := (1 + ρ) 2n+1 (ρ 2n+2 -1) (1 + ρ) 2n+2 -ρ 2n+2 , r 3 (ρ) := ρ 2n+1 (1 + ρ) 2n+2 -1 (1 + ρ) 2n+2 -ρ 2n+2 .
Using the above relation we eliminate m 2 from the coefficients of polynomial p(z). After this operation it factors

p(z) = (z -1)(z + 2n + 1)p 2 (z), where p 2 (z) = z 2 + bz + c. (3.7)
By a direct calculation we found that b = 2nλ, and

c b 2 = - 2n + 1 4n 2 , ( 3.8)
where

λ := m 3 + (m 1 + m 3 )ρ m 1 (1 + ρ) 2n+1 + m 3 ρ 2n+1 R(ρ), ( 3.9)
and Proof. We know that for an arbitrary k ∈ Z \ {-2, 0, 2} we have

R(ρ) := (ρ(1 + ρ)) 2n (1 + ρ) 2n+3 -1 -ρ 2n+3 (1 + ρ) 2n+1 (1 + ρ 2n+1 ) -ρ 2n+1 . ( 3
M (1) k ∪ M (2) k ⊂ M k , where M (1) k := p + k 2 p(p -1) | p ∈ Z , (3.11) and M
(2)

k := 1 2 k -1 k + p(p + 1)k | p ∈ Z . (3.12)
Moreover, for an arbitrary k ∈ Z \ {-2, 0, 2} numbers λ 1 = 1 and λ 2 = k -1 are elements of M

k . Now we assume that k = -2n is an negative integer, and that k = -2. We have to show that if λ 1 = λ > 0, and λ 2 = -(2n + 1)λ < 0, then λ = 1.

For k = -1 we have

M -1 := - 1 2 p(p -3) | p > 1, p ∈ N . (3.13)
This set contains just one positive element that is equal to one. Thus,

λ 1 , λ 2 ∈ M -1 , iff (λ 1 , λ 2 ) = (1, - 2
). This ends the proof for k = -1.

For k = -2n = -3, we have

M -3 = 6 i=1 M (i) -3 ,
where M

-3 := 25 24 - 1 6 (1 + 3p) 2 | p ∈ Z , M (3)
-3 := 25 24 - 3 32 (1 + 4p) 2 | p ∈ Z , M (4)
-3 := 25 24 - 3 50 (1 + 5p) 2 | p ∈ Z , M (5)
-3 := 25 24 - 3 50 (2 + 5p) 2 | p ∈ Z . (3.14) Now, assumption that n > 2 implies that 1 + 8p(1 + p) < 4, (6)
and this in turn forces that p = 0. However it is impossible because it gives a negative n. A contradiction proves our claim.

Proposition 3.5. For positive masses m 1 , m 2 and m 3 , the eigenvalue λ given by (3.9) cannot be equal to one.

Proof. Let us assume that λ = 1, and that masses m 1 , m 2 and m 3 are positive. Then from (3.9) we can determine m 3 and calculate the ratio m 2 /m 3 . After some simple algebra we obtain

m 2 m 3 = (ρ(1 + ρ)) 2n+1 -R(ρ) ((1 + ρ) 2n+1 -ρR(ρ) .
Taking into account the explicit form of R(ρ) given by (3.10) we obtain

m 2 m 3 = - r 3 (ρ) ρ = -ρ 2n (1 + ρ) 2n+2 -1 (1 + ρ) 2n+2 -ρ 2n+2 .
The above ratio is negative for a positive ρ. We have a contradiction with assumption that all masses are positive and this ends our proof.

Remark 3.6 Instead of the last Proposition we can use the following reasoning to conclude the proof of Theorem 1.1. If λ = 1, then matrix M r U ′′ r (e) has double eigenvalue l, and double eigenvalue -(2n + 1). Moreover, it is not diagonalisable, and it has two Jordan blocks of dimension two. It means that those eigenvalues correspond to a non-real Darboux point e. In this case, by the second part of Theorem 2.2, the system is not integrable.

Discussion and comments

For the partially reduced problem we can find other particular solutions. Namely, one can easily prove the following Darboux point c gives a particular solution of the three body problem corresponding to the triangular Lagrange solution. There are two good properties of this solution: we know it explicitly, and moreover we know explicitly the corresponding eigenvalues. Probably this is why it was used in [START_REF] Simon | On the meromorphic non-integrability of some problems in celestial mechanics[END_REF][START_REF] Morales-Ruiz | On the meromorphic non-integrability of some N-body problems[END_REF].

However, we cannot prove the non-integrability using only this solution. Let us consider as example the classical case 2n = 1. Conditions λ ± ∈ M -1 can be written in the following form -27Q = (p + 1)(p -1)(p -2)(p -4), (

for a certain integer p greater than one. As Q > 0, we have only one choice for p, namely p = 3, and this gives Q = 8/27. We can assume that the sum of masses is one. Then, m 1 and m 2 are independent parameters of the problem and condition Q = 8/27 defines a curve on the plane (m 1 , m 2 ). For all points lying on this curve the necessary conditions for the integrability are satisfied.

2 P

 2 .14) Thus, we can make the following canonical transformation q = CQ, p = BP, where C := (B -1 ) T , which transforms the Hamiltonian (2.1) into K(Q, P) := H(CQ, BP) = 1 T P + U(Q), where U(Q) := V(CQ).

Proposition 3 . 1 .

 31 Assume that 2n is a positive integer and masses m 1 , m 2 and m 3 are positive. Then equation M r U ′ r (e) = e has a solution e of the form e := 0, a(1

Proposition 3 . 3 .

 33 Assume that 2n is a positive integer, and e be the solution of M r U ′ r (e) = e given by Proposition 3.1. Then matrix M r U ′′ r (e) has eigenvalues (1, -(2n + 1), λ, -(2n + 1)λ), where λ > 0.

. 10 )Proposition 3 . 4 .

 1034 Evidently, we have λ > 0, and this finishes our proof. Assume that k = -2n is a negative integer and k = -2. Then λ ∈ M k and (k -1)λ ∈ M k if and only if λ = 1.

Proposition 4 . 1 .b 2 2 =

 412 Vectorc := α(a 1 , b 1 , a 2 , b 2 ) T , α 2(n+1) (a 1a 2 ) 2 + (b 1b 2 ) 2 = 1. (4.3)Moreover, matrix M r U ′′ r (c) has eigenvalues (-2n -1, 1, λ -, λ + ), whereλ ± := -n ± (n + 1) 1 -3Q, Q := 1 m 2 (m 1 m 2 + m 2 m 3 + m 3 m 1 ). (4.4)

Acknowledgements

This research has been partially supported by grant No. N N202 2126 33 of Ministry of Science and Higher Education of Poland and by EU funding for the Marie-Curie Research Training Network AstroNet.

Elements of sets M (i)

-3 have the following properties

, (s, 600) = 1.

( show that if

is not an element of M -3 . In fact using properties (3.15) we easily show that there is only one possibility λ 1 = 2/3. However, this gives λ 2 = -8/3. Thus, by properties (3.15)

-3 but it is not true, as one can check directly. This end the proof for k = -3.

For k = -2n = -4 we have

where

We can sort all elements of the sets

-4 := Thus we have

Just a direct inspection of (3.17) shows that λ 2 = -5λ ∈ M -4 and this ends the proof for k = -4

For k = -2n = -5 we have

where

(3.18) Proceeding in the way similar to the previous case we prove the statement for k = -5.

For k = -2n < -5, we have

k is a positive number, then

Hence, either p > 0 and p < 1 + 1 n , or p < 0 and p > 1 + 1 n .

As n > 2, the only possibility is p = 1, and this gives λ = 1.

k is a positive number, then 4nλ = 2n + 1 -4n 2 p(p + 1) > 0, for a non-negative integer p. The only possibility is p = 0, and this gives

We show that

k because λ 2 is not an integer. We show that λ 2 ∈ M

(2)

k , then there exists a non-negative integer p such that -(2n + 1) 2 = 2n + 1 -4n 2 p(p + 1).

From the above equation we find that n = 2 1 + 8p(1 + p) -3 .