text:
# Density and temperature dependence of nucleon-nucleon elastic cross section

## Abstract

The in-medium neutron-proton, proton-proton (neutron-neutron) elastic scattering cross sections ($`\sigma _{np}^{}`$, $`\sigma _{pp(nn)}^{}`$) are studied based on the effective Lagrangian of density dependent relativistic hadron theory in which the $`\delta `$\[$`a_0(980)`$\] meson is included. Our study shows that at low densities the $`\sigma _{np}^{}`$ is about $`34`$ times larger than $`\sigma _{pp(nn)}^{}`$ and at densities higher than the normal density the isospin effect is almost washed out. Because of coupling to $`\delta `$ meson the $`\sigma _{nn}^{}`$ and $`\sigma _{pp}^{}`$ are different in isospin asymmetric medium following the splitting of the proton and neutron mass. The isospin effect on the density dependence of the in-medium nucleon elastic cross section is dominantly contributed by the isovector $`\delta `$ and $`\rho `$ mesons. The temperature effect on the $`\sigma _{np}^{}`$ and $`\sigma _{pp(nn)}^{}`$ is studied. It is shown that the temperature effect is weaker compared with the density effect but it becomes obvious as density increases.

The rapid advance in nuclear reactions using rare isotopes has opened up several new frontiers in nuclear science LiB01 ; Pa97 ; Ye01 ; Ts01 ; Di01 ; So03 . In particular, the intermediate energy heavy rare isotopes currently available at National Superconducting Cyclotron Laboratory (NSCL/MSU) and the future new facilities RIA in USA and SIS-200 at GSI in Germany provide a unique opportunity to explore novel properties of dense neutron-rich matter that was not in reach in terrestrial laboratories before.

In order to study theoretically the neutron-rich nuclear collisions at intermediate energy within a microscopic transport model approach, both the isospin-dependent mean field and the two-body scattering cross sections should be introduced. The isospin dependent mean field has been studied extensively with both non-relativistic and relativistic theory. Concerning the nucleon-nucleon elastic cross sections (ECS), it is already known that up to hundreds MeV the free proton-neutron cross section is about $`23`$ time larger than that of proton-proton (neutron-neutron)’s Ch68 . However, it is not clear how the nuclear medium corrects the neutron-proton and proton-proton (or neutron-neutron) ECS. The isospin dependent in-medium ECS were calculated based on QHD-II model in which the isospin-vector channel was introduced through a coupling to the vector-isovector $`\rho `$ meson field Li00 . In this work we intend to introduce the scalar-isovector $`\delta [a_0(980)]`$ meson into the effective Lagrangian to calculate the in-medium two-body scattering elastic cross sections. In free NN potentials the $`\delta `$ meson is a standard and essential ingredient. The importance of the $`\delta `$ meson in nuclear matter with extreme neutron-to-proton ratios was conjectured in Hu96 ; Kubis97 ; Hu98 . It was pointed out that conventional mean-field models, neglecting the $`\delta `$ meson, are likely to miss an important contribution to the isospin degree of freedom which ought to be necessary for a proper description of strongly asymmetric matter. Recently, the role of the $`\delta `$ meson for asymmetric nuclear matter was studied by Liu02 . Some calculation results for the structure of exotic nuclei showed the importance of the inclusion of the $`\delta `$ meson for the stability conditions of drip-line nuclei Hofmann01 ; Leja01 . Therefore it would be worthwhile to study how the inclusion of the $`\delta `$ meson changes the medium effect on nucleon-nucleon cross sections. We take the same approach as in Mao94 ; Mao942 ; Li00  for the formalism and the effective lagrangian employed in the density dependent relativistic hadron theory for asymmetric nuclear matter and exotic nuclei given in Hofmann01 .

The effective Lagrangian density is taken as

$`L`$ $`=`$ $`\overline{\mathrm{\Psi }}[i\gamma _\mu ^\mu M_N]\mathrm{\Psi }+{\displaystyle \frac{1}{2}}_\mu \sigma ^\mu \sigma {\displaystyle \frac{1}{4}}F_{\mu \nu }F^{\mu \nu }`$
$`+{\displaystyle \frac{1}{2}}_\mu 𝜹^\mu 𝜹{\displaystyle \frac{1}{4}}L_{\mu \nu }L^{\mu \nu }{\displaystyle \frac{1}{2}}m_\sigma ^2\sigma ^2+{\displaystyle \frac{1}{2}}m_\omega ^2\omega _\mu \omega ^\mu {\displaystyle \frac{1}{2}}m_\delta ^2𝜹^2+{\displaystyle \frac{1}{2}}m_\rho ^2𝝆_\mu 𝝆\text{ }^\mu `$
$`+g_\sigma \overline{\mathrm{\Psi }}\mathrm{\Psi }\sigma g_\omega \overline{\mathrm{\Psi }}\gamma _\mu \mathrm{\Psi }\omega ^\mu +g_\delta \overline{\mathrm{\Psi }}𝝉\mathrm{\Psi }𝜹{\displaystyle \frac{1}{2}}g_\rho \overline{\mathrm{\Psi }}\gamma _\mu 𝝉\mathrm{\Psi }𝝆^\mu ,`$
where $`g_\sigma `$, $`g_\omega `$, $`g_\rho `$, and $`g_\delta `$ are density dependent and $`F_{\mu \nu }_\mu \omega _\nu _\nu \omega _\mu `$, $`L_{\mu \nu }_\mu 𝝆_\nu _\nu 𝝆_\mu `$.

The distribution functions for fermions ($`f(𝐩)`$) and antifermions ($`\overline{f}(𝐩)`$) respectively are $`f_i(𝐩)=1/(1+e^{(E_i^{}(𝐩)\mu _i^{})/T})`$, $`\overline{f}_i(𝐩)=1/(1+e^{(E_i^{}(𝐩)+\mu _i^{})/T})`$, where $`i`$ represents proton ($`i=p`$) or neutron ($`i=n`$) and $`E^{}=\sqrt{𝐩^2+M_{i}^{}{}_{}{}^{2}}`$. The effective nucleon mass $`M^{}`$ is given by

$$M^{}=M_0+\mathrm{\Sigma }_{H(\sigma )}(x,\tau )+\mathrm{\Sigma }_{H(\delta ^0)}(x,\tau ).$$
(1)
$`\mathrm{\Sigma }_{H(\sigma )}`$ and $`\mathrm{\Sigma }_{H(\delta )}`$ are the self-energy parts of nucleon contributed from $`\sigma `$ and $`\delta `$ meson, respectively. Because the self-energy $`\mathrm{\Sigma }_{H(\delta ^0)}(x,\tau )`$ has opposite sign for neutron and proton for isospin asymmetric medium the correction of the nuclear medium to proton mass and neutron mass from $`\delta `$ meson is of opposite sign. Thus the proton and neutron effective mass differ for isospin asymmetric systems. The effective chemical potential $`\mu ^{}`$ is

$$\mu _i^{}=\mu _i+\mathrm{\Sigma }_{H(\omega )}^0(x,\tau )+\mathrm{\Sigma }_{H(\rho ^0)}^0(x,\tau ).$$
(2)
$`\mathrm{\Sigma }_{H(\omega )}`$ and $`\mathrm{\Sigma }_{H(\rho )}`$ are the self-energy parts of nucleon contributed from $`\omega `$ and $`\rho `$ meson, respectively. Within the mean field approximation the equation of state can be written as

$`ϵ`$ $`=`$ $`{\displaystyle \underset{i=n,p}{}}2{\displaystyle \frac{d^3𝐩}{(2\pi )^3}E_i^{}(𝐩)[f_i(𝐩)+\overline{f}_i(𝐩)]}+U(\sigma )`$ (3)
$`+{\displaystyle \frac{1}{2}}\mathrm{\Sigma }_{H(\omega )}^0(x,\tau )\rho _B+{\displaystyle \frac{1}{2}}\mathrm{\Sigma }_{H(\rho ^0)}^0(x,\tau )\rho _{B3}{\displaystyle \frac{1}{2}}\mathrm{\Sigma }_{H(\delta ^0)}(x,\tau )\rho _{S3},`$
where the baryon and scalar densities $`\rho _B`$ and $`\rho _S`$ are

$$\rho _B=\underset{i=n,p}{}2\frac{d^3𝐩}{(2\pi )^3}[f_i(𝐩)\overline{f}_i(𝐩)],$$
(4)
$$\rho _S=\underset{i=n,p}{}2\frac{d^3𝐩}{(2\pi )^3}\frac{M_i^{}}{E_i^{}}[f_i(𝐩)+\overline{f}_i(𝐩)]$$
(5)
while $`\rho _{B}^{}{}_{3}{}^{}=\rho _{B}^{}{}_{p}{}^{}\rho _{B}^{}{}_{n}{}^{}`$ and $`\rho _{S}^{}{}_{3}{}^{}=\rho _{S}^{}{}_{p}{}^{}\rho _{S}^{}{}_{n}{}^{}`$.

The functional form of the density dependent coupling for $`\sigma `$, $`\omega `$, $`\delta `$, and $`\rho `$ meson taken from Ref. Hofmann01  is $`g_\alpha (\rho )=A_\alpha (1+B_\alpha (\rho /\rho _0+D_\alpha )^2)/(1+C_\alpha (\rho /\rho _0+E_\alpha )^2)`$, where $`\alpha `$ represents $`\sigma `$, $`\omega `$, $`\delta `$, or $`\rho `$, and the parameters $`A_\alpha `$, $`B_\alpha `$, $`C_\alpha `$, $`D_\alpha `$, and $`E_\alpha `$ are listed in Table 1.

By solving Eq. 1 self-consistently we obtain the effective proton and neutron mass. We illustrate the density dependence of the neutron and proton effective mass at $`T=0`$ MeV for $`Y_p=0.3`$, and $`Y_p=0.5`$ ($`Y_p=\rho _p/(\rho _p+\rho _n)`$) in Fig. 1 (a). The proton and neutron effective mass differ at $`Y_p=0.3`$ case due to the inclusion of the $`\delta `$ meson. The proton effective mass is larger than the neutron effective mass in neutron-rich nuclear medium. With the increase of density, the splitting of the proton and neutron mass increases, which will affect the in-medium ECS.

Fig. 1 (b) shows the temperature dependence of the neutron and proton effective mass for $`Y_p=0.3`$. The ratios $`M^{}(T)/M^{}(T=0)`$ for different cases are also shown in the upper-right plot. One can see from this figure that the effective nuclear mass first increases slightly until $`T=150`$ MeV and then decreases quickly as temperature increases further, which has been also observed in Refs. Mi02 ; Li96 . The sub-figure shows that increasing slope of $`M^{}(T)/M^{}(T=0)`$ depends on density as well as the species of nucleons. The behavior of the density and temperature dependence of the effective mass should influence the density and temperature dependence of in-medium ECS which will be discussed later.

With the same approach as in Mao94 ; Mao942 ; Li00  we obtain the expressions of $`\sigma _{pp(nn)}`$ and $`\sigma _{pp(nn)}`$ with inclusion of the individual $`\delta `$ term and the crossing terms of $`\delta `$-$`\sigma `$ $`\delta `$-$`\omega `$ as well as $`\delta `$-$`\rho `$. Because of the coupling to $`\delta `$ meson the collisions of $`npnp`$ and $`nppn`$ should be treated differently, which makes the calculations of in-medium ECS much more complicated. And for the same reason, $`\sigma _{nn}^{}`$ is not equal to $`\sigma _{pp}^{}`$ for isospin asymmetric systems.

Before calculating the ECS, we have to introduce the common used form factor of nucleon-meson-nucleon vertex $`F_i(t)=\mathrm{\Lambda }_i^2/(\mathrm{\Lambda }_i^2t)`$, where the subscripts $`i`$ denote the vertices for different meson-nucleon couplings and $`\mathrm{\Lambda }_i`$ is the cut-off mass for $`i`$ meson species. In this work, the values of $`\mathrm{\Lambda }_\sigma =1200`$ MeV, $`\mathrm{\Lambda }_\omega =808`$ MeV, (which are same as those chosen in BASS98 ), $`\mathrm{\Lambda }_\delta =1000`$ MeV and $`\mathrm{\Lambda }_\rho =800`$ MeV are used. We find there is very week influence on the results when taking $`\mathrm{\Lambda }_\delta =8001000`$MeV and $`\mathrm{\Lambda }_\rho =600800`$ MeV. Again, the $`F_i(t)`$ is different for $`npnp`$ and $`nppn`$ channels when $`M_n^{}M_p^{}`$.

The main contributions to the in-medium ECS come from the $`\sigma \sigma `$, $`\omega \omega `$ and the crossing term $`\sigma \omega `$. Fig. 2 (a) shows the contributions of $`\sigma \sigma `$, $`\omega \omega `$ and $`\sigma \omega `$ terms to $`\sigma _{pp(nn)}`$ and $`\sigma _{np}`$ at $`T=0`$ MeV, $`E_K=10`$ MeV for $`Y_p=0.5`$. The ratio between $`\sigma _{np}`$ and $`\sigma _{pp(nn)}`$ is also shown in the upper-right plot. From this figure one can find both $`\sigma _{pp(nn)}`$ and $`\sigma _{np}`$ have similar shape of density dependence and the $`\sigma _{np}`$/$`\sigma _{pp(nn)}`$ ratio is almost density independent.

The contributions from the $`\delta `$ and $`\rho `$ mesons and related crossing terms are shown in Fig. 3. In order to explore how $`\delta `$ and $`\rho `$ meson field affect the in-medium ECS, we show the individual contributions from each term. The sub-figure shows the sum of the $`7`$ terms contributing to $`\sigma _{pp(nn)}`$ and $`\sigma _{np}`$. This figure shows rich information about the contributions of each term to the in-medium ECS as well as the cancellation effect among the $`7`$ terms. Firstly, the magnitude of all individual contributions to the in-medium ECS decreases with density. This density dependence results from both the density dependence of the meson-nucleon couplings and the effective masses of protons and neutrons. Among the $`7`$ terms, the crossing terms of the $`\sigma `$ meson to $`\delta `$ and $`\rho `$ mesons provide the most important contributions and the crossing terms of $`\omega `$ to $`\delta `$ and $`\rho `$ mesons are the next important terms. Further, the cancellation among the contributions from different terms are very strong so the final results are the delicate balance of all contributed terms.

In Fig. 2 (b) we show the density dependence of $`\sigma _{pp}^{}`$, $`\sigma _{nn}^{}`$, and $`\sigma _{np}^{}`$ at $`T=0`$ MeV for $`Y_p=0.3`$ and $`Y_p=0.5`$. The $`\sigma _{pp(nn)}^{}`$ and $`\sigma _{np}^{}`$ for $`Y_p=0.5`$ are shown with thick line types while those for $`Y_p=0.3`$ are illustrated with light lines. One can see from the figure that $`\sigma _{pp}^{}`$, $`\sigma _{nn}^{}`$, and $`\sigma _{np}^{}`$ decrease with the increase of densities. This behavior is consistent with other theoretical calculations, for instance, Refs. LiGQ93 ; Li96 ; LiGQ2  and the phenomenological expression of $`\sigma ^{}=[1\alpha (\rho /\rho _0)]\sigma ^{free}`$. The most impressive thing shown in the figure is that the density dependence of $`\sigma _{np}^{}`$ and $`\sigma _{pp(nn)}^{}`$ is very different: At low densities, the $`\sigma _{np}^{}`$ is about $`34`$ times larger than $`\sigma _{pp(nn)}^{}`$ and then decreases quickly with density and finally the $`\sigma _{np}^{}/\sigma _{pp(nn)}^{}`$ ratio approaches to $`1`$ as density reaching about the normal density. The $`\sigma _{np}^{}/\sigma _{pp(nn)}^{}`$ ratio is shown in the small top-right plot. It means that at dense nuclear matter the isospin effect on in-medium ECS almost washes out. For $`Y_p=0.3`$ case, $`\sigma _{pp}^{}`$ and $`\sigma _{nn}^{}`$ differ. The value of $`\sigma _{pp}^{}`$ is larger than that of $`\sigma _{nn}^{}`$ because $`M_p^{}>M^{}>M_n^{}`$ (see Fig. 1 (a)). The splitting of $`\sigma _{pp}^{}`$ and $`\sigma _{nn}^{}`$ increases with density and the degree of isospin asymmetry of the medium following the splitting of the proton and neutron effective mass.

In Ref. Li00  the in-medium $`\sigma _{np}^{}`$ and $`\sigma _{pp(nn)}^{}`$ are calculated based on QHD-II model with the medium correction of the effective mass of $`\rho `$ meson taken into account. We found in Li00  that $`\sigma _{pp(nn)}^{}`$ depends on density obviously but $`\sigma _{np}^{}`$ depends on density weakly. The different results can be attributed to the inclusion of the $`\delta `$ mesons. In the mean field level, the contributions of the $`\delta `$ and $`\rho `$ mesons compensate each other in central potential, producing an effective isovector potential that is comparable in strength to the one obtained in calculations that include only the $`\rho `$ meson. While for the calculations of the in-medium cross sections, which is beyond the mean field, the direct compensation of the contributions from $`\delta `$ and $`\rho `$ mesons is lost. So we may not expect that the comparable results for the in-medium cross sections should be obtained for the cases of inclusion and not inclusion of the $`\delta `$ meson. The different results for the density dependence of $`\sigma _{pp(nn)}^{}`$ and $`\sigma _{np}^{}`$ obtained in two different models may provide us with a more strict check for the model.

Now let us investigate the temperature effect on $`\sigma _{pp(nn)}^{}`$ and $`\sigma _{np}^{}`$. Fig. 4 shows $`\sigma _{nn}^{}`$ (in (a)), $`\sigma _{pp}^{}`$ (in (b)), and $`\sigma _{np}^{}`$ (in (c)) at $`T=0`$, $`50`$, $`100`$, and $`150`$ MeV. In Fig. 4 one sees that the in-medium ECS increases with temperature and the influence of temperature on in-medium ECS increases with density. We should point out that this temperature effect is model dependent. With the model similar to that used in Ref. Li00  we find the temperature effect on the medium correction of nuclear mass and then on the ECS is weaker than the present results. One of the most important information contained in Fig. 4 is that the influence of temperature on in-medium cross sections is much weaker than that of density. The cause of this difference can be attributed to the different behavior of the density and temperature dependence of the nucleon effective mass. This tells us that one should consider the effect of medium density on the in-medium cross section more seriously than that of medium temperature.

In summary, in this work we have studied the density and temperature effects on in-medium elastic cross sections based on the density dependent relativistic hadron field theory in which the $`\delta `$ meson is included. Our results show that at low densities the $`\sigma _{np}^{}`$ is about $`34`$ times larger than $`\sigma _{pp(nn)}^{}`$ and then decreases quickly with density and finally the $`\sigma _{np}^{}/\sigma _{pp(nn)}^{}`$ ratio approaches to $`1`$ as densities reaching and exceeding the normal density. It means that the isospin effect on in-medium elastic cross section almost washes out at high densities. The analysis about the individual contributions from each term contributing to $`\sigma _{pp(nn)}^{}`$ and $`\sigma _{pn}^{}`$ indicates that the isospin effect on the density dependence of the in-medium elastic cross sections is dominantly caused by the delicate balance of the isovector $`\delta `$ and $`\rho `$ mesons. For isospin asymmetric medium, the $`\sigma _{nn}^{}`$ and $`\sigma _{pp}^{}`$ differ following the splitting of the proton and neutron mass due to the inclusion of the $`\delta `$ meson. Our calculation results show that the in-medium ECS increases with the temperature and the temperature effect of nuclear medium on $`\sigma _{np}^{}`$ and $`\sigma _{pp(nn)}^{}`$ is weaker compared with the density effect. The temperature effect enhances as density increases so the influence of temperature to in-medium ECS may have to be taken into account at higher nuclear density.

The work is supported by the National Natural Science Foundation of China under Nos. 10175093 and 10235030, Science Foundation of Chinese Nuclear Industry and Major State Basic Research Development Program under Contract No. G20000774, the Knowledge Innovation Project of the Chinese Academy of Sciences under Grant No. KJCX2-SW-N02, and the CASK.C. Wong Post-doctors Research Award Fund.


format: text

text:
# Photoproduction of 𝜼-mesic 3He

## Abstract

The photoproduction of $`𝜼`$-mesic 3He has been investigated using the TAPS calorimeter at the Mainz Microtron accelerator facility MAMI. The total inclusive cross section for the reaction $`𝜸`$3He$`\mathbf{}`$$`𝜼`$X has been measured for photon energies from threshold to 820 MeV. The total and angular differential coherent $`𝜼`$ cross sections have been extracted up to energies of 745 MeV. A resonance-like structure just above the $`𝜼`$ production threshold with an isotropic angular distribution suggests the existence of a resonant quasi-bound state. This is supported by studies of a competing decay channel of such a quasi-bound $`𝜼`$-mesic nucleus into $`𝝅^\mathrm{𝟎}`$pX. A binding energy of (-4.4$`\mathbf{\pm }`$4.2) MeV and a width of (25.6$`\mathbf{\pm }`$6.1) MeV is deduced for the quasi-bound $`𝜼`$-mesic state in 3He.

It has long been discussed whether the attractive strong interaction among mesons and nucleons may lead to the existence of mesic nuclei. Deeply bound pionic states are known to exist due to the superposition of the attractive electromagnetic and repulsive strong interaction. In the case of the neutral $`𝜼`$-meson an attractive strong $`𝜼`$N interaction would allow for the formation of $`𝜼`$-mesic nuclei. First predictions for such states were based on investigations of the $`𝜼`$-nucleon scattering length $`𝒂_{𝜼𝑵}`$. Its real part can be interpreted as a measure for the scattering of the initial particles while its imaginary part accounts for losses into other channels. First estimates of the $`𝜼`$-N scattering length have been extracted from coupled channel analyses, performed by Bhalerao and Liu Bhalerao and Liu (1985) in 1985. They found a scattering length $`𝒂_{𝜼𝑵}`$ of (0.27+i0.22) fm. The corresponding phase shifts have positive values which indicates an attractive potential. Based on these calculations, Liu and Haider Liu and Haider (1986) claimed that bound states between the $`𝜼`$ and a nucleus should be possible for nuclei with mass number A$`\mathbf{>}`$10.
In 1991, Ueda Ueda (1992) predicted the existence of a quasibound $`𝜼`$NN state with I=0, J=1. These theoretical studies have raised hopes of experimentally discovering quasi-bound states of $`𝜼`$ mesons in light nuclei. If they exist, one may expect them to be narrow in few-nucleon systems. The rapid slope of the near-threshold amplitude in pd$`\mathbf{}`$$`𝜼`$3He was interpreted by Wilkin Wilkin (1993) as an indirect evidence for the existence of $`𝜼`$-mesic nuclei. A recent experiment Sokol et al. (2001) has studied the possible photoproduction of $`𝜼`$-mesic Boron via the reaction $`𝜸\mathbf{+}^{\mathrm{𝟏𝟐}}𝑪\mathbf{}𝒑\mathbf{+}_𝜼^{\mathrm{𝟏𝟏}}𝑩\mathbf{}𝝅^\mathbf{+}\mathbf{+}𝒏\mathbf{+}𝑿`$ with an intermediate $`𝑺_{\mathrm{𝟏𝟏}}`$(1535) resonance which decays into $`𝝅^\mathbf{+}𝒏`$. The correlated energies and momenta of pion and neutron have been used as a signature for the reaction channel.
Concerning the lightest nuclei, photoproduction experiments on the deuteron Krusche et al. (1995); Weiss et al. (2001) and on 4He Hejny et al. (2002) have shown modifications of the near-threshold behaviour via a comparison with models based on an impulse approximation. Sibirtsev et al. Sibirtsev et al. (2002) took final state interactions of the N-N system as well as in the $`𝜼`$-N system into account in their calculation of the quasifree $`𝜼`$ photoproduction off 2H. The data were best reproduced for a scattering length $`𝒂_{𝜼𝑵}`$=(0.42 + i0.34) fm. Other analyses (e.g. Wilkin (1993); Kaiser et al. (1997); Nieves et al. (2001); Green and Wycech (1997)) of the scattering length yield varying results, especially for the real part of the scattering length. The $`𝜼`$-3He scattering length, extracted from photoproduction or hadroproduction, provides information on the possibility of a bound state of $`𝜼`$ and nucleus. In this Letter, the results for coherent $`𝜼`$ production off 3He and the $`𝜼`$-mesic decay channel into $`𝝅^\mathrm{𝟎}`$pX are presented in terms of a possible photoproduction of $`𝜼`$-mesic 3He.

The experiment was performed at the electron accelerator facility Mainz Microtron (MAMI) Ahrens et al. (1994). Photons were produced in a thin radiator foil and their energies were determined using the Glasgow tagged photon spectrometer Hall et al. (1996). The tagged photon energy range was 275-820 MeV with an energy resolution of $`\mathbf{}`$ 2 MeV. The reaction products were detected in the TAPS photon spectrometer Novotny et al. (1991); Gabler et al. (1994). TAPS consisted of six blocks each with 62 hexagonally shaped BaF2 scintillation detectors arranged in an 8x8 matrix and a forward wall with 138 BaF2 detectors in an upright 11x14 arrangement. The length of the crystals is 250 mm ($`\mathbf{}`$12 radiation lengths) with an inner diameter of 59 mm plus a cylindrical end cap of 25 mm length (54 mm diameter). The crystals were read out by photomultipliers. The six blocks were located in a horizontal plane around the target point at angles of 153, 104, 55, -54, -103, -152 with respect to the beam axis at distances of 55-59 cm. The forward wall was positioned at 0 at a distance of 55 cm. The solid angle covered by this setup is around 40% of the full solid angle Kotulla (2003). The target had an effective length of 115 mm and a diameter of 43 mm. It was filled with liquid 3He. More detail can be found in Mahmoud (1999).

For the identification of the $`𝜼`$ and $`𝝅^\mathrm{𝟎}`$ mesons the double photon decay channel was used. The two photon invariant mass was calculated from the momenta of the measured photons. The experimental resolution is 21 MeV (FWHM) for $`𝝅^\mathrm{𝟎}`$ detection and 60 MeV for the $`𝜼`$ meson. The identification of protons is achieved via the characteristic relation between deposited energy and time of flight, as well as pulse shape analysis Novotny et al. (1996). The calibration of the proton energy was performed by calculating the missing mass for the reaction $`𝜸`$3He$`\mathbf{}`$$`𝜼`$pd. Energy losses of the protons and quenching effects in the detector have been corrected for such that the missing mass corresponds to that of the proton. Random coincidences in the tagging spectrometer have been taken into account by subtracting background events before and after the prompt time coincidence peak.
The inclusive $`𝜼`$ cross section was deduced from the rate of $`𝜼`$ events divided by the number of target nuclei per cm2, the photon flux, the branching ratio of $`𝜼`$ into 2$`𝜸`$, and the detection and analysis efficiency. These efficiencies were determined by Monte Carlo simulations using the GEANT3 code. The photon flux was measured by counting the number of scattered electrons in the tagger. The missing energy, Emiss, for the coherent reaction $`𝜸`$3He$`\mathbf{}`$$`𝜼`$3He was obtained by subtracting the calculated $`𝜼`$ energy from two body kinematics from the measured $`𝜼`$ energy in the $`𝜼`$3He center of momentum system. Coherent events should be centered at $`𝑬_{𝒎𝒊𝒔𝒔}`$= 0 MeV while break-up events are distributed at negative values of the missing energy.

The fraction of coherent $`𝜼`$ production (thin solid line in fig. 1) and quasifree $`𝜼`$ production (dashed line) was obtained by fitting simulations for both production mechanisms to the experimental data. The simulation for the quasifree production includes the elementary production cross section on the proton folded with the Fermi momentum distribution. The resulting coherent $`𝜼`$ cross section is plotted in fig. 2. The peak-like behaviour of the near threshold cross section can be interpreted as a first sign for a modified $`𝜼`$ production process. To investigate this behaviour further, the angular distribution for the two energy ranges 600-625 MeV and 625-745 MeV are plotted in the lower part of fig. 2. While the latter distribution agrees with the expected angular dependence (d$`𝝈`$/d$`𝛀`$$`\mathbf{}`$F2(q2)), the distribution at threshold is unexpectedly isotropic. Here, F2(q2) is the 3He form factor depending on the momentum transfer q. The observations are consistent with the assumption of an $`𝜼`$-mesic nucleus which isotropically decays into the coherent $`𝜼`$ channel.
There are different theoretical efforts for describing the experimental data. A plane wave impulse approximation (PWIA) calculation by Kamalov and Tiator Kamalov et al. (1993); Tiator et al. (1994) does not deal with final state interaction (FSI) effects (dotted line in fig. 2). FSI are explicitly taken into account in the calculations of Shevchenko and collaborators Shevchenko et al. (2003, 2002) The model employs the finite rank approximation (FRA) which excludes virtual excitations of the nucleus during the interaction with the $`𝜼`$ meson. The calculation is performed in such a way as to reproduce the real part of the scattering length to be 0.75 fm (dashed and dot-dashed lines in fig. 2). A third calculation has been performed by Fix and Arenhövel Fix and Arenhövel (2003). In contrast to the assumptions of Shevchenko which neglect target excitations, a two-body quasi particle separation is applied to the $`𝜼`$-3N problem. The analysis obtains values for the $`𝜼`$-3He scattering length that do not allow a bound state of $`𝜼`$ and nucleus. However, the existence of a virtual pole near threshold is claimed - implying a quasi-resonance with positive binding energy that modifies the threshold behaviour of the coherent cross section. The cross section shows a sharp cusp effect at threshold but stays a factor of 2 below the measured data (solid line in fig. 2).
Earlier measurements of the $`𝜼`$ production off 3He Mayer et al. (1996); Bilger et al. (2002) in a different initial state - the proton induced reaction on the deuteron - showed a similar energy dependence when dividing the different phase space factors out and thus comparing the average amplitude as a function of the $`𝜼`$ momentum in the CM frame (fig. 3). Only the two data points between 0.6 and 0.8 fm-1 seem to be significantly lower in the photoproduction data compared to the proton induced case. Also shown in fig. 3 is an optical model fit to the data below 0.5 fm-1.
Additional information about a possible $`𝜼`$-mesic state can be deduced by studying another decay channel for an $`𝜼`$-mesic nucleus. The quasi-bound system will decay into $`𝜼`$3He for energies above the $`𝜼`$ production threshold. Since the phase space is limited, the momentum of the decaying $`𝜼`$ meson will preferentially be low and hence there is a high probability for the $`𝜼`$ to be captured by a nucleon into an S11(1535) resonance. This resonance has some 50% branching ratio into $`𝜼`$N but also roughly the same probability to decay into $`𝝅`$N. The advantage of the latter decay mode is that the $`𝜼`$ meson is acting as an intermediate particle only. Only the lighter $`𝝅`$ meson is emitted in the final state. Hence the reaction can also take place below the $`𝜼`$ production threshold. Below threshold, the decay into $`𝜼`$3He is not possible and the decay proceeds exclusively into the $`𝝅`$N channel.

The signature for the decay of an $`𝜼`$-mesic nucleus into $`𝝅`$N is given by correlated pairs of $`𝝅^\mathrm{𝟎}`$ and p. With a low momentum $`𝜼`$ in the intermediate state, the decaying S11 will also be at low momentum. Thus, the relative angle of $`𝝅^\mathrm{𝟎}`$ and proton has to be near 180 in the $`𝜸`$-3He CM system. The triangles in the left panel of fig. 4 show the excitation function for relative angles between 170 and 180 as a function of $`𝑾`$=$`\sqrt{\mathrm{𝟐}\mathbf{}𝑬_𝜸\mathbf{}𝒎_{{}_{}{}^{\mathrm{𝟑}}𝑯𝒆}\mathbf{+}𝒎_{{}_{}{}^{\mathrm{𝟑}}𝑯𝒆}^\mathrm{𝟐}}\mathbf{}𝒎_𝒅\mathbf{}𝑩_{{}_{}{}^{\mathrm{𝟑}}𝑯𝒆}`$. Here, the total CM energy is reduced by the mass $`𝒎_𝒅`$ of the residual nucleus and the 3He binding energy of $`𝑩_{{}_{}{}^{\mathrm{𝟑}}𝑯𝒆}`$= -5.5 MeV (the $`𝜼`$-$`{}_{}{}^{\mathrm{𝟑}}𝑯𝒆`$ binding energy for a structure at resonance mass $`𝑾_𝑹`$ would then be given as $`𝑩_𝜼\mathbf{=}𝑾_𝑹\mathbf{}\mathbf{(}𝒎_𝜼\mathbf{+}𝒎_𝒑\mathbf{)}`$). The reaction channel is dominated by the quasifree $`𝝅^\mathrm{𝟎}`$ production. This background channel is determined by subtracting the scaled $`𝝅^\mathrm{𝟎}`$p excitation function for relative opening angles of 150 to 170. A remaining structure just below the production threshold for $`𝜼`$ mesons may indicate the formation for an $`𝜼`$-mesic 3He nucleus. Its peak position is $`𝑾_𝑹`$=1483 MeV ($`\mathbf{\pm }`$ 5 MeV) with a width $`𝚪`$ = 39 MeV ($`\mathbf{\pm }`$ 21 MeV). The experimental resolution is given by the tagger spectrometer resolution of about 2 MeV. The statistical significance of the structure $`𝝈\mathbf{=}\frac{𝑺}{\sqrt{\mathbf{(}𝚫𝑺\mathbf{)}^\mathrm{𝟐}\mathbf{+}\mathbf{(}𝚫𝑩𝑮\mathbf{)}^\mathrm{𝟐}}}`$ (S is the number of signal counts and BG the corresponding background under the peak) is $`\mathbf{}`$3.5$`𝝈`$. fig. 5 shows a comparison of both decay channels with a simultaneous fit of both channels. The result, $`𝑾_𝑹`$=1481.2 MeV ($`\mathbf{\pm }`$ 4.2 MeV) and $`𝚪`$ = 25.3 MeV ($`\mathbf{\pm }`$ 6.1 MeV), is consistent with the results given in fig. 4 but reduces the error bars on the extracted parameters. The different phase space factors for the two decay channels were taken into account by assuming a single intermediate S11(1535) resonance neglecting multiscattering processes. The fit function is a Breit-Wigner distribution $`𝒇_{𝑩𝑾}\mathbf{(}𝑾\mathbf{)}\mathbf{=}\frac{\frac{\mathrm{𝟏}}{\mathrm{𝟒}}𝚪^\mathrm{𝟐}}{\mathbf{(}𝑾\mathbf{}𝑾_𝑹\mathbf{)}^\mathrm{𝟐}\mathbf{+}\frac{\mathrm{𝟏}}{\mathrm{𝟒}}𝚪^\mathrm{𝟐}}`$ folded with the corresponding phase space factors

$`𝑩_{𝝅\mathbf{,}𝜼}\mathbf{(}𝑾\mathbf{)}\mathbf{=}𝒃_{𝝅\mathbf{,}𝜼}\mathbf{}\sqrt{\frac{\mathbf{(}𝑾^\mathrm{𝟐}\mathbf{}\mathbf{(}𝒎_𝒑\mathbf{}𝒎_{𝝅\mathbf{,}𝜼}\mathbf{)}^\mathrm{𝟐}\mathbf{)}\mathbf{}\mathbf{(}𝑾^\mathrm{𝟐}\mathbf{}\mathbf{(}𝒎_𝒑\mathbf{+}𝒎_{𝝅\mathbf{,}𝜼}\mathbf{)}^\mathrm{𝟐}\mathbf{)}}{\mathrm{𝟒}𝑾^\mathrm{𝟐}}}\mathbf{.}`$

Here, $`𝒃_{𝝅\mathbf{,}𝜼}`$ are constants that fix the partial widths for the decay of the $`𝑺_{\mathrm{𝟏𝟏}}`$ into $`𝝅`$ and $`𝜼`$ to 75 MeV at the pole mass of 1535 MeV, respectively. Further, $`𝒎_𝒑`$ is the proton mass and $`𝒎_{𝝅\mathbf{,}𝜼}`$ are the meson masses.
In summary, we have investigated the photoproduction of $`𝜼`$-mesic 3He by studying two possible decay channels - the coherent $`𝜼`$ production and the decay into $`𝝅^\mathrm{𝟎}`$p. Both decay channels show indications for photoproduction of $`𝜼`$-mesic 3He. In the $`𝜼`$ channel, a resonance-like behaviour of the coherent cross section associated with an isotropic angular distribution in the threshold regime is observed. In contrast, a strong forward peaking is expected for coherent $`𝜼`$ production from form factor considerations. In the $`𝝅^\mathrm{𝟎}`$p channel, correlated $`𝝅^\mathrm{𝟎}`$-proton pairs with relative angles near 180 in the CM system have been observed that give rise to a peak-like structure at energies slightly below the $`𝜼`$ production threshold. The extracted resonance parameters (binding energy=(-4.4 $`\mathbf{\pm }`$ 4.2) MeV, full width=(25.6 $`\mathbf{\pm }`$ 6.1) MeV) are consistent with expectations for $`𝜼`$-mesic nuclei.
It is a pleasure to acknowledge inspiring discussions with W. Cassing and C. Wilkin. We thank the accelerator group of the Mainz Microtron MAMI, as well as the other technicians and scientists of the Institut für Kernphysik at the Universität Mainz for the outstanding support. This work was supported by DFG SPP 2034, SFB 221, SFB 443, the U.K. Engineering and Physical Sciences Research Council, and Schweizerischer Nationalfonds.


text:
# Diffusion in stellar interiors: critical tests of three numerical methods

## 1 Introduction

For many physical processes in stellar interiors as heat transfer, transport of angular momentum, mixing of chemical elements through time dependent convection, semiconvection or turbulence, a diffusion equation needs to be solved. Diffusion is thus involved in numerous interesting astrophysical problems. For instance, Vauclair and Charbonnel (Vau98 (1998)) discuss the importance of selective diffusion in low–metallicity stars and the consequences for lithium primordial abundance; Charbonnel (Cha95 (1995)) has shown how rotational diffusion in low mass stars may lead to the destruction of 3He and thus to new insights on the evolution of the abundance of this cosmological element; most of the results on the effects of rotation in stellar interiors are based on the resolution of a diffusion equation for both the transport of the chemical species and of the angular momentum (Chaboyer & Zahn Cha92 (1992); Zahn Za92 (1992); Talon & Zahn Ta97 (1997); Denissenkov et al. Deni99 (1999); Heger et al. He20 (2000);Maeder & Meynet MM01 (2001); Meynet & Maeder MM00 (2003) ). Diffusion plays also a key role in our understanding of the structure and the cooling of white dwarfs (cf. Kawaler Ka95 (1995)).

For solving numerically the diffusion equation, different methods are available. Quite generally, they can be classed into two main categories: the finite differences methods and the finite elements methods. Depending on how the time discretisation is performed, one distinguishes three subclasses in each of these two main categories, namely the subclasses of the explicit, implicit and of the “Crank–Nicholson” type methods. Thus, at least six different methods are available. Here we are interested in modelling the long term evolution of stars (secular evolution). This immediately prevents us of using explicit methods which requires the use of much too short time steps (see below). Among the four remaining methods, the “Crank–Nicholson” finite elements method was not tested in view of the results obtained by the “Crank–Nicholson” finite differences method (see Sect. 5.4). Thus we shall focus our attention on three of them namely

* the implicit finite elements method,
* the “Crank–Nicholson” finite differences method,
* the implicit finite differences method.

Descriptions of these methods as well as a general discussion of their respective advantages and weaknesses can be found for instance in Press et al. (Teu92 (1992)) for the finite differences methods and in Zienkiewicz & Taylor (zi00 (2000)) for the finite elements methods. If these considerations can help in choosing the best method in general, we think that very valuable complementary information can be gained by applying the methods to the resolution of realistic astrophysical cases, where the different timescales involved are very specific. This is what we propose here and this is the main aim of this work. Let us recall that, as said in Press et al (Teu92 (1992)) “… differencing partial differential equations is an art as much as a science”, and that a numerical scheme, which from a theoretical point of view does appear as very robust can show some weaknesses when applied to a real particular situation. This is precisely why it is necessary to test realistic astrophysical situations. This allows us to make a new and appropriate analysis to support the choice of the best method to handle diffusion problems.

The interested reader will also find in this paper a detailed description of the discretized equations for resolving the diffusion of chemical elements and of the angular momentum in stars. We also propose a new way to estimate the diffusion coefficient in the interface between a convective and a radiative zone.

In Sect. 2 we briefly recall the physical problem to be resolved. The three numerical methods are described in details in Sect. 3. The estimate of the diffusion coefficient at the interface between a convective and a radiative zone is presented in Sect. 4. Sect. 5 is devoted to various tests and comparisons. Section 6 summarizes the main results.

## 2 The physical problem

Diffusion is a process by which components of a mixture move from one part of a system to another as a result of random motion. Let us briefly recall how the diffusion velocity is defined (see e.g. Battaner ba96 (1996)). In a multicomponent fluid, one defines the mean velocity $`v_0`$ of the mixture as

$$v_0=\frac{\underset{i}{}m_in_i<v_i>}{_im_in_i},$$

where

$$<v_i>=\frac{1}{n_i}_pf_iv_id\tau _p,$$

is the mean velocity of the $`i`$th component, $`n_i`$ is the partial number density of particles $`i`$, $`m_i`$ their mass, $`v_i`$ their velocity and $`f_i(x,p,t)`$ the probability distribution function, which gives the probability to find a particle $`i`$ at the position $`x`$, with the momentum $`p`$, at the time $`t`$; $`\tau _p`$ is the phase space volume for the particle’s momentum. One can wonder, why in the expression for $`v_0`$, the weighting factor is the partial mass density and not the partial number density. The question could be formulated in another way: when is the mean velocity of the mixture zero ? When the net flux of the number of particles is zero ($`_in_i<v_i>=0`$), or when the net flux of the mass is zero ($`_im_in_i<v_i>=0`$) ? Obviously when the net flux of the mass is zero, otherwise the centre of gravity of the system is moving and thus the mean velocity of the mixture is not zero.

The peculiar velocity of a particle is defined as

$$V_i=v_iv_0.$$

The mean value of $`V_i`$ is

$$<V_i>=\frac{1}{n_i}_pf_iV_id\tau _p=<v_i>v_0.$$
(1)
This quantity is called the diffusion velocity of the $`i`$th component. It can be shown very easily, using the expression for $`v_0`$ that the different diffusion velocities must compensate one another, i.e.

$$\underset{i}{}m_in_i<V_i>=0.$$
(2)
There is thus no net mass flux associated to diffusion. In that respect, in a star, the changes of chemical composition due to diffusion can be treated very similarly as those due to the nuclear reactions. Like the nuclear reactions, diffusion modifies the mean molecular weight (and thus the hydrostatic structure of the star), but does not induce any net mass flux.

In the context of the Boltzmann microscopic theory, it is possible to deduce from the equations of motion for the different components, expressions for the $`<V_i>`$’s (see Battaner ba96 (1996)). We shall not repeat these developments here. Instead we consider that diffusion is equivalent to a displacement of particles whose velocities may be related to the diffusive coefficient, $`D`$, by the expression (see also Sect. 3.2)

$$<V_i>=\frac{D}{X_i}X_i,$$

where $`X_i`$ is the mass fraction of element $`i`$. The minus sign results from the fact that the velocity is directed in the opposite direction to the mass fraction gradient. This expression automatically satisfies the condition seen above that $`_im_in_i<V_i>=0`$, indeed

$$\underset{i}{}m_in_i<V_i>=\rho D\left(\underset{i}{}X_i\right)=0,$$

since $`_iX_i=1`$ and the relation between $`n_i`$ and $`X_i`$ is

$$n_i=\frac{\rho X_i}{A_im_H},$$

where $`\rho `$ is the density, $`A_i`$ the atomic weight expressed in units of proton mass $`m_H`$.

We are interested in resolving the diffusion equation in spherically symmetric systems. Let us consider a one-dimensional problem, where particles $`i`$ may diffuse along the radial direction $`r`$. Thus

$$<V_i>=\frac{D}{X_i}\frac{X_i}{r}.$$
(3)
The variation of the number of particles in the element of volume $`\upsilon `$ and of surface $`S`$, due only to diffusion, is

$$\frac{}{t}_\upsilon n_id\upsilon =_Sn_i<V_i>dS.$$
(4)
From the divergence theorem and the expression of $`<V_i>`$ just seen above one obtains

$$\frac{n_i}{t}=\frac{1}{r^2}\frac{}{r}(r^2n_i<V_i>)=\frac{1}{r^2}\frac{}{r}\left(\rho r^2D\frac{\left(\frac{n_i}{\rho }\right)}{r}\right).$$
(5)
Replacing $`n_i`$ by its expression as a function of $`X_i`$ in Eq. (5) gives

$$\frac{}{t}(\rho X_i)=\frac{1}{r^2}\frac{}{r}(\rho r^2D\frac{X_i}{r}).$$
(6)
Thus Eq. (6) results from the expressions for the $`<V_i>`$ (Eq. 3) and from the continuity equation for the $`n_i`$’s (Eq. 4).

Sometimes, Eq. (6) is written as below

$$\frac{}{t}(\rho c_i)=\frac{1}{r^2}\frac{}{r}(\rho r^2D\frac{c_i}{r}),$$
(7)
with partial concentration $`c_i`$ instead of mass fraction, without indication on the precise meaning of what concentration means. Partial concentration in mass is equivalent to the mass fraction, while partial concentration in number is equal to $`n_i/n`$ where $`n=_in_i`$. Replacing $`c_i`$ in Eq. (7), by $`n_i/n`$ and expressing $`n_i`$ as a function of $`X_i`$, one obtains

$$\frac{}{t}(\rho \mu X_i)=\frac{1}{r^2}\frac{}{r}\left(\rho r^2D\frac{}{r}(\mu X_i)\right),$$
(8)
where $`\mu =\rho /(nm_H)`$ is the mean molecular weight of the ions. This equation is identical to Eq. (6) only when $`\mu `$ remains constant as a function time and is constant as a function of $`r`$. Thus Eq. (8) is equivalent to Eq. (6) only when applied to minor constituents, the abundances of which do not affect $`\mu `$ and when $`\mu `$ has no gradient. In more general cases, Eq. (8) is not equivalent to Eqs. (5) and (6).

This can also be seen in a slightly different way. When, in Eq. (7), the $`c_i`$’s are identified with the number fractions, one obtains

$$\frac{}{t}\left(\rho \frac{n_i}{n}\right)=\frac{1}{r^2}\frac{}{r}\left(r^2\rho \frac{n_i}{n}w_i\right),$$
(9)
where we have introduced a velocity $`w_i=\frac{D}{n_i/n}\frac{}{r}\left(\frac{n_i}{n}\right)`$. This velocity is not a diffusive velocity in the sense that $`_im_in_iw_i=nDm_H\mu /r`$ is not equal to zero except in zones where $`\mu `$ is constant. Eq. (9) expresses the change of $`\rho n_i/n=\mu n_im_H`$ in an element of volume resulting from the transport of the quantity $`\rho n_i/n`$. But physically, this is not $`\rho n_i/n`$ which diffuses but the particles themselves. As above, there are two conditions for Eq. (9) to be equivalent to Eqs. (5) and (6): 1) the particles which diffuse must have a very small abundance, so small that their diffusion does not affect $`\mu =\rho /(nm_H)`$; 2) and $`\mu `$ does not vary with the radius. We conclude thus that when the diffusive transport is described by an equation of the form given in Eq. (7), the concentrations are equivalent to mass fractions. The identification with number fractions results in unphysical description of the process except in very particular situations.

Starting from Eq. (6) and summing over all the chemical species $`i`$, one obtains

$$\frac{\rho }{t}=\frac{1}{r^2}\frac{}{r}\left(\rho r^2D\frac{}{r}(1)\right)=0.$$

This is quite consistent with the fact that the diffusive velocities must compensate each other. Thus the density can be put outside from the time derivative in Eq. (6):

$$\rho \frac{X_i}{t}|_{m_r}=\frac{1}{r^2}\frac{}{r}\left(\rho r^2D\frac{X_i}{r}\right),$$
(10)
where $`m_r`$ is the langrangian mass coordinate. This is the equation we shall numerically resolve in this paper. The conditions at the centre and the surface are:

$$\frac{X_i}{r}|_{m_r=0}=\frac{X_i}{r}|_{m_r=M}=0,$$
(11)
where $`M`$ is the total mass of the star.

The equation expressing the diffusion of the angular momentum is (see e.g. Endal & Sofia En78 (1978))

$$\rho \frac{(r^2\mathrm{\Omega })}{t}|_{m_r}=\frac{1}{r^2}\frac{}{r}\left(\rho r^4D^{}\frac{\mathrm{\Omega }}{r}\right),$$
(12)
where $`\mathrm{\Omega }`$ is the angular velocity and $`D^{}`$ the diffusion coefficient for the angular momentum. In our numerical experiments, the radii do not change with time111Let us note that in more realistic stellar models, the radii of course vary as a function of time, however time steps can be chosen sufficiently small for considering them to be constant during one time step. The same can be said for the diffusion coefficient and other structure variables as for instance the density.. Eq. (12) then becomes

$$\rho r^2\frac{\mathrm{\Omega }}{t}|_{m_r}=\frac{1}{r^2}\frac{}{r}\left(\rho r^4D^{}\frac{\mathrm{\Omega }}{r}\right)=\mathrm{div}\left(\rho r^2D^{}(\mathrm{\Omega })\right).$$
(13)
The conditions at the centre and at the surface are

$$\frac{\mathrm{\Omega }}{r}|_{m_r=0}=\frac{\mathrm{\Omega }}{r}|_{m_r=M}=0.$$
(14)
Let us notice that in general the transport equation for the angular momentum may contain other terms expressing the advection of angular momentum by meridional circulation and the effects of magnetic braking. As these last two terms bring no specific problems with respect to diffusion, we do not consider them here. However, the advection of angular momentum needs a particularly careful treatment as well (cf. Meynet & Maeder MMV (2000)).

## 3 The three tested methods

### 3.1 The implicit finite elements method

The detailed description of the implicit finite elements method used to resolve the diffusion equation for the chemical species is presented in Schatzman et al. (Sch81 (1981)) (see the appendix by Glowinsky & Angrand). We present here the procedure for the case of the angular momentum diffusion. The basic idea of this method is to decompose the unknown function, here $`\mathrm{\Omega }(m_r,t)`$, as a linear combination of well chosen independent functions. This decomposition can be performed in many different ways. We adopt here the same decomposition as in Schatzman et al. (Sch81 (1981)). Of course, the conclusions concerning the ability of this method to provide physical solutions will only refer to this particular choice.

We decompose the star in $`K`$ shells, with the langrangian mass coordinate of the ith mesh point being $`m_i`$ ($`m_i`$ is the mass inside the sphere of radius $`r_i`$). In the following all the quantities with an indice $`i`$ are evaluated at the mesh point $`i`$. The shells are numbered from 1 at the surface to $`K`$ at the centre. Let us introduce $`K`$ functions $`b_i(m_r)`$ defined by

$$\begin{array}{cccc}& & \frac{m_rm_{i1}}{m_im_{i1}}\hfill & \mathrm{if}m_r[m_i,m_{i1}],\hfill \\ & \hfill & & \\ b_i(m_r)=& \hfill & \frac{m_rm_{i+1}}{m_im_{i+1}}\hfill & \mathrm{if}m_r[m_{i+1},m_i],\hfill \\ & \hfill & & \\ & & 0\hfill & \mathrm{if}m_r[m_{i+1},m_{i1}].\hfill \end{array}$$

The function $`b_i`$ is equal to one at $`m_r=m_i`$, is equal to zero at $`m_{i+1}`$ and $`m_{i1}`$ and varies linearly as a function of $`m_r`$ inbetween. By multiplying Eq. (13) by each of the functions $`b_i(m_r)`$, one obtains $`K`$ equations. The integration over the whole volume of the star, $`V`$, of each of these $`K`$ equations gives

$$_V\rho r^2\frac{\mathrm{\Omega }}{t}b_idV=_V\mathrm{div}\left(\rho r^2D^{}\mathrm{𝐠𝐫𝐚𝐝}(\mathrm{\Omega })\right)b_idV.$$
(15)
Using the general relations:

$$\mathrm{div}(a𝐯)=a\mathrm{div}(𝐯)+\mathrm{𝐠𝐫𝐚𝐝}(a)𝐯,$$
(16)
$$_V\mathrm{div}(𝐯)dV=_S𝐯d𝐒,$$
(17)
where $`a`$ is a scalar, $`𝐯`$ a vector and $`S`$ the surface of the volume $`V`$, one obtains that the right hand term of Eq. (15) becomes

$$_S\rho r^2D^{}\mathrm{𝐠𝐫𝐚𝐝}(\mathrm{\Omega })b_id𝐒_V\rho r^2D^{}\mathrm{𝐠𝐫𝐚𝐝}(\mathrm{\Omega })\mathrm{𝐠𝐫𝐚𝐝}(b_i)dV.$$

The integral on the surface is null, since $`\mathrm{𝐠𝐫𝐚𝐝}(\mathrm{\Omega })=\mathrm{𝟎}`$ on the surface. Thus one has $`K`$ equations of the type

$$_V\rho r^2\frac{\mathrm{\Omega }}{t}b_idV+_V\rho r^2D^{}\frac{\mathrm{\Omega }}{r}\frac{b_i}{r}dV=0.$$
(18)
Using $`\mathrm{d}m_r=\rho \mathrm{d}V=4\pi r^2\rho \mathrm{d}r`$, one obtains

$$_0^M\left[r^2\frac{\mathrm{\Omega }}{t}b_i+r^2D^{}(4\pi r^2\rho )^2\frac{\mathrm{\Omega }}{m_r}\frac{b_i}{m_r}\right]dm_r=0.$$
(19)
The functions $`b_i`$ constitute a set of independent functions, therefore the unknown function $`\mathrm{\Omega }(m_r,t)`$ can be expressed as a linear combination of $`b_i`$’s. One can write

$$\mathrm{\Omega }(m_r,t)=\underset{j=1}{\overset{K}{}}\mathrm{\Omega }_jb_j(m_r),$$
(20)
where $`\mathrm{\Omega }_j=\mathrm{\Omega }(m_j,t)`$. Equation (20) simply says that to obtain $`\mathrm{\Omega }(m_r,t)`$, with $`m_r[m_{i+1},m_i]`$, one simply interpolates linearly as a function of mass between $`\mathrm{\Omega }_{i+1}`$ and $`\mathrm{\Omega }_i`$.

Expressing $`\mathrm{\Omega }(m_r,t)`$ as indicated in Eq. (20), Eq. (19) becomes

$$\underset{j}{}A_{ij}\frac{\mathrm{\Omega }_j}{t}+\underset{j}{}B_{ij}\mathrm{\Omega }_j=0.$$
(21)
where

$$A_{ij}=_0^Mr^2b_ib_jdm_r,$$
(22)
$$B_{ij}=_0^Mr^2D^{}(4\pi r^2\rho )^2\frac{b_i}{m_r}\frac{b_j}{m_r}dm_r.$$
(23)
One has used the facts that the $`\mathrm{\Omega }_j`$’s do not depend on $`m_r`$ and the $`b_j`$’s do not depend on time. When $`|ji|>1`$, there is no overlap between the functions $`b_i`$ and $`b_j`$, thus

$$A_{ij}=0,\mathrm{for}|ji|>1.$$
(24)
For a given value of $`i`$ only, $`A_{i,i1}`$, $`A_{i,i}`$ and $`A_{i+1,i}`$ are different from zero. The same is true for the $`B_{ij}`$. Let us estimate $`A_{i,i1}`$:

$$A_{i,i1}=_{m_i}^{m_{i1}}r^2\frac{m_rm_i}{m_{i1}m_i}\frac{m_rm_{i1}}{m_im_{i1}}dm_r,$$
(25)
one approximates $`r^2`$ by $`0.5\left(r^2(m_{i1})+r^2(m_i)\right)`$, then, by trivial integration, one obtains

$$A_{i,i1}=0.5\left(r^2(m_{i1})+r^2(m_i)\right)\frac{m_im_{i1}}{6}.$$
(26)
The other matrix elements are obtained in the same way:

$$A_{i+1,i}=0.5\left(r^2(m_{i+1})+r^2(m_i)\right)\frac{m_{i+1}m_i}{6},$$
(27)
$$A_{i,i}=2(A_{i,i1}+A_{i+1,i}).$$
(28)
At the centre and at the surface, one has

$$\begin{array}{cc}A_{K,K}\hfill & =0.5r^2(m_{K1})\frac{m_{K1}}{3},\hfill \\ A_{K,K1}\hfill & =0.5A_{K,K},\hfill \\ A_{2,1}\hfill & =0.5A_{1,1},\hfill \\ A_{1,1}\hfill & =0.5\left(r^2(1)+r^2(2)\right)\frac{m_1m_2}{3}.\hfill \end{array}$$

Let us set $`\overline{D}=r^2(4\pi r^2\rho )^2D^{}`$ and $`\overline{D}_{i,i1}`$ an appropriate mean value of $`\overline{D}`$ between the shells $`i`$ and $`i1`$ (see Sect. 4), then the $`B_{i,j}`$ becomes

$$B_{i,i1}=\frac{\overline{D}_{i,i1}}{m_im_{i1}},$$
(29)
$$B_{i+1,i}=\frac{\overline{D}_{i+1,i}}{m_{i+1}m_i},$$
(30)
$$B_{i,i}=B_{i,i1}B_{i+1,i}.$$
(31)
At the centre and at the surface, one obtains

$$\begin{array}{cc}B_{K,K}\hfill & =\frac{\overline{D}_{K,K1}}{m_{K1}},\hfill \\ B_{K,K1}\hfill & =B_{K,K},\hfill \\ B_{2,1}\hfill & =B_{1,1},\hfill \\ B_{1,1}\hfill & =\frac{\overline{D}_{2,1}}{m_1m_2}.\hfill \end{array}$$

Adopting an implicit discretisation in time, we have for Eq. (21)

$$\underset{j}{}A_{i,j}\left(\frac{\mathrm{\Omega }_j^{n+1}\mathrm{\Omega }_j^n}{\mathrm{\Delta }t}\right)+\underset{j}{}B_{i,j}\mathrm{\Omega }_j^{n+1}=0,$$
(32)
or

$$\underset{j}{}\left(\frac{A_{i,j}}{\mathrm{\Delta }t}+B_{i,j}\right)\mathrm{\Omega }_j^{n+1}=\underset{j}{}\frac{A_{i,j}}{\mathrm{\Delta }t}\mathrm{\Omega }_j^n,$$
(33)
where $`\mathrm{\Omega }_j^n`$ is equal to $`\mathrm{\Omega }(m_j,t^n)`$. This linear system of equations is then solved by standard procedure.

### 3.2 The Crank–Nicholson finite differences method

Let us consider in a first step the case of the diffusion of the chemical elements. Multiplying Eq. (10) by $`4\pi r^2`$ and integrating with respect to the spatial coordinate $`r`$ from $`r\mathrm{\Delta }r/2`$ to $`r+\mathrm{\Delta }r/2`$, one obtains

$$\mathrm{\Delta }m\frac{X_i}{t}|_r=\left(4\pi \rho r^2D\frac{X_i}{r}\right)|_{r\mathrm{\Delta }r/2}^{r+\mathrm{\Delta }r/2},$$
(34)
where $`\mathrm{\Delta }m=_{r\mathrm{\Delta }r/2}^{r+\mathrm{\Delta }r/2}4\pi r^2\rho 𝑑r`$. The quantity $`\mathrm{\Delta }r`$ has been chosen sufficiently small for $`\frac{X_i}{t}`$ to remain constant over the spatial range $`\mathrm{\Delta }r`$ around $`r`$. This integration reduces the equation to a first order differential equation, which simply expresses the fact that the time variation of the mass of element $`i`$ in a spherical shell is equal to the difference between the mass of element $`i`$ which enters in the shell and the mass of this same element which goes out from the shell. As indicated in Sect. 2, one can interpret the quantity,

$$\frac{D}{X_i}\frac{X_i}{r}=D\frac{\mathrm{ln}X_i}{r}=<V_i>,$$
(35)
as a diffusion velocity, $`<V_i>`$, for the element $`i`$ (let us recall that the minus sign appears because, when the spatial abundance gradient is positive, the velocity is directed inward). In that case, Eq. (34) writes

$$\mathrm{\Delta }m\frac{X_i}{t}|_r=4\pi \rho r^2X_i<V_i>|_{r\mathrm{\Delta }r/2}^{r+\mathrm{\Delta }r/2},$$

which makes clearly appear the difference of the mass fluxes between the two borders of the shell. This is a more elementary description of diffusion.

One can also derive Eq. (35) from the continuity equation. Indeed supposing that the bulk gas velocity is zero, that $`\rho /t=0`$ (stationary situation), and that there is no sink/source terms of matter, the continuity equation becomes

$$\rho \frac{X_i}{t}|_r=\frac{1}{r^2}\frac{}{r}(\rho r^2<V_i>X_i).$$
(36)
By equating the right hand sides of Eqs. (10) and (36), one obtains for $`<V_i>`$ Eq. (35) above.

Let us now consider the discretized abundance profile of element $`i`$ sketched in Fig. 1. The mass of element $`i`$ removed from point $`k+1`$ and added to point $`k`$ per unit time is (see Eq. 34):

$$4\pi r_{k+1/2}^2\rho _{k+1/2}D_{k+1/2}\frac{X_i}{r}|_{k+1/2},$$
(37)
where the physical quantities with a lower subscript $`k+1/2`$ are estimated at the midpoint between the mesh points $`k+1`$ and $`k`$. For instance, one takes $`r_{k+1/2}=0.5(r_k+r_{k+1})`$. For the diffusion coefficient, we use the procedure exposed in Sect. 4. A similar expression can be written for the mass of element $`i`$ removed from point $`k`$ and added to point $`k1`$. Let us set $`X_{i,k}`$, the mass fraction of element $`i`$ evaluated at the mesh point $`k`$. Thus one has,

$`\mathrm{\Delta }m_k{\displaystyle \frac{(X_{i,k})}{t}}=`$ (38)
$`4\pi r_{k+1/2}^2\rho _{k+1/2}D_{k+1/2}{\displaystyle \frac{X_i}{r}}|_{k+1/2}`$
$`+4\pi r_{k1/2}^2\rho _{k1/2}D_{k1/2}{\displaystyle \frac{X_i}{r}}|_{k1/2},`$
where $`\mathrm{\Delta }m_k=_{r_{k+1/2}}^{r_{k1/2}}4\pi r^2\rho 𝑑r.`$ Let us now discretize this equation as a function of time. The left hand term of Eq. (38) can be written

$`\mathrm{\Delta }m_k{\displaystyle \frac{(X_{i,k})}{t}}\mathrm{\Delta }m_k{\displaystyle \frac{X_{i,k}^{n+1}X_{i,k}^n}{\mathrm{\Delta }t}},`$ (39)
where the superscript $`n`$ indicates that the quantity is estimated at the time $`t^n`$. The Crank–Nicholson method consists in evaluating the right hand term at the same time as the left hand term, i.e. at time $`t^{n+1/2}`$. The system is centered in time and thus second order accurate in time. One obtains

$`{\displaystyle \frac{X_{i,k}^{n+1}X_{i,k}^n}{\mathrm{\Delta }t}}=`$ (40)
$`{\displaystyle \frac{\sigma _{k+1/2}}{\mathrm{\Delta }m_k}}D_{k+1/2}{\displaystyle \frac{\frac{X_{i,k+1}^{n+1}+X_{i,k+1}^n}{2}\frac{X_{i,k}^{n+1}+X_{i,k}^n}{2}}{r_{k+1}r_k}}`$
$`+{\displaystyle \frac{\sigma _{k1/2}}{\mathrm{\Delta }m_k}}D_{k1/2}{\displaystyle \frac{\frac{X_{i,k}^{n+1}+X_{i,k}^n}{2}\frac{X_{i,k1}^{n+1}+X_{i,k1}^n}{2}}{r_kr_{k1}}},`$
where one has introduced the quantity $`\sigma _{k+1/2}`$ defined by

$$\sigma _{k+1/2}=4\pi r_{k+1/2}^2\rho _{k+1/2}.$$
(41)
Let us define the quantity $`[k+1,k]`$ by

$$[k+1,k]=\frac{1}{2}\sigma _{k+1/2}\frac{D_{k+1/2}\mathrm{\Delta }t}{r_{k+1}r_k},$$
(42)
and let us separate the abundances evaluated at time $`t^n`$ from those estimated at time $`t^{n+1}`$. Eq. (40) then becomes

$`{\displaystyle \frac{[k,k1]}{\mathrm{\Delta }m_k}}X_{i,k1}^{n+1}+\left(1{\displaystyle \frac{[k,k1]}{\mathrm{\Delta }m_k}}{\displaystyle \frac{[k+1,k]}{\mathrm{\Delta }m_k}}\right)X_{i,k}^{n+1}`$
$`+{\displaystyle \frac{[k+1,k]}{\mathrm{\Delta }m_k}}X_{i,k+1}^{n+1}=`$
$`{\displaystyle \frac{[k,k1]}{\mathrm{\Delta }m_k}}X_{i,k1}^n+\left(1+{\displaystyle \frac{[k,k1]}{\mathrm{\Delta }m_k}}+{\displaystyle \frac{[k+1,k]}{\mathrm{\Delta }m_k}}\right)X_{i,k}^n`$
$`{\displaystyle \frac{[k+1,k]}{\mathrm{\Delta }m_k}}X_{i,k+1}^n.`$ (43)
Following a similar line of reasoning, one obtains for the equations at the center and at the surface (see Fig. 1),

$`{\displaystyle \frac{[K,K1]}{\mathrm{\Delta }m_K}}X_{i,K1}^{n+1}+\left(1{\displaystyle \frac{[K,K1]}{\mathrm{\Delta }m_K}}\right)X_{i,K}^{n+1}=`$
$`{\displaystyle \frac{[K,K1]}{\mathrm{\Delta }m_K}}X_{i,K1}^n+\left(1+{\displaystyle \frac{[K,K1]}{\mathrm{\Delta }m_K}}\right)X_{i,K}^n,`$ (44)
$`{\displaystyle \frac{[2,1]}{\mathrm{\Delta }m_1}}X_{i,2}^{n+1}+\left(1{\displaystyle \frac{[2,1]}{\mathrm{\Delta }m_1}}\right)X_{i,1}^{n+1}=`$
$`{\displaystyle \frac{[2,1]}{\mathrm{\Delta }m_1}}X_{i,2}^n+\left(1+{\displaystyle \frac{[2,1]}{\mathrm{\Delta }m_1}}\right)X_{i,1}^n,`$ (45)
where $`\mathrm{\Delta }m_K=_0^{r_{K1/2}}4\pi r^2\rho 𝑑r`$ and $`\mathrm{\Delta }m_1=_{r_{11/2}}^{r_1}4\pi r^2\rho 𝑑r`$. Eqs. (43), (44) and (45) constitute a system of linear equations whose unknowns are the $`X_{i,k}^{n+1}`$. It is solved by using classical methods of tridiagonal matrix inversion.

The discretized equations describing the diffusion of the angular momentum are obtained in a similar way. The final result are equations similar to Eqs. (43) (44) and (45) with $`X_i`$ replaced by $`\mathrm{\Omega }`$, $`\mathrm{\Delta }m_k`$ replaced by $`\mathrm{\Delta }m_kr_k^2`$, $`\sigma _{k+1/2}`$ replaced by $`\sigma _{k+1/2}r_{k+1/2}^2`$ and $`\sigma _{k1/2}`$ replaced by $`\sigma _{k1/2}r_{k1/2}^2`$. Of course the diffusion coefficient must be the one describing the diffusion of the angular momentum.

Let us emphasize here that this method is not fully implicit and thus, one can expect that some sort of Courant’s condition will limit its domain of validity (see Sect. 5).

### 3.3 The implicit finite differences method

In the implicit finite differences method, the right hand term of Eq. (38) is estimated at time $`t^{n+1}`$ (in the explicit method the right hand term would be estimated at time $`t^n`$). In this case, one obtains

$`{\displaystyle \frac{X_{i,k}^{n+1}X_{i,k}^n}{\mathrm{\Delta }t}}=`$ (46)
$`{\displaystyle \frac{\sigma _{k+1/2}}{\mathrm{\Delta }m_k}}D_{k+1/2}{\displaystyle \frac{X_{i,k+1}^{n+1}X_{i,k}^{n+1}}{r_{k+1}r_k}}`$
$`+{\displaystyle \frac{\sigma _{k1/2}}{\mathrm{\Delta }m_k}}D_{k1/2}{\displaystyle \frac{X_{i,k}^{n+1}X_{i,k1}^{n+1}}{r_kr_{k1}}},`$
Let us define the bracket terms $`[k+1,k]`$ by

$$[k+1,k]=\sigma _{k+1/2}\frac{D_{k+1/2}\mathrm{\Delta }t}{r_{k+1}r_k},$$

doing so, we obtain, separating the abundances at time $`t^n`$ from those estimated at time $`t^{n+1}`$,

$`{\displaystyle \frac{[k,k1]}{\mathrm{\Delta }m_k}}X_{i,k1}^{n+1}+\left(1{\displaystyle \frac{[k,k1]}{\mathrm{\Delta }m_k}}{\displaystyle \frac{[k+1,k]}{\mathrm{\Delta }m_k}}\right)X_{i,k}^{n+1}`$ (47)
$`+{\displaystyle \frac{[k+1,k]}{\mathrm{\Delta }m_k}}X_{i,k+1}^{n+1}=X_{i,k}^n.`$
Following a similar line of reasoning, one obtains for the equations at the center and at the surface,

$$\frac{[K,K1]}{\mathrm{\Delta }m_K}X_{i,K1}^{n+1}+\left(1\frac{[K,K1]}{\mathrm{\Delta }m_K}\right)X_{i,K}^{n+1}=X_{i,K}^n,$$
(48)
$$\frac{[2,1]}{\mathrm{\Delta }m_1}X_{i,2}^{n+1}+\left(1\frac{[2,1]}{\mathrm{\Delta }m_1}\right)X_{i,1}^{n+1}=X_{i,1}^n.$$
(49)
One can easily check that the system of discretized equations (both in the cases of the Crank–Nicholson scheme and in the implicit finite differences method) conserves the integrated quantities of the elements over the mass of the star. The system also keeps equal to one the sum of the $`X_i`$ and does not induce any diffusion when the chemical gradients are flat.

For obtaining the equations for the angular momentum, we apply the same recipe as indicated at the end of the previous section. The system of equations describing the transport of the angular momentum also conserves the total angular momentum.

## 4 Estimate of $`D_{k1/2}`$

Let us consider two mesh points as in Fig. 2. Let us suppose that the diffusion coefficient $`D_{k1}`$ is the coefficient for the whole region between $`r_{k1}`$ and $`r_c`$. Similarly the diffusion coefficient $`D_k`$ is the coefficient for the whole region between $`r_c`$ and $`r_k`$. When both mesh points are in a radiative or convective region, the radius $`r_c`$ is equal to $`(r_{k1}+r_k)/2`$, otherwise $`r_c`$ is taken as the radius of the limit between the radiative and the convective zone.

Over the path $`f\mathrm{\Delta }r`$ (see Fig. 2), the particles have an average diffusive velocity $`V_k`$ (see Eq. 35) and over the path $`(1f)\mathrm{\Delta }r`$ they have a diffusive velocity $`<V_{k1}>`$. The total time for going from $`k`$ to $`k1`$ is equal to

$$\mathrm{\Delta }t=f\mathrm{\Delta }r/<V_k>+(1f)\mathrm{\Delta }r/<V_{k1}>.$$
(50)
The mean velocity $`<V_{k1/2}>`$ over the whole interval is given by $`\mathrm{\Delta }r/\mathrm{\Delta }t`$, thus one has

$$<V_{k1/2}>=\frac{1}{\frac{f}{<V_k>}+\frac{(1f)}{<V_{k1}>}}$$

$$=\frac{<V_k><V_{k1}>}{f<V_{k1}>+(1f)<V_k>}.$$

Now, the diffusive coefficient is proportional to the diffusive velocity (see Eq. 35). This suggests the way to compute the diffusive coefficient between two shells:

$$D_{k1/2}=\frac{D_{k1}D_k}{fD_{k1}+(1f)D_k}.$$
(51)
This expression is more physical than

$$D_{k1/2}=(D_{k1}+D_k)/2,$$

that simply averages the diffusion coefficients. Equation (51) implies that if $`D_k>>D_{k1}`$ and $`(1f)`$ is of the order of 0.5 then $`D_{k1/2}D_{k1}`$. As expected, the smallest diffusion coefficient governs the diffusion between the two mesh points. The simple algebraic mean would give $`D_{k1/2}D_k/2`$ a much greater diffusion coefficient, which is not physically justified. One sees also that when $`D_{k1}=D_k=D`$ then $`D_{k1/2}=D`$ whatever the value of $`f`$. In Sect. 5 below we shall illustrate by a numerical example the importance of correctly evaluating the diffusion coefficient at the interface between a convective and a radiative zone.

When the finite elements method is used, we need to evaluate the mean value between two mass shells of the diffusion coefficient, $`D`$, multiplied by a factor $`\alpha `$, which for the transport of the chemical elements is $`(4\pi r^2\rho )^2`$ and for that of the angular momentum is $`r^2(4\pi r^2\rho )^2`$ (see Sect. 3.1). In that case, one adopts the following procedure:

$$\alpha _{k1/2}D_{k1/2}=0.5\left(\alpha _k+\alpha _{k+1}\right)\frac{D_{k1}D_k}{fD_{k1}+(1f)D_k}.$$

## 5 Tests and comparisons

### 5.1 Initial configuration and Courant’s condition

In this section, we study how the different methods described above behave in a very simple configuration. Let us consider a uniform density sphere of 1 M, with a radius of 1 R, composed of two chemical elements $`X`$ and $`Y`$. The initial distribution of these two elements in the sphere is a quasi step function. The variation as a function of the radius of the abundance of one of these two elements, let us call it $`X`$, at the beginning of the computation is shown in Fig. 3 by the dotted line, (cf. also the dotted lines in Figs. 5 to 8). At the center, its initial value is 0.99, in the outer regions $`X=0.01`$. The abundance $`Y`$ is simply $`1X`$. The diffusion coefficient is set equal to $`10^{14}`$ cm2 s-1 in the interior region (i.e. for a radius $`r`$ inferior to 0.7963 R or a mass inferior to 0.5049 M) and to 1.5 106 cm2 s-1 in the outer zone. Such values for the diffusion coefficient imply that the interior has a very small mixing timescale (of the order of $`R^2/D`$ $``$ 1 yr) and therefore will be always homogenized by mixing, while the outer one will have a much longer mixing timescale (of the order of 4.26 Myr). Such a situation is quite analog to a convective core surrounded by a radiative envelope in a star. In the following we shall define the interior region as “the core”, and the outer one as “the envelope”. In our numerical experiments we keep the diffusion coefficient constant as a function of time. The sphere is decomposed in 100 shells of equal mass.

From these data one expects that the whole star will be completely mixed in less than about 107 yr and that the final abundance of element X in the homogeneous star will be 0.5048. From this initial structure, on can also estimate the “Courant condition”. Let us recall that the “Courant condition” imposes a superior limit $`\mathrm{\Delta }t_c`$ to the time step for an explicit method to be stable (see e.g. Press et al. Teu92 (1992)). To estimate this limit, one has to compute for each element and for each mass shell, the time required for the element to diffuse through the width of the shell i.e.

$$\mathrm{\Delta }t_X=\frac{\mathrm{\Delta }r}{v_X}=\frac{\mathrm{\Delta }r}{\frac{D}{X}\left|\frac{X}{r}\right|}.$$

Then one has to take the smallest value of all these times. In a discretized form, $`\mathrm{\Delta }t_c`$ may be written

$$\mathrm{\Delta }t_c=\mathrm{Min}\left(\frac{(r_ir_{i+1})^2}{D_{i+1/2}\left|\frac{X_iX_{i+1}}{\frac{(X_i+X_{i+1})}{2}}\right|}\right).$$

For the simple initial structure considered here, the diffusion time through a shell takes a non infinite value only at the interface between the core and the envelope. In this case, one has that $`r_ir_{i+1}`$ is equal to 0.000199 R, $`D_{i+1/2}`$, estimated from Eq. (51) with $`f`$ taken equal to 0.5, is equal to 3 106 cm2 s, $`X_{i+1}X_i`$ = 0.98 and $`\frac{X_i+X_{i+1}}{2}`$ is equal to 0.5. This gives

$$\mathrm{\Delta }t_c1\mathrm{yr}.$$

Let us stress that implicit methods are generally stable for any time step, and thus are not limited by the Courant condition exposed above (cf. Press et al. Teu92 (1992)). In the following we shall test this point.

### 5.2 Comparisons of the methods

Let us begin by illustrating the importance of correctly estimating the diffusion coefficient at the interface between a convective and a radiative zone. On Fig. 3, the resulting distributions of the element $`X`$ inside the star after 1000 yr is shown. The time step used is $`\mathrm{\Delta }t`$ = 10 yr The continuous line shows the results obtained using Eq. (51) for $`D_{k1/2}`$. The dashed line represent the solution obtained using a simple algebraic mean. We see that the results are significantly different, in particular the algebraic mean tends to slightly increase the convective core. This may have important consequences when such increases are repeatedly applied over the whole evolution of a star. The use of Eq. (51) which results from physical considerations is thus recommended, and this the one we used in the numerical experiments we now describe.

We have computed the diffusion of the chemical elements (and of the angular momentum, see Sect. 5.7) for different durations $`\tau `$ of the period during which the diffusion operates and for different time steps $`\mathrm{\Delta }t`$. In Fig. 4 the set of values ($`\mathrm{log}\mathrm{\Delta }t`$, $`\mathrm{log}\tau `$) explored are indicated. For each couple of values, we performed computations with the three methods described above, namely the implicit finite elements method, the Crank–Nicholson finite differences method and the implicit finite differences method. Most of the results are displayed in Figs. 5 to 8. On these figures, the results obtained for increasing durations, using the same time step, are ordered horizontally (from left to right), while the results obtained for the same duration, with increasing time steps, are arranged vertically (from top to bottom). Since the Courant time step is equal to one year, the time step $`\mathrm{\Delta }t`$ expressed in years gives directly the time step in units of the Courant time step.

### 5.3 The implicit finite elements method

The solutions obtained with the implicit finite elements method are shown in Figs. 5 and 6 (the dashed–dotted lines). This method gives reasonable results in all the cases explored in this work except for small values of $`\tau `$. We note that for $`\tau <1000`$ yr, whatever is the time step, the solution leads to negative values of $`X`$ just above the “convective core”. If the integration is performed over a sufficiently long time, the system eventually reaches a reasonable solution, even when instabilities appear at earlier time (see for instance in Figs. 5 and 6 the evolution when $`\tau `$ increases with $`\mathrm{\Delta }t`$ = 0.5 yr). One notes that the star is completely mixed for $`\tau >10^7`$, i.e. for durations more than twice the mixing timescale of the envelope (about 4.3 Myr, see Sect. 5.1), which is not very satisfactory. The final abundance in the homogeneous star is on the other hand equal to that expected ($`0.5`$).

As expected from an implicit scheme, reasonable solutions are obtained even if the time steps are much greater than the time step given by the Courant condition. Inspection of the results for $`\tau 1000`$ yr show in general a great stability of the solution with respect to the choice of the time step. The solutions obtained with the implicit finite elements method are in general less mixed that those obtained with the implicit finite differences method (see Figs. 7 and 8).

### 5.4 The Crank–Nicholson finite differences method

This method is a kind of compromise between a fully implicit and an explicit scheme. In order to better understand what happens here, let us briefly recall a few generalities about implicit and explicit schemes (see Press et al Teu92 (1992), p. 838–842).

Explicit finite differences schemes are stable only if the time steps satisfy the Courant condition. Such methods are not suited for the computation of the secular evolution of stars. Indeed, we are interested in modeling the evolution of features with spatial scales of the order of the radius of the star $`R`$, much greater than the distance between two mesh points $`\mathrm{\Delta }r`$. If we are limited to time steps satisfying the Courant condition, we will need of the order of $`R^2/(\mathrm{\Delta }r)^2`$ steps before anything interesting begins to happen. We would like instead use timesteps of the order of $`R^2/D`$ or, maybe, for purpose of accuracy, somewhat smaller.

With such great timesteps, it is no long possible to describe accurately what happens at small spatial scales. However at small scale, the differencing scheme must do “something stable, innocuous, and perhaps not too physically unreasonable” write Press et al. (Teu92 (1992)). One possibility is to use a Crank–Nicholson differencing scheme. One of its main property is to let small–spatial–scale features maintain their initial amplitudes. In that case, according to Press et al. (Teu92 (1992)), the evolution of the larger–scale features, in which we are interested in, take place superposed with a kind of “frozen in” (though fluctuating) background of small–scale features. This is what happens in our numerical experiments. Indeed looking at the continuous lines in Figs. 5 and 6, one sees that the amplitude of the initial step in chemical composition is more or less maintained, although with great fluctuations.

The method gives reasonable solutions for not too big time steps. More precisely, physical solutions are obtained when $`\mathrm{\Delta }t`$ in years is inferior to $`10\sqrt{\tau }`$, or when $`\tau `$ in years is inferior to 100 yr (see Fig. 4). This restriction on the time step is reminiscent of a kind of Courant’s condition. This is not so surprising given that the Crank–Nicholson scheme is a mixture of both an implicit (not submitted to Courant condition) and an explicit method (submitted to Courant condition).

Let us finally note that in its domain of validity in the log $`\mathrm{\Delta }t`$ versus log $`\tau `$ plane this method gives the same solution as the implicit finite differences method, and since, it is centered in time, this scheme is second–order accurate in time.

### 5.5 The implicit finite differences method

Another possibility for imposing an “inoccuous” behaviour to the small–scale features is to adopt an implicit method. Such schemes drive small–scale features to their equilibrium form, i.e. imposes that at small scales

$$\frac{}{r}\left(\rho r^2D\frac{X_i}{r}\right)0.$$

This can be seen from Eq. (46) with $`\mathrm{\Delta }t\mathrm{}`$.

The solutions obtained with the implicit finite differences method are shown in Figs. 7 and 8 (continuous lines). The results show in general a great stability of the solution with respect to the choice of the time step. Only when the time step is of the same order of magnitude as $`\tau `$ can we notice some differences (compare for instance the continuous line for $`\tau `$ = 107 yr and $`\mathrm{\Delta }t`$ = 106 yr with the case $`\tau `$ = 107 yr and $`\mathrm{\Delta }t`$ = 107 yr in Fig. 8).

This method leads to a completely mixed star after a time less than 107 in agreement with the estimate made in Sect. 5.1. The final abundance of the element $`X`$ in the homogeneous star is also equal to that expected.

Although, in that case, the differencing scheme is only first–order accurate in time, it has the advantage over the Crank–Nicholson finite differences method and the implicit finite elements method to propose a physical solution (i.e. without negative abundance values) for all the cases investigated here (see Fig. 4). Moreover, as seen just above, it predicts a timescale for the star to be completely homogenized in agreement with the usual analytical estimate (see Sect. 5.1). In that respect this method does appear as the best one.

### 5.6 Conservation of the sum of the abundances

Let us now have a look on the ability of the different schemes to conserve the sum of the abundances (in mass fraction). In each shell and at any time step, one should have that $`X`$+$`Y`$=1. In Table 1, we indicate for all the cases where the three methods give physical solutions, the maximum over the star of the quantity $`\mathrm{\Delta }\chi =1X(k)Y(k)`$ at the end of the computation. As expected we see that in general $`\mathrm{\Delta }\chi `$ becomes greater when $`\tau `$ increases. The Crank–Nicholson finite differences method appears to give, in most of the cases, the smallest values for $`\mathrm{\Delta }\chi `$ (inferior to 10-4 for the cases considered). This is likely related to the fact that this method is second order accurate in time. The implicit finite elements method gives in general the highest values (in the worst case $`\mathrm{\Delta }\chi `$ is of the order of 10-3), while the implicit finite differences method gives in general results inbetween. Thus the implicit finite differences method enables to avoid unphysical solutions and keeps reasonably well the sum of the abundances equal to one.

### 5.7 Conservation of the angular momentum

We have performed similar tests and comparisons for the diffusion of the angular momentum. We started from an initial configuration where the “convective core” defined in Sect 5.1 has an angular velocity equal to 1$``$10-5 sec-1 and the radiative envelope an angular velocity equal to 1$``$10-6 sec-1. The other variables were taken as described in Sect. 5.1. In such situation the Courant time step, which writes

$$\mathrm{\Delta }t_c=\mathrm{Min}\left(\frac{(r_ir_{i+1})^2}{D_{i+1/2}\left|\frac{\mathrm{\Omega }_i\mathrm{\Omega }_{i+1}}{\frac{(\mathrm{\Omega }_i+\mathrm{\Omega }_{i+1})}{2}}\right|}\right),$$

is equal to 1.24 year, not very different from the one we obtained for the diffusion of the chemical species. It is thus not surprising that the results we obtained are qualitatively similar to those presented in Figs. 5 to 8. The zones where unphysical solution are encountered are the same as those presented in Fig. 4. We show in Table 1 the value of the quantity $`\mathrm{\Delta }B`$ which is equal to

$$\mathrm{\Delta }B=\frac{B_{\mathrm{final}}B_{\mathrm{initial}}}{B_{\mathrm{initial}}},$$

where $`B_{\mathrm{initial}}`$, $`B_{\mathrm{final}}`$ are the total angular momentum of the star at the beginning, respectively at the end of the computation. The relative error on the total angular momentum is much greater than the error on the sum of the abundances. This arises because the total angular momentum is an integrated value over the whole system, implying that the errors can accumulate. The sum of the abundances, instead, is evaluated locally in a given shell.

We see that the finite differences method appear to better conserve the total angular momentum. The Crank–Nicholson finite differences method gives the best results, followed by the implicit finite differences method which reaches the same level of accuracy in most cases. Thus the implicit finite differences method that we recommended on the base of the results obtained for the diffusion of the chemical species gives also very satisfactory results for the diffusion of the angular momentum.

## 6 Conclusion

From the above numerical experiments, it appears that the implicit finite differences method seems to be the most robust one, giving physical solutions in all the cases studied here spanning more than nine orders of magnitude in $`\tau `$ and $`\mathrm{\Delta }t`$. It has moreover the following characteristics: 1) it reduces the problem to a first order differential equation, 2) it enables an easy and clear interpretation of what happens physically in the system, 3) it is quite easy to implement in a code, 4) it conserves reasonably well the sum of the abundances at each mesh point as well as the total angular momentum.

On the basis of the new tests and analysis made here, we can recommend this method to resolve the diffusion equation in stellar interiors. Of course many more tests could be performed by changing the initial conditions. We restrained our discussion here on a case with a very sharp gradient in order to test the different numerical methods in some extreme conditions. Adopting an initially shallower gradients will tend to make the things much more easier for all three methods and they would give identical results.

The diffusion coefficient between two mesh points has to be evaluated correctly in order to obtain reliable results. The mean diffusion coefficient is equal neither to the algebraic nor to the geometric mean. Its expression is given in Eq. (51).