Démonstration :
On rappelle la définition de Jn(z) :
et le développement classique de 1/[1-e.cosE] -1 , fonction paire périodique de moyenne nulle vaut: 
formula_60
avec formula_61 
car formula_62
Série entière de l'excentricité.
C'est encore Lagrange qui trouve la solution en inventant pour l'occasion son théorème d'inversion des fonctions holomorphes ; et Laplace donnera le rayon de convergence : mais Cauchy, pas content du tout, fonde la théorie des séries analytiques pour résoudre ce problème épineux, qui verra son aboutissement avec les travaux de Puiseux.
avec formula_64 et D := opérateur dérivée.
C'est l'application du théorème d'inversion de Lagrange.
Cas des comètes : formula_65.
Le premier à se confronter au problème est Horrocks, puis surtout Halley (1705), pour les calculs sur sa comète d'excentricité e = 0,9673.
Il faut modifier légèrement la solution de Barker (e = 1). Et Bessel(1805) résout ce cas, mais pour e > 0.997
Gauss (1809) s'illustra en donnant une belle solution pour 0,2 < e < 0,95
Autant dire que le voisinage de (0,95 ; 0,98) est fertile en problèmes, en cas d'itération !
Calcul numérique.
Les calculs via les intégrateur symplectiques exigent de rester toujours en butée du nombre de digits, dans le moindre coût de calcul.
Depuis 300 ans, on cherche la « meilleure » méthode. Elle reste à trouver !
Bien sûr, cela dépend beaucoup du doublet (M,e), M compris entre 0 et Pi et de e, surtout quand e est voisin de 1. 
Nijenhuis (1991) adopte la méthode de Mikkola (1987) qui est la méthode de Newton d'ordre 4, en choisissant « adéquatement » le germe Eo en fonction du doublet (M,e).
Il est clair que dans les calculs numériques, le volume de calculs est essentiel, autant que le nombre de décimales, vu l'instabilité du système solaire évaluée à un coefficient de Liapunov de 10^(t/5Myr). On se heurte à une muraille exponentielle : difficile d'aller plus loin que 25 Myr, même avec un traitement 128 bits.
Ce sont ces calculs (astronomiques... mais informatisés) qui tournent sur les machines de l'IMCCE-Paris. Le calcul de l'ensoleillement terrestre à la latitude 65°Nord, I(65,t) est calculé et on essaie d'en déduire la corrélation avec le climat passé : l'échelle géologique jusqu'au Néogène (25M ans) en est déduite(échelle géologique Gradstein 2004). Prochaine étape prévue : les 65 M ans.
Histoire des sciences.
Avant Kepler, l'équation est déjà étudiée ! bien sûr, pas pour le même problème, mais pour la même équation : 
c'est le problème de la réduction des coordonnées locales aux cordonnées géocentriques : il faut réduire la correction de parallaxe. Habash al Hasib s'y est déjà attaqué.
Avant 1700, il y a déjà beaucoup de tentatives : Kepler naturellement, Curtz (1626), Niele, Bouillau (1645, 1657), Seth Ward (1653), Paganus (1657), Horrebow (1717), Cassini (1669), Newton (1665?), Wren (1658), Wallis (1659)... De toutes, celle de Jeremiah Horrocks (1638) est de plus grande beauté. Cf le Colwell, déjà cité.
compléments.
En 1770, Lagrange trouve les deux séries, mais le changement des termes dans les séries le laisse perplexe. 1821 : Cauchy enfin ! Sitôt après, 1824, Bessel (1784-1846)fera une étude extensive de "ses" fonctions , déjà apparues en 1703 dans une lettre de jean Bernouilli à Leibniz. Daniel Bernouilli fait la théorie du mode propre de la corde suspendue et introduit Jo(x); Euler généralisant a besoin des In(x) , les bessel-modifiées.
On montre que E(Horrocks) = M + e/1sin M +e²/2 sin 2M +e³/3 sin3M +... et E(H) -E = 1/6 .e³sin³M ; pas si mal!
la fonction étant croissante , on "tire" trop bas avec x0 (F(x) est négatif), trop haut avec x1 (F(x) est positif) : alors la racine est entre les deux et on prend la corde.
x = a ch u et y = b sh u ; r = a(1-e ch u)
On pose 1/cos g = ch u et tg g = sh u 
soit g = gudermannien (u) = gd(u) = 2 arctg(exp u) -Pi/2.
En attente , les perturbations , pour faire de mon mieux.
les perturbations du mvt de Kepler sont parmi les plus "dures" car il y a la dégénérescence banale de SO(3), mais aussi la dégénéréscence de SO(4) pour les états liés : du coup il faut comprendre la structure de la sphère S3 dont on sait qu'elle se retourne comme un gant ou peut se transformer en une foliation torique de Hopf, etc. Comment la perturbation agit sur chacun de ces aspects est encore à inventorier, même si on en connaît pas mal sur le sujet, en particulier gràce aux travaux de Poincaré, KAM, Mather, etc. Il est vrai que le niveau est plus élevé ici, puisqu'il s'agit de problèmes le plus souvent non intégrables.
à la manière directe : Danjon-Pollard-Duriez.
la perturbation F est installée au temps t=0 , avec OMo et Vo donnés , càd Lo,Eo et eo données et passage au périgée donné.On appellera ko la direction de Lo, et uo = OMo/ro, et formula_66 pour compléter le trièdre, dont le vecteur-rotation instantanée sera formula_67 ( v signifiera donc vecteur-vitesse). Sept équations sont bien compréhensibles :
Moins évidente est la variation de l'anomalie moyenne :
Il en résulte les équations de Gauss.
équations de Gauss.
le quintuplet [a,e,i,formula_72]s'en déduit projeté sur le reférentiel initial et final :
Et il reste encore dM/dt à écrire !
Comme de plus il faut projeter l perturbation sur la base initiale et la base finale , l'interprétation est sévère.
heureusement, la perturbation dérive souvent d'un potentiel : cela simplifie l'écriture et la compréhension de ces 6 équations, sur lesquelles il faut se pencher qq temps pour les assimiler.
pertinence des équations de Gauss.
demandée ici, pour "souffler un peu" : le cours est construit ainsi ! ne rien faire que l'on ne puisse refaire ou retenir ! Pour retenir, il faut manipuler et croiser les équations jusqu'à ce que cela devienne "machinal" et au fond "intuitif" . Donc la question posée est : en quoi les 6 équations précédentes vous semblent-elles pertinentes ?
Perturbation de Kepler : effet Stark classique.
Si à la force newtonienne vient se rajouter une petite force F, la trajectoire va être légèrement perturbée. Néanmoins si F est parallèle au vecteur excentricité, la symétrie ne sera pas entièrement détruite.
Il convient de prendre les bonnes coordonnées pour traiter ce problème. Comme on sait traiter le mouvement keplerien en système de coordonnées paraboliques, il faut évidemment en profiter.
Mais si F devient trop grand, il apparaît clairement que l'atome va pouvoir s'ioniser plus facilement. 
En mécanique quantique cela sera encore plus évident via l'effet tunnel, conduisant à l'ionisation Stark, fragilisant surtout les atome de Rydberg.
Mouvement d'Euler à 2 centres d'attraction.
Euler a vite compris que la composante du vecteur excentricité permettait d'intégrer le problème à 2 soleils fixes et une planète. Cela s'opère grâce à un système de coordonnées bifocales.
Vinti s'est fait le promoteur de cette méthode : ébauche
Mouvement si Terre-galette (Béletskii).
Beletskii a fait remarquer que le problème d'Euler pouvait s'appliquer à un Soleil légèrement allongé de forme cigare. Par prolongation analytique, avec des masses « imaginaires », il a proposé une interprétation simple du mouvement d'un satellite terrestre sous l'action perturbante du bourrelet (le terme J2(P2(cos(theta)/r³) dans le potentiel gravitationnel. On retrouve les effets décrits dans satellite artificiel.
Perturbation de Kepler par planète proche : Terre & Lune.
Ce problème est ardu : Newton disait que cela lui donnait mal à la tête.
Il a fallu attendre Clairaut (1741) pour avoir une première théorie de la Lune.
Aujourd'hui avec les miroirs posés sur la Lune (Apollo et Lunakhod), on peut comparer la théorie analytique à celle numérique. La précision théorique des LLR (laser lunar range: tir laser vers la Lune) est de quelques centimètres. La théorie analytique comprend plusieurs milliers de termes, mais donne aussi une précision de quelques mètres.
à compléter (séminaire Laskar du 09/03/06).
Perturbation de Kepler par planète lointaine : Terre & Jupiter.
Là, le problème est plus facile . L'essentiel de la méthode consiste en une méthode variable rapide- variable lente, due à Legendre, puis Gauss.
à compléter.
Perturbation de Kepler et symétries.
Bien sûr, chaque fois qu'un système possède une symétrie continue, le théorème de Noether donne une intégrale première, ce qui permet d'éliminer une variable de l'espace des phases.
Comment s'opère cette réduction ? 
Le livre de Cordani, celui de Marsden & Ratiu expliquent cette réduction.
Enfin, le problème garde toujours sa symétrie symplectique : il faudra expliquer comment fonctionnent les intégrateur symplectique (Laskar & Robutel, Celestial Mechanics, 2001,80, 39-62). 
Applications.
Elles sont innombrables : 
Perturbations du mouvement de Kepler.
C'est évidemment essentiel.
Pour les satellites artificiels, il faut tenir compte de la forme non sphérique de la Terre , ET de toutes les autres petites perturbations ( pression de radiation du Soleil sur les panneaux solaires, action de gravité différentielle de la Lune et du Soleil, etc.
Pour l'astronome , il y a essentiellement deux problèmes : 
À l'heure actuelle, les programmes de calculs peuvent envisager de traiter (sur un temps pas "trop grand") le mouvement de l'ensemble des planètes. On sait depuis peu que Pluton n'est pas une vraie planète. Ceci dit, le mouvement des planètes sur des échelles de qq 10^6 années commence à être sensible aux conditions initiales (la Terre est un cas particulier car la Lune vient stabiliser son inclinaison et son excentricité).
Pour le programme Galileo (le GPS européen), la précision sur le positionnement de la constellation de satellites artificiels est assez impressionnante(inférieure au centimètre).
insert provisoire:atome d'hydrogène.
Cet article suit l'article atome d'hydrogène.
La résolution de l'équation de Schrödinger, écrite en coordonnées polaires, se découple des variables (formula_79) et conduit à une équation à une dimension en r, appelée équation radiale de Leibniz-Schrödinger, puisque ce n'est jamais que la célèbre équation de Leibniz de 1685 traduite en mécanique quantique.
Mais l'équation de Schrödinger (1926) peut se résoudre autrement comme Pauli l'a montré en 1925 !
Équation radiale.
L'équation radiale 1D de Leibniz-Schrödinger s'écrit pour r>0:
avec E valeur propre négative ,
et S(r) s'annulant "vite" à l'infini, et S(0) =0 :il s'agit donc d'un problème aux limites dit de Sturm (par opposition à un problème aux conditions initiales, dit de Cauchy); de plus formula_80.
[On reconnaît dans formula_81 la barrière centrifuge de Leibniz (l entier positif) (l=0 correspond à L =0 ; le problème classique n'a pas de correspondant simple en mécanique quantique, encore que ...)].
Il SUFFIT de chercher la fonction d'onde formula_82 en coordonnées sphériques sous la forme :
où les Y(l,m) sont les fonctions harmoniques sphériques. On appelle ce procédé courant dans les équations aux dérivées partielles, la séparation des variables. Souvent, on appelle R(r) := S(r)/r , la partie radiale de la fonction d'onde.
Harmoniques sphériques.
Il n'y a rien de mystérieux (et surtout rien à voir avec la MÉCANIQUE quantique) dans ce qui semble être un tour de passe-passe. L'étude en électrostatique classique de l'opérateur Laplacien conduit à ces mêmes fonction Y(l,m) , appelées harmoniques sphériques, qui sont des fonctions usuelles dès que la symétrie sphérique entre en jeu. L'entier relatif m ne peut prendre que 2l+1 valeurs, de m = -l à m = +l , l étant un entier positif. 
Ce sont ces harmoniques sphériques qui "quantifient" le problème sphérique par les deux nombres quantiques l et m (comme il est usuel dans tout problème de Sturm, dit "aux limites", des équations différentielles), ces deux entiers l et m qui auront tant d'importance dans l'étude de l'atome à N électrons et donc de la Classification périodique.
Dorénavant, nous n'indiquerons plus ce facteur dit de normalisation.
formula_86 ;
formula_87;
formula_88;
soit la base {x,y,z} dite orbitales formula_89, formula_90, formula_91
formula_92 ;
formula_93; et avec moins , 2i yz ;
formula_94;
formula_95;
Soit la base {3z^2-r^2, xz, zy, yx, x^2-y^2) dont chaque fonction est de laplacien nul.
soit la base { z(5z^2-3r^2), x(5z^2-3r^2), y(5z^2-3r^2),zxy,z(x^2-y^2),x(x^2-y^2), y(x^2-y^2)}dont chaque fonction est de laplacien nul.
C'est ce théorème qui est sans arrêt utilisé pour la théorie de l'atome d'hydrogène.
En chimie ,on représente souvent les fonctions 1/r^(l+1) . Pl comme les harmoniques sphériques des orbitales l ; parfois on prend leur carré; etc. 
Dans l'atome à N électrons pour N< 119, l< 5 : donc cela suffit au physicien de l'atome, qui leur a donné des noms et des représentations mnémotechniques diverses. Ne pas oublier que l'on peut combiner à volonté ces fonctions, pour former ce que les chimistes appellent des orbitales hybridées du sous espace propre du niveau d'énergie En( en particulier les fameuses orbitales paraboliques de Kleinert).
Multiplicité (2l+1).
Le nombre quantique l est appelé nombre quantique azimutal (on voit qu'il joue, par son terme l(l+1), le même rôle que le carré du moment cinétique, L², en mécanique classique). Évidemment l'équation radiale a ramené le mouvement à UNE seule dimension, la variable radiale, avec la fonction S(r) qui doit s'annuler en r=0 (n'oublions pas c'est S(r)/r qui intervient ) et qui doit être de carré sommable sur l'intervalle r>0 . 
On aura donc des valeurs propres de cette équation linéaire, dépendant donc de l , formula_97 , mais pas de m (on dit que la multiplicité de la valeur propre est : 2l+1 ; en physique & chimie on dit : il y a dégénérescence du multiplet égale à 2l+1). 
Le nombre quantique m s'appelle nombre quantique magnétique, car sous l'effet d'un champ magnétique (effet Zeeman) l'énergie dépend alors de la valeur de m, et l'on voit une multiplicité de niveaux d'énergie, d'où la dénomination .
Enfin le nombre k , entier positif, s'appelle nombre quantique radial et donne le nombre de nœuds (k pour knots !) de S(r) pour r > 0 .
Comme la spectroscopie est née un siècle avant la mécanique quantique, la tradition est restée d'appeler le nombre quantique azimutal l par des lettres latines : 
l= 1 -> s ; l=2-> p ; l=3 -> d ; l=4 -> f et ensuite g, h .
Résultat final.
Au final, on trouve une énergie E(l,m,k) indépendante de m, soit E(l,k), mais, de façon incroyable (sauf pour Pauli), ne dépendant que de la somme l+k-1 = n , qui doit être un entier positif, et appelé nombre quantique principal.
C'est la fameuse équation déjà trouvée par Bohr en 1913:
Il y a ce qu'on appelait une dégénérescence accidentelle, avant l'introduction par Pauli en mécanique quantique du vecteur invariant de Runge Lenz.
La multiplicité, g, du niveau d'énergie En est donc :
pour l variant de 0 à n-1 et
pour m variant de -l à +l
formula_98 .
Et, compte-tenu du spin (1/2) de l'électron ,g vaut le double , soit 2.n²
Inutile d'aller plus loin pour décrire la classification périodique, la configuration de l'élément Z= 119 est celle d'un alcalin :
formula_99 soit Kr(Z=36) puis,
formula_100 Rn(Z=86) , puis
formula_101 Uuo(Z=118), 
puis 8s.
Sur les 64 orbitales de la couche S, n= 8 , on n'a besoin de connaître que l'orbitale (8s): ce calcul requiert impérativement la mécanique quantique relativiste , car les électrons (1s²) de la première couche sont soumis à des vitesses non négligeables devant c .
De même, la configuration de l'élément Z= 121 est Uuo,(8s²,5g), la sous-couche 5g pouvant contenir jusqu'à 2*9 =18 électrons.
Ce faisant, on obtiendra ainsi tous les niveaux d'énergie des éléments ET des séries isoélectroniques, ce qui permettra de décrire certains traits de la classification périodique.
Auparavant, on va finir le raisonnement de Schrödinger (1926) ; puis on reviendra sur celui, plus subtil, de Pauli (1925).
Équation radiale-réduite et Polynômes de Laguerre.
Si l'on revient à l'équation radiale de Leibniz-Schrödinger, on peut démontrer que pour r voisin de zéro, S(r) varie comme r^(l+1) , et que pour r très grand, S" + 2E S = 0 . 
Il est courant de poser 2E = -1/n² et donc S(r) varie comme exp (- r/n) à l'infini : n pour l'instant n'est qu'un réel!
Alors le dernier changement de fonction inconnue est logiquement l'essai suivant qui se révèle fructueux : S(r) = r^(l+1).exp(-r/n).g(r) ; mais on s'aperçoit qu'en changeant la variable r en s : = 2r/n l'équation s'arrange mieux :
L'équation radiale-réduite devient :
s f"(s) + (2l+2-s) f'(s) + (n-l-1) f(s) = 0 , avec g(r) = f(2r/n) = f(s) 
Les matheux et Schrödinger ont reconnu cette équation immédiatement (?) : elle conduit à la fonction hypergéométrique dégénérée de Kummer, qui conduit aux polynômes de Laguerre, ssi n-l-1 est un entier positif : donc n est un entier positif et l = 0, 1 , 2 , .. ,n-1 . Et le nombre k est simplement k = n+l-1.
Infeld-Hull et la "factorisation".
dans RevModPhys 23,1951,21-68 , on constate que la méthode des opérateurs d'échelle était bien connue à l'époque (cf aussi Durand, CRAS1950,230,273):
L'idée est classique : 
soit A = 1/2 -a/r -d/dr et B = 1/2 - b/r +d/dr en unités "bien choisies".
A et B sont opérateurs sur les fonctions de carrés sommables sur [0, infty[.
Ils sont opérateurs conjugués pour a = b .
et l'équation de Leibniz-S s'écrit :
A(l+1)B(l+1) Snl = (n-l-1)Snl/r
En multipliant par Snl et en sommant il apparaît immédiatement que n-1> l ;
et B S = 0 pour l = n-1 d'où la valeur de S "circulaire" :
S(r) = r^n .exp (-r/2)
Qq calculs permettent de trouver que 
S(n+1, l) = r A(n) S(n,l)
S(n-1,l) = rB(n) S(n,l) .1/[(n-1-l)n+l)]
et toutes sortes de relations sur les polynômes de Laguerre.
Noter aussi que l'équation du second ordre peut s'écrire , comme assez souvent : 
K(n,l) S(n, l-1) = A S(n,l)
K(n,l) S(n,l) = B S(n, l-1) (Durand p 449)
k+1)<r^k> -2n(2k+1)<r^(k-1)> +[(2l+1)²-k²]<r^(k-2)> = 0
Certaines se trouvent dans atome d'hydrogène
Invariant de Runge Lenz, quantique.
Champ coulombien.
les valeurs propres de l'énergie ne dépendent pas séparément de n et de l , mais seulement de n , entier positif, qui de ce fait est appelé nombre quantique principal de couche (avec n= 1 -> couche K , n=2 -> couche L ..). 
Ce fait, très exceptionnel pour l'énergie, ne sera plus vrai pour un potentiel V(r) quelconque, même voisin de -e²/r. Il convient donc de ne pas trop s'y attacher, sauf si l'on veut s'expliquer cette dégénérescence (anciennement appelée dégénérescence accidentelle), via le raisonnement de Pauli.
vecteur excentricité quantique.
Le vecteur excentricité (cf. mouvement keplerien et invariant de Runge Lenz)vaut :