We study the plasma thermalization and acceleration process in collisionless magnetic reconnection by using a two-dimensional, particle-in-cell numerical simulation and discuss the plasma mixing process of cold lobe ions into the plasma sheet. We find that the plasma thermalization timescale is longer than the dynamic timescale of magnetic reconnection and that non-Maxwellian ion velocity distribution functions produced during the evolution of magnetic reconnection play an important role on the dynamical structure of the plasma sheet. The ion velocity distribution functions are characterized by four class of distributions: (1) anisotropic, high-speed ion beams along the magnetic field line in the plasma sheet boundary layer, (2) two counter-streaming ions along the magnetic field line inside the plasmoid and around the edge of magnetic diffusion region, (3) nongyrotropic, dumbbell-like ions near the X-type region, and (4) the thermal distribution function downstream of slow shocks. We compare the behavior of those ion velocity distribution functions obtained by our reconnection simulation with the distribution functions observed by the Geotail satellite. Most of ion dynamics observed in the Earth's magnetotail can be well understood by our kinetic reconnection simulation. We find that the plasma mixing between the meandering ions accelerated around the X-type region and the cold ions convected directly from the lobe without passing through the X-type region plays a significant role on the formation of non-Maxwellian ions. ¿ 1998 American Geophysical Union