We have computed ion velocity distributions for the auroral ionosphere with the use of a Monte Carlo computation scheme that includes a detailed description of the interaction between ions and neutrals. This means that both the speed and angular dependence of the interaction are accounted for in a self-consistent manner. We have used the model to deal with several types of interactions all at once and to describe quantitatively how the velocity distribution evolves in the presence of a mixture of neutral constituents. The model also contains an ad hoc procedure to deal with the effect of chemistry, which is important for periods during which N+2 ions are abundant. We have found that the shape of NO+ ion velocity distributions reaches an asymptotic limit for electric fields reaching about 125 mV/m, that the shape of O+ and N+2 ion velocity distributions may be very difficult to model analytically for very large electric field strengths, but that conversely, analytical theories should do an adequate job for fields up to 100 or 150 mV/m for these two ion species. We have, among other things, presented our results in terms of an ion temperature anisotropy, showing how this anisotropy, or the partition parameters that are sometimes used to describe it, change with electric field strength. We have found in the process that the anisotropy of O+ ions is noticeably greater than what has been assumed to be the case in the past. Finally, we have found that the shape of the distribution is sensitive to the presence of a secondary neutral constituent. For example, the presence of just 25% atomic oxygen is sufficient to give the O+ distribution a detectably greater toroidal character. ¿ American Geophysical Union 1992