On reformule ci-après le système des
équations de Burgers pour les $x_i$ et $\mu$.
Le principe de la méthode numérique, d\^u à Burgers \cite{b69},
Iben \& MacDonald \cite{id85},
consiste à remarquer que la résolution d'un système
linéaire permet d'écrire la vitesse de diffusion $w_i$ de la particule $i$
sous la forme:
\[w_i=v_i+\sum_j b_{ij}\ddx {x_j}x.\]
Iben \& MacDonald  \cite{id85} n'ont pas tenu compte des
``residual heat flow vector''
qui furent introduits par la suite par Cox \etal{cgk89} et repris
par Thoul \etal{tbl94}.

On utilise les notations suivantes:
\begin{itemize}
\item $c$ nombre d'ions, il y a $c+1$ espèces de particules,
\ie\ ions+ électrons,
\item $E$ champ électrique,
\item $e$ charge de l'électron, en unité électrostatique,
 et indice de la particule "électron",
\item $g_{\rm G}$ accélération due à la gravité,
\item $g_\Omega$ accélération centrifuge,
\item $g_{\rm e}$ accélération effective,
\item $g_{\rm Ri}$ accélération radiative sur la particule $i$,
\item $g_i$ accélération due à la gravité effective
et aux forces radiatives sur $i$,
\item $g_{\rm t}$ accélération totale,
\item $g_{j,i}$ poids statistique du fondamental du niveau d'ionisation
$j$ de l'élément chimique $i$,
\item $m_e=\nu_e m_u$ masse de l'électron,
\item $m_i \equiv \nu_im_u$ masse de l'élément chimique $i$,
\item $m_{ij}\equiv\frac{m_im_j}{m_i+m_j}$ masse réduite des particules
 d'indices$i$ et $j$,
\item $m_u=1/N_0$ masse atomique unité, inverse du nombre d'Avogadro,
\item $n_e$ nombre d'électrons par unité de volume,
\item $n_i$ nombre d'ions de l'élément chimique $i\not=e$
par unité de volume,
\item $n_{j,i}$ nombre d'ions de l'élément chimique $i$
par unité de volume, dans l'état d'ionisation $j=0,\ldots,Z_i$,
\item $r$ rayon, variable d'espace eulérienne,  
\item $T,\ P,\ \rho$ température, pression gaz parfait, densité,
\item $t$ temps,
\item ${\cal V}=\dsdx rt$ vitesse d'entrainement du fluide, ${\cal V}\sim 0$
(hypothèse de l'équilibre quasi-statique),
\item $w_i$ vitesse de diffusion de la particule $i$ par rapport à la vitesse
moyenne du fluide,
\item $X_i$ proportion en masse de l'abondance de l'élément chimique $i$,
$\sum_iX_i\equiv1$,
\item $x_i\equiv X_i/\nu_i$ proportion par mole de l'abondance de
l'élément chimique $i$,
\item $x_{j,i}$ taux d'ionisation du niveau d'ionisation
$j$ de l'ion $i$,
\item $Z_i$ charge de la particule $i$, $Z_e=-1$,
\item $\bar{Z_i}$ charge moyenne de l'élément $i$ avec ionisation,
\item $\eta$ paramètre de dégénérescence,
\item $\mu$ poids moléculaire moyen,
\item $\nu_e$ "masse atomique" de l'électron, $\nu_e=m_e/m_u$,
\item $\nu_i$ masse atomique de l'élément chimique $i$,
\item $\chi_{j,i}$ potentiel d'ionisation du niveau d'ionisation
$j$ de l'élément $i$.
\end{itemize}

\PB{Les notations ne sont pas totalement cohérentes en ce qui concerne
 les électrons.}

\subsection{Charge moyenne des ions}\label{sec:saha}
On utilise l'équation de Saha (Cox \& Guili \cite{cg68}, formule 15-30)
dans laquelle on limite les fonctions de partition aux poids statistiques
des niveaux fondamentaux:
\begin{equation}\label{eq:saha0}
\frac{n_{j-1,i}}{n_{j,i}}=\frac{g_{j-1,i}}{g_{j,i}}\exp\left(\eta-\Delta\mu
+\frac{\chi_{j,i}}{kT}\right)\equiv\varphi_{j,i},\ j=1,\ldots,Z_i.
\end{equation}
Le taux de dégénérescence $\eta$ vérifie:
%\begin{equation}\label{eq:ne}
\[F_{\frac12}(\eta)=\frac{n_e}{4\pi}\left(\frac{h^2}{2m_ekT}\right)^{\frac32},\]
%\end{equation}
$F_{\frac12}(\eta)$ est la fonction de Fermi-Dirac (Clayton \cite{c68},
formule 2-57). $\Delta\mu$ est une correction numérique introduite par
Eggleton \etal{eff73} pour simuler l'ionisation de pression et éviter la
recombinaison à haute température (Kippenhahn \& Weigert \cite{kw91} \S
14.6):
\[\Delta\mu=n_e c_F \left(\frac{a_0}{\bar Z}\right)^3(1+\frac{20\chi_0}{kT}),
\ a_0=\frac{h^2}{4\pi^2m_ee^2},\ \bar Z=X+2Y+8Z,\ c_F=15,\ \chi_0=13.6.\]
%Cette correction est utilisée dans la routine {\tt saha\_eff}.

\medskip

\PB{$\mu$ n'est pas, dans ce paragraphe, le poids moléculaire moyen.}

\medskip

L'expérience des calculs a montré que le formalisme d'Eggleton
\etal{eff73}, même avec la modification de Proffit \& Michaud \cite{pm91}
formule (4), donnait des taux d'ionisation moyens avec des paliers
et des gradients beaucoup plus marqués que ceux obtenus avec une meilleure
prise en compte de l'ionisation de pression. Ainsi, au centre du Soleil,
le formalisme d'Eggleton \etal{eff73} donne un taux d'ionisation de 100\% pour
Fe XXV alors que le calcul plus précis de Gabriel \cite{g97} ne donne que
85\%.

Pour tenir compte de l'ionisation de pression et d'éviter la recombinaison,
la routine {\tt saha} de \cesam2k utilise une correction
numérique dérivée
de l'analyse de Clayton \cite{c68} p. 140-145. Le potentiel d'ionisation est
réduit lorsque la distance moyenne entre les ions devient de l'ordre de
grandeur de la longueur de Debye:
\[R_{\rm D}=\sqrt{\frac{kT}{4\pi e^2\rho
N_0\zeta}},\quad\zeta=\sum_{i\not=e}\bar{Z_i}(\bar{Z_i}+1)x_i\]
et se réduit\footnote{$je$ est la charge de l'ion.} à:
\[\chi^\prime_{j,i}=\max(0,\chi_{i,j}-\chi_{i,j}^{\rm C}),
\quad\chi_{i,j}^{\rm C}=\frac{je^2}{R_{\rm D}}.\]
Cette correction appliquée seulement sur le potentiel d'ionisation s'est
avérée insuf\-fisante pour éviter la recombinaison à haute
température,
en particulier, de l'ion HII. La routine {\tt saha} utilise un ajustement
{\em numérique} agissant à la fois sur le potentiel d'ionisation et le
quotient des fonctions de partition celles-ci devant être tronquées
de façon cohérente
avec la diminution des potentiels d'ionisation (\cf Clayton \cite{c68} figure
2-16, p. 146) ce qui n'est pas possible dans la présente situation puisque
les fonctions de partition sont limitées aux poids statistiques des niveaux
fondamentaux. L'artifice consiste en une diminution progressive de 1 à 0
%{\em à la fois}, du potentiel d'ionisation et 
du rapport des poids
statistiques des niveaux fondamentaux $g_{j-1,i}/g_{j,i}$, dès que la
quantité:
\[x=\frac{\chi_{i,j}}{\chi^\prime_{j,i}}-1\]
devient inférieure à une valeur $p$ fixée, pour le
niveau d'ionisation $j=1,\ldots,Z_i$ de l'élément chimique $i$.
L'expérience des calculs a
montré que la valeur $p=4$ permettait de retrouver, pour des modèles
d'étoiles de la séquence principale de $1M_\odot$ et
$1.4M_\odot$ respectivement, des taux d'ionisation voisins de ceux obtenus avec un calcul plus
exact, mais beaucoup plus lourds, de la correction de pression d'ionisation.
La fonction de raccordement $f(x)$ utilisée est un morceau de cubique sur
$[0,p]$, de pente nulle en $x=0$ et $x=p$:
\[f(x)=0\ {\rm si}\ x\leq 0,\ \left(\frac xp\right)^2
\left[-2\left(\frac xp\right)+3\right]\ {\rm si}\ x\in[0,p],\ f(x)=1,\ {\rm si}
\ x\geq p.\]
Avec ce formalisme approximatif, l'équation de Saha \ref{eq:saha0} devient~: 
\begin{equation}\label{eq:saha1}
\frac{n_{j-1,i}}{n_{j,i}}=\frac{g_{j-1,i}}{g_{j,i}}f(x)\exp\left(\eta
+\frac{\chi_{j,i}}{kT}	%f(x)
\right)\equiv\varphi_{j,i},\ j=1,\ldots,Z_i.
\end{equation}
On obtient ainsi l'ionisation totale du niveau $j-1$ lorsque  
$g_{j-1,i}/g_{j,i}f(x)=0$ \ie\ dès que 
$\chi_{i,j}\leq \chi_{i,j}^{\rm C}$.

\medskip

Le nombre de particules $n_i,\ i\not= e$ est:
\begin{eqnarray}
n_i&=&\sum_{j=0}^{Z_i}n_{j,i}=n_{Z_i,i}\left(1+\sum_{j=0}^{Z_i-1}
\frac{n_{j,i}}{n_{Z_i,i}}\right)=n_{Z_i,i}\left(1+\sum_{j=0}^{Z_i-1}
\prod_{k=j+1}^{Z_i}\varphi_{k,i}\right)\nonumber \\
&=&n_{Z_i,i}\left(1+\frac{n_{0,i}}{n_{Z_i,i}}+\frac{n_{1,i}}{n_{Z_i,i}}
+\ \ldots\ +\frac{n_{Z_i-1,i}}{n_{Z_i,i}}\right)\nonumber \\
&=&n_{Z_i,i}\left(1+\varphi_{1,i}\varphi_{2,i}\varphi_{3,i}
\ \ldots\ \varphi_{Z_i,i}+
\varphi_{2,i}\varphi_{3,i}\ \ldots\ \varphi_{Z_i,i}+\ \ldots\ +\varphi_{Z_i,i}
\right)\nonumber \\
&=&n_{Z_i,i}\left(1+\varphi_{Z_i,i}
\left(1+\varphi_{Z_i-1,i}\left(1+\varphi_{Z_i-2,i}
\left(1+\ldots\varphi_{1,i}\right)\right)\right)\right)\nonumber
\end{eqnarray}
et comme:
\[\varphi_{j+1,i}=\frac{n_{j,i}}{n_{j+1,i}}=\frac{n_{j,i}}{n_i}
\frac{n_i}{n_{j+1,i}}=\frac{x_{j,i}}{x_{j+1,i}},\]
on a:
\begin{eqnarray}
x_{Z_i,i}\equiv\frac{n_{Z_i,i}}{n_i}&=&\frac1{\left(1+\varphi_{Z_i,i}
\left(1+\varphi_{Z_i-1,i}\left(1+\varphi_{Z_i-2,i}
\left(1+\ldots\varphi_{1,i}\right)\right)\right)\right)},\nonumber \\
x_{j,i}&=&\varphi_{j+1,i}\,x_{j+1,i},\ j=Z_i-1,\ldots,0.\nonumber
\end{eqnarray}
La charge moyenne de l'élément chimique $i$ et le nombre
d'électrons par mole et
par unité de volume s'écrivent respectivement:
\[\bar{Z_i}\equiv\sum_{j=1}^{Z_i}jx_{j,i},
\ x_e=\sum_{j\not= e}\bar{Z_i}x_i,
\ n_e=\rho N_0x_e=\rho N_0\sum_{j\not=e}\bar{Z_i}x_i.\]
On a aussi:
\[X_i=\frac{n_i\nu_im_u}{\sum_jn_j\nu_jm_u}=\frac{n_i\nu_i}{\rho N_0},
\ x_i=\frac{X_i}{\nu_i}=\frac{n_i}{\rho N_0}.\]


\subsection{L'équation de diffusion}
Pour la particule $i$ l'équation de diffusion s'écrit (Iben \& Mc Donald
\cite{id85}):
\[\dsdx{n_i}t=-\frac 1{r^2}\dsdx{ }r(r^2n_iw_i)
+\left(\dsdx{n_i}t\right)_{\rm nucl.}\]
et les équations de Burgers \cite{b69} Eq. (18.4) et (18.5), voir aussi Cox
\etal{cgk89}, Thoul \etal{tbl94}:
\begin{eqnarray}\label{eq:bu1}
\ddx{P_i}r-\frac{\rho_i}\rho\ddx Pr &-& n_i\bar{Z_i} eE
= \sum_j K_{ij}(w_j-w_i) +
\sum_j K_{ij} z_{ij}\frac{m_jr_i-m_ir_j}{m_i+m_j} \\
\frac52 n_ik\ddx Tr = &-& \frac52 \sum_j K_{ij} z_{ij}
\frac{m_j(w_j-w_i)}{m_i+m_j}-\frac25K_{ii}z_{ii}^{\prime\prime}r_i \nonumber \\
&-&r_i\sum_{j\not= i}K_{ij}\frac{3m_i^2+m_j^2z_{ij}^\prime\label{eq:bu2}
+\frac45m_im_jz_{ij}^{\prime\prime}}{(m_i+m_j)^2} \\
&+&\sum_{j\not= i}K_{ij}\frac{m_im_j(3+z_{ij}^\prime
-\frac45z_{ij}^{\prime\prime})}{(m_i+m_j)^2} r_j.\nonumber
\end{eqnarray}
Les équations de conservation dynamique et statique de la masse et de la
charge ont respectivement pour expression:
\begin{equation}\label{eq:cons}
\sum_i\nu_in_iw_i=0,\ e\sum_i\bar{Z_i}n_iw_i=0,\ \sum_i\bar{Z_i}n_i=0,
\ \sum_iX_i=\sum_i\nu_ix_i=1.
\end{equation} 
La densité partielle $\rho_i$ de la particule $i$
est définie par:
\begin{equation}\label{eq:roi}
\rho_i=\rho X_i=\rho x_i\nu_i=n_i\nu_im_u,\ \sum_i \rho_i=\rho
=\sum_in_i\nu_im_u,
\end{equation} 
et:
\[n_i=\rho N_0\frac{X_i}{\nu_i}=\rho N_0 x_i.\]
Le second terme du membre de gauche de l'équation (\ref{eq:bu1}) correspond
à la
force externe qui s'exerce sur la particule $i$, à savoir, la gravité
corrigée de l'accélération radiative.
\subsection{Intégrales de collision}\index{intégrales de collision}
Les coefficients de résistance du flux de chaleur
$z_{ij},z_{ij}^\prime, z_{ij}^{\prime\prime}$ sont définis à
partir des intégrales de collision $\Omega^{(kl)}_{ij}$,
 \cf Paquette \etal{ppfm86} Eq. (23) -- (25):
\[z_{ij}=1-\frac25\frac{\Omega^{(12)}_{ij}}{\Omega^{(11)}_{ij}},
\ z_{ij}^\prime=2.5-\frac25\frac{5\Omega^{(12)}_{ij}-\Omega^{(13)}_{ij}}
{\Omega^{(11)}_{ij}},
\ z_{ij}^{\prime\prime}=\frac{\Omega^{(22)}_{ij}}{\Omega^{(11)}_{ij}},\]
qui sont déterminées à partir de tabulations.
Paquette \etal{ppfm86} ont écrit $\Omega^{(kl)}_{ij}$
sous la forme, Eq. (65):
\[
\Omega^{(kl)}_{ij}\equiv F^{(kl)}_{ij}\epsilon_{ij},
\ \epsilon_{ij}\equiv\frac{e^4}4\sqrt{\frac\pi{2k^3}}
\ \frac{\bar{Z_i}^2\bar{Z_j}^2}
{\sqrt{m_{ij}}}\ \frac1{\sqrt{T^3}}.
\]
Pour $(kl)=(11),\ (12),\ (13),\ (22)$ et selon
qu'il s'agit d'un potentiel attractif, \ie\ électron/ion,
ou répulsif, \ie\ ion/ion ou électron/électron,
Paquette \etal{ppfm86} ont calculé et tabulé
les coefficients d'interpolation par splines naturelles de
$\ln F^{(kl)}_{ij}$ en fonction de la variable
$\psi_{ij}\equiv \ln\left[\ln\left(1+\gamma_{ij}^2\right)\right]$ avec:
\[\gamma_{ij}\equiv \frac{4kT\lambda}{\bar{Z_i}\bar{Z_j}e^2},
\ \lambda\equiv\max\left\{\left(\frac{kT}{4\pi e^2\sum_in_i\bar{Z_i}} \right)^{\frac12},
\left(\frac3{4\pi\sum_{\rm ions}n_i} \right)^{\frac13}\right\}.
\]

%Pour simplifer,
%Cox \etal{cgk89} et Thoul \etal{tbl94} utilisent les approximations
%de Muchmore \cite{m84}:
%\[z_{ij}=0.6,\ z_{ij}^\prime=1.3,\ z_{ij}^{\prime\prime}=2.0\]

Les quantités $K_{ij}$ sont les coefficients de résistance, ils
représentent les effets des collisions entre les particules
de types $i$ et $j$, ils ont pour expression
\cf Michaud \& Proffitt \cite{mp93}:
\begin{equation}\label{eq:kij}
K_{ij}=\frac{16}3 n_i n_j m_{ij}\Omega^{(11)}_{ij}.
\end{equation}
%Ces auteurs donnent une forme simplifiée pour
%$\Omega^{(11)}_{ij}$ valable pour la séquence principale:
%\[\Omega^{(11)}_{ij}=\epsilon_{ij}C_{ij},\ C_{ij}=
%\frac{\ln\left[1+\exp(1.2\ln\Lambda_{ij})\right]}{1.2},\]
%\[\ln\Lambda_{ij}=-19.26-\ln \bar{Z_i}\bar{Z_j}
%-\frac12\ln\rho-\frac12\ln\left[1+\left(\frac{X+1}2\right)\right]+\frac32\ln T,\]
%où $X$ est l'abondance en masse de l'hydrogène.

\subsection{Equations de Burgers pour les $x_i$ et $\mu$.}
En remplaçant la charge de chaque ion par sa charge moyenne tenant compte
de son taux d'ionisation, le poids moléculaire moyen est:
\[\mu^{-1}=\sum_{j\not= e}(1+\bar{Z_j})\frac{X_j}{\nu_j}=
\sum_{j\not= e}(1+\bar{Z_j})x_j.\] Pour les dérivations on utilisera:
\[d\mu=-\mu^2\sum_{j\not= e}(1+\bar{Z_j})dx_j,
\ \dsdx\mu{x_k}=-\mu^2(1+\bar{Z_k}),\ k\not=e.\]
Avec l'hypothèse du gaz parfait, pour chaque espèce $i$, les pressions
totale et partielles s'écrivent:
\[P=\frac{\rho{\cal R}T}\mu+P_{\rm R}=\sum_i P_i,\]
\[ P_i=n_ikT+P_{i\rm R}=
\rho x_i N_0kT+P_{i\rm R}
=\frac{P\mu}{{\cal R}T}x_iN_0kT+P_{i\rm R}=P\mu x_i+P_{i\rm R},\]
puisque la constante des gaz parfaits vérifie ${\cal R}\equiv kN_0$. On a donc:
\[n_ik=\frac PT\mu x_i.\]
%On introduit une correction $\beta$ qui permet, dans une certaine mesure, de
%tenir compte des effets non idéaux et de la pression de radiation:
%\[\beta=\frac P{\cal P}\leq 1\]
%c'est le rapport entre la pression "gaz parfait" et celle donnée
%par l'équation d'état $\cal P$.
%[{\sl Il doit être possible de tenir compte ici de la pression
%de radiation.}]
En posant, pour les électrons:
%\[x_e\equiv\frac23\frac{F_{\frac32}(\eta)}{F_{\frac12}(\eta)}
%\frac{n_e}{\rho N_0}=\frac23\frac{F_{\frac32}(\eta)}{F_{\frac12}(\eta)}
%\frac{\sum_{i\not= e}\bar{Z_i}n_i}{\rho N_0}=
%\frac23\frac{F_{\frac32}(\eta)}{F_{\frac12}(\eta)}\sum_{i\not= e}
%\bar{Z_i}x_i,\ \dsdx{x_e}{x_k}=\frac23\frac{F_{\frac32}(\eta)}
%{F_{\frac12}(\eta)}\bar{Z_k},\ k\not= e,\]

\[x_e\equiv\frac{n_e}{\rho N_0}=\frac{\sum_{i\not= e}\bar{Z_i}n_i}{\rho N_0}=
\sum_{i\not= e}\bar{Z_i}x_i,\ \dsdx{x_e}{x_k}=\bar{Z_k},\ k\not= e,\]
%en utilisant la valeur de $n_e$ tirée de Eq.~(\ref{eq:ne}),
comme $P_{e\rm R}=0$, on retrouve bien la pression électronique~: 
%avec dégénérescence (Clayton \cite{c68}, formule 2-57)~: 
%\begin{eqnarray}
%P\mu x_e&=&\frac{\rho{\cal R}T}\mu\mu\frac23\frac{F_{\frac32}(\eta)}
%{F_{\frac12}(\eta)}\frac{n_e}{\rho N_0}=
%kT\frac23\frac{F_{\frac32}(\eta)}{F_{\frac12}(\eta)}F_{\frac12}(\eta)
% \frac{4\pi kT}{3h^3}(2m_ekT)^{\frac32}\nonumber \\
%&=&\frac{8\pi kT}{3h^3}(2m_ekT)^{\frac32}F_{\frac32(\eta)}=P_e.\nonumber
%\end{eqnarray}
\[P\mu x_e=\frac{\rho{\cal R}T}\mu\mu \frac{n_e}{\rho N_0}=n_ekT=P_e.\]
L'hypothèse d'équilibre quasi-statique et $\rho_i=\rho\nu_i x_i$
(Eq. \ref{eq:roi}) entrainent:
\[\frac{\rho_i}\rho\ddx Pr=x_i\nu_i\rho g_{\rm t},\]
et:
\[n_ik\ddx Tr=\frac PT\mu x_i\frac TP\dsdx{\ln T}{\ln P}
\ddx{P}r=\mu x_i\nabla\rho g_{\rm t},\]
d'où:
\begin{eqnarray}
\ddx{P_i}r&=&\ddx{ }rP\mu x_i+\ddx{P_{i\rm R}}r=\mu x_i\ddx Pr+
\ddx{P_{i\rm R}}r+\mu P\ddx{x_i}r+P x_i\ddx{\mu}r
\nonumber \\
&=&x_i\rho\mu g_{\rm t}+\rho_i g_{{\rm R}i}+\mu P\ddx{x_i}r-P x_i\mu^2\sum_{j\not= e}
(1+\bar{Z_j})\ddx{x_j}r\nonumber \\
&=&x_i\rho(\mu g_{\rm t}+\nu_i g_{{\rm R}i})+P\mu\left(\ddx{x_i}r-x_i\mu\sum_{j\not= e}
(1+\bar{Z_j})\ddx{x_j}r\right)\nonumber
\end{eqnarray}
%Remarquons que le point essentiel de la théorie de Burgers est que
%l'équilibre quasi-statique est {\em globalement}
%réalisé {\em et non pas} pour chaque espèce prise individuellement:
%\[\ddx{P_i}r\not=\rho_i(g_{\rm e}+g_{{\rm R}i})\]
On a par ailleurs:
\[n_i\bar{Z_i}eE=\rho N_0 x_i\bar{Z_i}eE,
%\  n_ik\ddx Tr=\rho N_0 x_i k\frac TP \nabla \rho g_{\rm e}=
%\mu x_i\nabla\rho g_{\rm e},\]
%\[
\ \sum_i\nu_ix_iw_i=0,\ \sum_i\bar{Z_i}x_iw_i=0.\]
On introduit les relations précédentes dans les
équations de Burgers~\eq{eq:bu1},~\eq{eq:bu2}
et de conservation~\eq{eq:cons}.
Les coefficients de résistance $K_{ij}$ sont de la forme:
\[K_{ij}=\frac{16}3 n_i n_j m_{ij}\Omega^{(11)}_{ij}=\frac{16}3(\rho N_0)^2
x_ix_jm_{ij}\Omega^{(11)}_{ij}=k_{ij}x_ix_j\rho^2,
\ k_{ij}\equiv\frac{16N_0^2}3m_{ij}\Omega_{ij}^{(11)}.\]
En fonction des $x_i$ et de $\mu$ les équations de Burgers s'écrivent:
\begin{eqnarray}
&&x_i\rho\left(\mu g_{\rm t}-\nu_i\left(g_{\rm t}-g_{{\rm R}i}\right)\right)+P\mu\left(\ddx{x_i}r-
x_i\mu\sum_{j\not=e}(1+\bar{Z_j})\ddx{x_j}r\right)-\rho N_0 x_i\bar{Z_i}eE
\nonumber \\
&&=\rho^2\left\{x_i\sum_j k_{ij}x_j(w_j-w_i) -
x_i\sum_j k_{ij} x_j z_{ij}\frac{m_ir_j-m_jr_i}{m_i+m_j}\right\},\ i=1,\ldots,c,
\nonumber\\
&&-x_i\rho\nabla\mu g_{\rm t}=\rho^2\left\{x_i\sum_j k_{ij} x_jz_{ij}
\frac{m_j(w_j-w_i)}{m_i+m_j}+x_i\frac4{25}k_{ii}x_iz_{ii}^{\prime\prime}r_i
\right.  \nonumber \\
&&+\frac25r_ix_i\sum_{j\not= i}k_{ij}x_j\frac{3m_i^2+m_j^2z_{ij}^\prime
+\frac45m_im_jz_{ij}^{\prime\prime}}{(m_i+m_j)^2}\nonumber  \\
&&\left. -\frac25x_i\sum_{j\not= i}k_{ij}x_jr_j\frac{m_im_j(3+z_{ij}^\prime
-\frac45z_{ij}^{\prime\prime})}{(m_i+m_j)^2}\right\},\ i=1,\ldots,c,e, \nonumber
\end{eqnarray}
ou encore, en posant:
\[\nu_{i\rm R}\equiv\nu_i\left(1-\frac{g_{{\rm R}i}}{g_{\rm t}}\right),
\ \dsdx{\nu_{i\rm R}}{x_i}=-\frac{\nu_i}{g_{\rm t}}
\left(\dsdx{g_{{\rm R}i}}{x_j}-\frac{g_{{\rm R}i}}{g_{\rm t}}\dsdx{g_{\rm t}}{x_j}
\right)
\]
puis arrangements et division de chaque membre par $x_i\rho^2$:
\begin{eqnarray}
\frac{g_{\rm t}}\rho(\mu-\nu_{i\rm R}) &+&\frac{{\cal R}T}\rho
\left(\frac1{x_i}\ddx{x_i}r-\mu\sum_{j\not= e}(1+\bar{Z_j})
\ddx{x_j}r\right)=\sum_{j\not= i} k_{ij}x_j(w_j-w_i) \nonumber \\
&-&\sum_{j\not= i} k_{ij} x_j z_{ij}\frac{m_ir_j-m_jr_i}{m_i+m_j}
+\frac{\bar{Z_i}eN_0E}\rho,\ i=1,\ldots,c,
\label{eq:BU1}\\
-\frac{\mu\nabla g_{\rm t}}\rho&=&\sum_{j\not= i} k_{ij} x_jz_{ij}
\frac{m_j(w_j-w_i)}{m_i+m_j}+  \nonumber \\
&+&r_i\left(\frac4{25}k_{ii}x_iz_{ii}^{\prime\prime}+\frac25
\sum_{j\not= i}k_{ij}x_j\frac{3m_i^2+m_j^2z_{ij}^\prime
+\frac45m_im_jz_{ij}^{\prime\prime}}{(m_i+m_j)^2}\right)-\nonumber  \\
&-&\frac25\sum_{j\not= i}k_{ij}x_jr_j\frac{m_im_j(3+z_{ij}^\prime
-\frac45z_{ij}^{\prime\prime})}{(m_i+m_j)^2},\ i=1,\ldots,c,e.
\label{eq:BU2}
\end{eqnarray}
Les données sont la température $T$, la densité $\rho$,
les gravités $g_{\rm t}$ et $g_{{\rm R}i}$, le gradient $\nabla$, les abondances
des ions $x_i,\ i=1,\ldots,c$. Des abondances on déduit
le nombre d'électrons et le poids moléculaire moyen.

\medskip

\PB{Faute de mieux, le poids moléculaire moyen, les pressions partielles,
ne peuvent pas être cohérents avec ceux utilisés dans
l'équation d'état retenue par \cesam2k.}

\medskip

\noindent Les inconnues sont les $c+1$ vitesses de diffusion $w_1,\ldots,w_c,w_e$,
les $c+1$ flux de chaleur $r_1,\ldots,r_c,r_e$ et le champ électrique $E$,
soit $2c+3$ inconnues. On ajoute aux $2c+1$
\eq{eq:BU1} et~\eq{eq:BU2}, les deux équations de conservation:
\begin{equation}\label{eq:BU3}
\sum_{i=1}^c\nu_ix_iw_i=0,\ \sum_{i=1}^e\bar{Z_i}x_iw_i=0.
\end{equation}
Les~\eq{eq:BU1},~\eq{eq:BU2} écrites pour les $x_i$ et $\mu$ et
\eq{eq:BU3} forment le système linéaire:
\[A\omega=\gamma+G\,D_x,\ G\equiv\frac{{\cal R}T}\rho\Delta,\]
on a noté:
\[\omega\equiv(w_1,\ldots,w_c,w_e,r_1,\ldots,r_c,r_e,E)^T,\]
\[D_x\equiv(\ddx{x_1}r,\ldots,\ddx{x_c}r,0,\ldots,0)^T,\]
\[\gamma\equiv\frac{g_{\rm t}}\rho\left(\mu-\nu_{1\rm R},\ldots,
\mu-\nu_{c\rm R},-\mu\nabla,\ldots,-\mu\nabla,0,0\right)^T.\]
\[\dsdx{}{x_j}\left(\frac{\mu g_{\rm t}\nabla}\rho\right)=\frac{g_{\rm t}\nabla}\rho
\dsdx\mu{x_j}+\frac{\mu\nabla}\rho\dsdx{g_{\rm t}}{x_j}+
\left\{\frac{\mu g_{\rm t}}\rho\dsdx\nabla{x_1}-
\frac{\mu g_{\rm t}\nabla}{\rho^2}\dsdx\rho{x_1}\right\}_{j=1},\]
\[\dsdx{}{x_j}\left(\frac{g_{\rm t}(\mu-\nu_{i\rm R})}\rho\right)=
\frac{(\mu-\nu_{i\rm R})}\rho\dsdx{{g_{\rm t}}}{x_j}+\frac{g_{\rm t}}\rho
\left(\dsdx\mu{x_j}-\dsdx{\nu_{i\rm R}}{x_j}\right)+
\left\{\frac{g_{\rm t}}{\rho^2}(\mu-\nu_{i\rm R})\dsdx\rho{x_1}\right\}_{j=1}.\]


\begin{eqnarray}
\Delta&\equiv&{\rm\,Diag}\left(\frac1{x_1},\frac1{x_2},\ldots,\frac1{x_c}%,\frac1{x_e}
,0,\ldots,0\right)-\mu{\cal M},\nonumber\\
{\cal M}&\equiv&\left(\begin{array}{cccccccccc}
(1+\bar Z_1)&(1+\bar Z_2)&(1+\bar Z_3)&\ldots&(1+\bar Z_c)&0&0 &\ldots&0&0\\
(1+\bar Z_1)&(1+\bar Z_2)&(1+\bar Z_3)&\ldots&(1+\bar Z_c)&0&0 &\ldots&0 &0\\
\ddots  &\ddots&\ddots &\ddots&\ddots  &\ddots&\ddots&\ddots&\ddots&\ddots\\
(1+\bar Z_1)&(1+\bar Z_2)&(1+\bar Z_3)&\ldots&(1+\bar Z_c)&0&0 &\ldots&0 &0\\
(1+\bar Z_1)&(1+\bar Z_2)&(1+\bar Z_3)&\ldots&(1+\bar Z_c)&0&0&\ldots&0&0\\
0&0&0&\ldots&0&0&0&\ldots&0&0\\
\ddots  &\ddots&\ddots &\ddots&\ddots  &\ddots&\ddots&\ddots&\ddots&\ddots\\
0&0&0&\ldots&0&0&0&\ldots&0&0\\
\end{array}\right).\nonumber
\end{eqnarray}
La structure de la matrice $A$, de dimension $[2c+3\times2c+3]$,
est la suivante:
\[\left(\begin{array}{ccc}
A_w & A_r& A_e\\
A_c&0&0
\end{array}\right),\]
les dimensions des sous matrices  sont respectivement $A_w[2c+1\times c+1]$,
$A_r[2c+1\times c+1]$, $A_e[2c+1\times 1]$ et $A_c[2\times c+1]$. En notant:
\[ p_{ij}\equiv k_{ij}x_j,
\ q_{ij}\equiv\frac{k_{ij}x_jz_{ij}m_j}{m_i+m_j}=
p_{ij}\frac{z_{ij}m_j}{m_i+m_j}, \]
\[A_w\equiv\left(\begin{array}{cccccccc}
-\sum_{j\not= 1} p_{1j} & p_{12} & p_{13} &\ldots&p_{1c} &p_{1e} \\
p_{21} &-\sum_{j\not= 2} p_{2j}& p_{23} &\ldots&p_{2c} &p_{2e} \\
\ddots&\ddots&\ddots&\ddots&\ddots&\ddots\\
p_{c1} &p_{c2} & p_{c3} &\ldots&-\sum_{j\not= c} p_{cj}&p_{ce} \\
%p_{e1} &p_{e2} & p_{e3} &\ldots&p_{ec} &-\sum_{j\not= e} p_{ej} \\
-\sum_{j\not= 1} q_{1j} & q_{12} & q_{13} &\ldots&q_{1c} &q_{1e}\\
q_{21} &-\sum_{j\not= 2} q_{2j}  &q_{23} &\ldots&q_{2c} &q_{2e}\\
\ddots&\ddots&\ddots&\ddots&\ddots&\ddots\\
q_{c1} &q_{c2} & q_{c3} &\ldots&-\sum_{j\not= c} q_{cj}&q_{ce}\\
%q_{e1} &q_{e2} & q_{e3} &\ldots&q_{ec} &-\sum_{j\not= e} q_{ej}
\end{array}\right),\]
avec les notations:
\begin{eqnarray}
\bar{q_{ij}}&\equiv&-\frac{k_{ij}x_jz_{ij}m_i}{m_i+m_j}=
-q_{ij}\frac{m_i}{m_j},\nonumber \\
d_i&\equiv&\frac4{25}k_{ii}x_iz_{ii}^{\prime\prime}+\frac25
\sum_{j\not= i}k_{ij}x_j\frac{3m_i^2+m_j^2z_{ij}^\prime
+\frac45m_im_jz_{ij}^{\prime\prime}}{(m_i+m_j)^2},\nonumber \\
f_{ij}&\equiv&-\frac25k_{ij}x_j\frac{m_im_j(3+z_{ij}^\prime
-\frac45z_{ij}^{\prime\prime})}{(m_i+m_j)^2},\nonumber
\end{eqnarray}
\[
A_r\equiv\left(\begin{array}{cccccccc}
\sum_{j\not= 1} q_{1j}&\bar{q_{12}}&\bar{q_{13}}&\ldots&
\bar{q_{1c}}&\bar{q_{1e}} \\
\bar{q_{21}}&\sum_{j\not= 2}q_{2j}&\bar{q_{23}}&\ldots&
\bar{q_{2c}}&\bar{q_{2e}} \\
\ddots&\ddots&\ddots&\ddots&\ddots&\ddots\\
\bar{q_{c1}}&\bar{q_{c2}}&\bar{q_{c3}}&\ldots&\sum_{j\not= c}
 q_{cj}&\bar{q_{ce}}\\
%\bar{q_{e1}}&\bar{q_{e2}}&\bar{q_{e3}}&\ldots&
%\bar{q_{ec}}&\sum_{j\not= e}q_{ej}\\
d_1 &f_{12}&f_{13}&\ldots&f_{1c}&f_{1e}\\
f_{21}&d_2   &f_{23}&\ldots&f_{2c}&f_{2e}\\
\ddots&\ddots&\ddots&\ddots&\ddots&\ddots\\
f_{c1}&f_{c2}&f_{c3}&\ldots&d_c   &f_{ce}\\
f_{e1}&f_{e2}&f_{e3}&\ldots&f_{ec}&d_e
\end{array}\right),\]
\[A_e\equiv\frac{N_0e}\rho (\bar Z_1,\ldots,\bar Z_c,0,\ldots,0)^T,\]
\[A_c\equiv\left(\begin{array}{cccccccc}
\nu_1x_1&\nu_2x_2&\ldots&\nu_cx_c&0&0&\ldots&0 \\
\bar Z_1x_1&\bar Z_2x_2&\ldots&\bar Z_cx_c&\bar Z_ex_e&0&\ldots&0 \\
\end{array}\right).\]
La solution s'écrit:
\[\omega=A^{-1}\gamma+A^{-1}G\, D_x=v+B\, D_x\Rightarrow,\]
ou encore, pour les ions, $i=1,\ldots,c$:
\[w_i=v_i+\sum_{j=1}^cb_{ij}\ddx{x_j}r,\ i=1,\ldots,c,\]
les $v_i$ sont les $c$ premiers coefficients de $v=A^{-1}\gamma$,
les $b_{ij},\ (i,j=1,\ldots,c)$, sont les coefficients des $c$
premières lignes et colonnes de $B=A^{-1}G$.
Pour l'ion $i$ l'équation de diffusion devient:
\begin{eqnarray}
\dsdx{n_i}t&=&-\frac 1{r^2}\dsdx{ }r(r^2n_iw_i)+
\left(\dsdx{n_i}t\right)_{\rm nucl.}=\nonumber \\
\dsdx{\rho N_0 x_i}t&=&-\frac1{r^2}
\dsdx{ }r\left\{ N_0r^2 \rho x_i\left(\sum_{j=1}^cb_{ij}\dsdx{x_j}r+v_i\right)
\right\}+\left(\dsdx{\rho N_0 x_i}t\right)_{\rm nucl.}\nonumber
\end{eqnarray}
ou encore, en notant $\dsdt{\star}\equiv\dsdx\star t+{\cal V}\bullet\nabla\star$
l'opérateur de dérivation
lagrangien "en suivant le mouvement":
\[\dsdt{x_i}=-\frac{4\pi}{4\pi r^2\rho}
\dsdx{ }r\left\{ \frac{4\pi r^4\rho^2}{4\pi r^2\rho}
 x_i\left(\sum_{j=1}^cb_{ij}\dsdx{x_j}r+v_i\right)\right\}
+\left(\dsdt{x_i}\right)_{\rm nucl.}\ \Longrightarrow\]
\[\dsdt{x_i}=-\dsdx{ }m\left\{ 4\pi r^2 \rho \left(4\pi r^2 \rho x_i
\sum_{j=1}^cb_{ij}\dsdx{x_j}m +v_ix_i\right)\right\}
+\left(\dsdt{x_i}\right)_{\rm nucl.}.\]
On a utilisé l'équation de continuité:
\[0=\dsdx\rho t+\nabla\bullet\rho {\cal V}=\dsdx\rho t+\rho\nabla\bullet
{\cal V}+
{\cal V}\bullet\nabla\rho=x_i\dsdx\rho t+x_i\rho\nabla\bullet
{\cal V}+x_i{\cal V}\bullet\nabla\rho,\]
et la conservation des particules d'espèce $i$:
\[\nabla\bullet x_i\rho {\cal V}=0\Rightarrow\rho {\cal V} \bullet\nabla x_i=
-x_i\rho\nabla\bullet {\cal V}-x_i {\cal V}\bullet\nabla\rho,\]
\[\Rightarrow \dsdx{\rho x_i}t=x_i\dsdx\rho t+\rho\dsdx{x_i}t
-x_i\dsdx\rho t-x_i\rho\nabla\bullet {\cal V}-x_i {\cal V}\bullet\nabla\rho
=\rho\left(\dsdx{x_i}t+{\cal V} \bullet\nabla x_i\right)=\rho\dsdt{x_i}.\]
En posant:
\[D_i\equiv\left(d_{i1},\ldots,d_{ic},\right)
\equiv-\left(x_ib_{i1},\ldots,x_ib_{ic},\right),\ i=1,\ldots,c,\]
on obtient la forme de l'équation de diffusion utilisée dans
\cesam2k\footnote{Avec la variable d'espace Lagrangienne $m$ les opérateurs
$\dsdt{ }$ et $\dsdx{ }t$ sont équivalents.}:
\[\dsdt{x_i} = -\dsdx {F_i}m + \Psi_i,\ F_i=-4\pi R^2\rho\left(4\pi R^2\rho
D_i\bullet\nabla\!_{\rm m}{\cal X}-v_i x_i\right)\]
où $\nabla\!_{\rm m}{\cal X}$ est le vecteur:
\[\left( \dsdx{x_1}m,\ldots,\dsdx{x_{\rm c}}m\right)^{\rm T}.\]

\subsection{Jacobien.}
\cesam2k intégrant l'équation de diffusion avec un schéma implicite,
il est nécessaire de calculer les dérivées:
\[ \dsdx{v_i}{x_k},\ \dsdx{d_{ij}}{x_k}=-x_i\dsdx{b_{ij}}{x_k}-b_{ij}\delta_{ik}
,\ i,j,k=1,\ldots,c,\ \delta_{ik}=1\ {\rm si}\ i=k,\ 0\ {\rm sinon}.\]
En dérivant:
\[A^{-1}A=I\Rightarrow\dsdx{A^{-1}}{x_k}A+A^{-1}\dsdx A{x_k}=0\Rightarrow
\dsdx{A^{-1}}{x_k}=-A^{-1}\dsdx A{x_k}A^{-1},\]
d'où:
\[\dsdx{A^{-1}\gamma}{x_k}=\dsdx{A^{-1}}{x_k}\gamma+A^{-1}\dsdx\gamma{x_k}=
-A^{-1}\dsdx A{x_k}A^{-1}\gamma+A^{-1}\dsdx\gamma{x_k}=
A^{-1}\left(\dsdx\gamma{x_k}-\dsdx A{x_k}v\right),\]
et de façon similaire:
\[\dsdx B{x_k}=A^{-1}\left(\dsdx G{x_k}-\dsdx A{x_k}B\right).\]

\subsection{Equation de diffusion du moment angulaire}\label{sec:eqdima}
L'équation de diffusion du moment angulaire par unité de masse
$\frac23r^2\Omega$ s'écrit:\index{moment angulaire}
\[
\dsdt{r^2\Omega}=\frac\partial{\partial m}
\left(16\pi^2r^4\rho^2{\cal D}_\Omega
\dsdx{r^2\Omega}m\right)
\]
où ${\cal D}_\Omega$ est le coefficient de diffusion
du moment angulaire. En tenant compte de l'hypothèse de contraction
quasi-statique \ie\ $\dsdx rt~\sim0$, on obtient:
\[
r^2\dsdx\Omega t=\frac\partial{\partial m}\left[16 \pi^2 r^4\rho^2
{\cal D}_\Omega\left(2r\Omega\dsdx rm+r^2\dsdx\Omega m\right)\right]=
-\dsdx{F_\Omega}m-2r\Omega\dsdx rt.
\]
avec:
\[
F_\Omega\equiv -8\pi r^3\rho{\cal D}_\Omega\left(\Omega+2\pi r^3\rho
\dsdx\Omega m\right)
\]
équation de même structure que l'équation de diffusion pour les
espèces chimiques~\eq{eq:diff}.
