         Astron. Astrophys. 317, 193-202 (1997)

## 4. Method of calculation

The previous set of equations shows that the numbers of direct and reverse charge-exchange processes are of the same order, and that, due to the reverse charge-exchange, the calculation of the oxygen atoms distribution depends on the oxygen ions distribution. As a consequence, both distributions have to be evaluated self-consistently.

To do so, in replacement of the two integro-differential Eqs. (1) and (2), we use an iterative method to solve alternatively the four terms of equation system (3) on the one hand (step 1) and the continuity equation for O ions (4) on the other (step2). More specifically, being given an initial oxygen ions distribution, the number density, the velocity and the temperature of the oxygen atoms as well as the mass source for oxygen ions (Eq. (4)) are calculated with a Monte-Carlo method with trajectory splitting, simulating the flow of a Maxwellian population of O atoms through an interface where it interacts with protons, photons and gravity. Then the resulting ions source is used to calculate numerically the new ion density distribution (step 2), which is reinjected in the next Monte-Carlo computation,etc... The process is stopped when both distributions have converged.

The number density of oxygen ions at the inner boundary is not preliminarily fixed and is here a result of the calculation. Numerical experiments with Eq. (4) show that this parameter has a non- negligible influence in only a very small region near the Sun.

The Monte-Carlo method with splitted trajectories used here is similar to the one build by Malama (1991) for the heliospheric interface neutral H - proton coupling. The computational region, shown in Fig. 1, is divided into cells defined by 50 radial distances and 40 angles. Oxygen and hydrogen atom fluxes at the external boundary surface are characterized by maxwellian distribution for the same temperature and bulk velocity as the neutral H and the protons. More specifically, for each velocity vector w making an angle teta with the normal to the boundary surface, the infinitesimal flux out from the surface corresponding to this velocity vector is identical to the infinitesimal flux corresponding to a maxwellian distribution entering the surface. where is the inner normal to the boundary surface, is the Maxwell distribution function of the A-atoms in the unperturbed LISM. A equals O for O atoms and H for H atoms.

During the Monte-Carlo simulation, each O or H atom is followed along its trajectory, which is influenced by gravity (for O atom) or gravity plus radiation pressure (for H atom). If it is photoionized, this atom disappears from the calculation and, for O, a new oxygen ion is created. For O, when it experiences a charge exchange with a proton, a new oxygen ion is created and the O atom disappears from the calculation. For H, when it experiences a charge exchange with an oxygen ion, a new oxygen atom is created, with a velocity selected in the distribution of the local plasma, and its trajectory is computed again. For an H atom, when it experiences a charge exchange with a proton, a new H atom is created, with a velocity selected in the distribution of the local plasma, and its trajectory is computed again. The H atoms are necessary to obtain the inner sources of O atoms, which are mainly due to charge exchange of H atoms with O ions. For the H atoms-protons charge-exchange the following parameters are used:  is the charge exchange cross section, with being the relative atom- proton velocity measured in centimeters per second, , are constants, , (Maher & Tinsley, 1977 ). The photoionization rate on the Earth's orbit equals (Banks, Kockarts, 1973). The ratio of radiation pressure to gravitation is here equal to 0.55. For a good statistical accuracy, splitting procedures are used (Malama, 1991). Density and sources of O atoms are calculated within each cell of a grid.

© European Southern Observatory (ESO) 1997