Evolution of the virial ratio $-K/U$ in $N$-body simulation

Systems that are capable of reaching equilibrium will do so if given enough time. This is a curious question; I guess you mean what is the mechanism by which gravitational potential energy and kinetic energy are redistributed? The answer to that is mainly two-body gravitational interactions

You started your system way out of equilibrium. What happens is that it collapses, gives lots of kinetic energy to the particles and "overshoots" the equilibrium (a bit like a lightly damped suspension spring). After that it returns towards equilibrium, though looks like it does a further mild subvirial expansion before settling down. I would guess you would find that the timescale for this oscillatory behaviour is the "crossing time" (roughly the freefall timescale in this case), although it settles into an equilibrium suspiciously quickly - for 500 particles one expects a relaxation time of $sim 10$ crossing times (NB very unrealistic for a galaxy).

What do you mean by "the system"? If you mean the remaining bound component then yes it should stay in equilibrium and roughly at constant total energy (since the particles that evaporate have an energy of $sim 0$). But maybe what you are asking about is gravothermal collapse. Since the total energy is some fixed fraction of $M^2/R$, then as mass evaporates, the radius must decrease as $M^{-2}$ and the density goes up as $M^5$. I would think that the continual contraction makes the bound part slightly subvirial, but then after 10-20 relaxation times you will end up with some numerical catastrophe at the centre of the simulation as the density heads towards infinity.