Electron-electron scattering (EES) affects the shape of the energy distribution function of charge carriers. This effect has to be considered, for instance, in the modeling of hot carrier solar cells, or in the physics-based modeling of hot carrier degradation of semiconductor devices. The widely used Boltzmann transport equation becomes nonlinear if EES is included in the scattering operator. Traditionally, the numerical solution of that nonlinear equation requires additional approximations and comes at high computational cost. In this work we resort to a two-particle formulation. The kinetic equation for the two-particle distribution function is linear and can be solved by Monte Carlo (MC) methods. While deterministic methods suffer from the curse of high dimensionality, the efficiency of the MC method degrades only little when the dimension of the underlying phase space is increased. For the solution of the stationary transport problem a two-particle MC algorithm has been developed. In stationary simulations of bulk silicon no visible effect of EES on the distribution function and consequently on the electron mobility is observed. In stationary device simulations, however, EES causes an enhanced high energy tail relative to the thermal tail. Currently, methods for statistical enhancement at high energies are under development. The transient transport problem is addressed by an ensemble MC algorithm. The mean energy of hot carriers is observed to relax faster in the presence of EES, whereas the cold carriers experience a temporary energy increase.