Calcul \(2s\) pour \(H_2^+\) : intégrales analytiques et petites distances

Oui, on peut faire nettement mieux qu’un calcul « 100 % numérique ».

1. Les intégrales sont en grande partie analytiques

Les orbitales hydrogénoïdes \(2s\) sont des fonctions très simples :

\phi_{2s}(\mathbf r) = \frac{1}{4\sqrt{2\pi}\,a_0^{3/2}} \left( 2-\frac{r}{a_0} \right) e^{-r/(2a_0)}.

Les intégrales à deux centres contiennent donc essentiellement :

des polynômes en \(r_A\) et \(r_B\) ;
des exponentielles ;
les facteurs coulombiens \(1/r_A\) ou \(1/r_B\).

En coordonnées elliptiques prolates,

\xi=\frac{r_A+r_B}{R}, \qquad \eta=\frac{r_A-r_B}{R},

on a :

r_A=\frac{R}{2}(\xi+\eta), \qquad r_B=\frac{R}{2}(\xi-\eta),

et l’élément de volume devient :

d\tau = \frac{R^3}{8} (\xi^2-\eta^2)\, d\xi\,d\eta\,d\varphi.

Les fonctions exponentielles se simplifient fortement :

e^{-r_A/(2a_0)} e^{-r_B/(2a_0)} = e^{-R\xi/(2a_0)}.

L’intégration sur \(\eta\) porte alors sur des polynômes et peut souvent être effectuée analytiquement. Il reste au maximum des intégrales unidimensionnelles très régulières en \(\xi\), voire des expressions fermées formées de polynômes et d’exponentielles.

Ainsi :

\(S(R)\) peut être obtenu analytiquement ;
les intégrales de Coulomb et de résonance peuvent être réduites analytiquement ou semi-analytiquement ;
il n’est pas nécessaire de conserver l’intégration cylindrique bidimensionnelle actuelle.

2. La difficulté aux petites distances n’est pas seulement l’intégration

Même avec des intégrales parfaitement exactes, la formule

E_u = \frac{H_{AA}-H_{AB}}{1-S}

devient numériquement délicate lorsque \(R\to0\), car :

S\to1,

donc :

1-S\to0,

et simultanément :

H_{AA}-H_{AB}\to0.

Le quotient possède une limite finie, mais son calcul naïf soustrait deux nombres presque égaux au numérateur et au dénominateur.

C’est une perte de précision par annulation catastrophique.

3. La branche antisymétrique ne disparaît pas lorsque \(R\to0\)

La combinaison normalisée est :

\psi_u = \frac{\phi_A-\phi_B} {\sqrt{2(1-S)}}.

Pour une petite distance \(R\), on peut développer :

\phi_A-\phi_B \simeq R\frac{\partial\phi_{2s}}{\partial z}.

En même temps :

1-S\simeq C R^2.

Ainsi, le facteur \(R\) du numérateur est compensé par celui provenant de la normalisation :

\sqrt{1-S}\propto R.

La limite de \(\psi_u\) est donc une fonction finie, proportionnelle à :

\frac{\partial\phi_{2s}}{\partial z}.

Cette limite possède une symétrie de type \(p_z\). Le passage sous la branche \(g\) peut donc avoir une véritable signification dans cette base : à courte distance, la branche antisymétrique ne ressemble plus simplement à « deux \(2s\) qui s’opposent », mais tend vers une fonction angulaire différente.

Il faut néanmoins rappeler qu’à \(R\to0\), le système physique tend vers l’ion unicentrique \(\mathrm{He}^+\), alors que notre base conserve des orbitales d’hydrogène de charge nucléaire \(Z=1\). La base minimale n’est donc pas optimale dans cette limite.

4. Ce que je propose pour le nouveau calcul \(2s\)

Je proposerais une version en trois niveaux de robustesse.

A. Coordonnées elliptiques prolates

Calculer les intégrales directement avec :

1\leq\xi<\infty, \qquad -1\leq\eta\leq1.

Cela élimine les singularités mobiles de la grille cylindrique.

L’intégration sur \(\varphi\) est immédiate et celle sur \(\eta\) pourra être faite analytiquement.

Il ne restera qu’une quadrature régulière en \(\xi\), par exemple après le changement :

\xi = 1+\frac{t}{1-t}, \qquad 0\leq t<1.

B. Calcul direct des combinaisons normalisées

Au lieu de calculer séparément :

H_{AA}-H_{AB}

puis de diviser par \(1-S\), on peut calculer directement :

E_u = \langle\psi_u|H|\psi_u\rangle,

avec :

\psi_u = \frac{\phi_A-\phi_B} {\sqrt{2(1-S)}}.

Cela ne supprime pas entièrement la difficulté, mais évite certaines soustractions effectuées après l’intégration.

C. Développement asymptotique pour les petites distances

Sous une valeur comme :

R<0{,}5a_0,

on pourra utiliser des séries :

1-S = c_2R^2+c_4R^4+\cdots,

et :

H_{AA}-H_{AB} = d_2R^2+d_4R^4+\cdots.

Alors :

E_u(R) = \frac{d_2+d_4R^2+\cdots} {c_2+c_4R^2+\cdots},

sans soustraction de nombres presque égaux.

5. Vérification de la réalité de l’inversion

Pour être certain du résultat, la nouvelle version devrait afficher en supplément :

\Delta E(R) = E_u(R)-E_g(R).

Le croisement est donné par :

\Delta E(R_c)=0.

Un petit panneau ou un affichage numérique donnerait :

la distance du premier croisement ;
la distance du minimum local de \(E_g\) ;
la valeur de \(1-S\) ;
éventuellement une alerte lorsque la base devient presque linéairement dépendante.

On pourrait aussi comparer trois précisions de quadrature :

N_\xi=80,\qquad120,\qquad180,

et vérifier que la position du croisement ne bouge pratiquement plus.

Proposition concrète

Je ferais maintenant :

molecule_h2p_02_02_2s_elliptique_SI.py

avec :

coordonnées elliptiques prolates ;
intégration analytique en \(\eta\) lorsque possible ;
quadrature unidimensionnelle en \(\xi\) ;
précision renforcée aux petites distances ;
courbes jusqu’à environ \[ R=0{,}05a_0 ; \]
tracé supplémentaire de \[ E_u-E_g ; \]
indication précise des croisements et minima ;
conservation de l’interface actuelle.

C’est la meilleure étape suivante : elle permettra de déterminer si le minimum vers \(5{,}7\,\text{Å}\) et l’inversion à plus courte distance sont stables, tout en montrant jusqu’où la base minimale \(2s_A,2s_B\) reste physiquement interprétable.