  We develop a new 3+1 dimensional Monte Carlo cascade solving the kinetic
on-shell Boltzmann equations for partons including the inelastic gg <-> ggg
pQCD processes. The back reaction channel is treated -- for the first time --
fully consistently within this scheme. An extended stochastic method is used to
solve the collision integral. The frame dependence and convergency are studied
for a fixed tube with thermal initial conditions. The detailed numerical
analysis shows that the stochastic method is fully covariant and that
convergency is achieved more efficiently than within a standard geometrical
formulation of the collision term, especially for high gluon interaction rates.
The cascade is then applied to simulate parton evolution and to investigate
thermalization of gluons for a central Au+Au collision at RHIC energy. For this
study the initial conditions are assumed to be generated by independent
minijets with p_T > p_0=2 GeV. With that choice it is demonstrated that overall
kinetic equilibration is driven mainly by the inelastic processes and is
achieved on a scale of 1 fm/c. The further evolution of the expanding gluonic
matter in the central region then shows almost an ideal hydrodynamical
behavior. In addition, full chemical equilibration of the gluons follows on a
longer timescale of about 3 fm/c.
