Systems with long-range (LR) forces, for which the interaction potential decays with the interparticle distance with an exponent smaller than the dimensionality of the embedding space, remain an outstanding challenge to statistical physics. The internal energy of such systems lacks extensivity and additivity. Although the extensivity can be restored by scaling the interaction potential with the number of particles, the non-additivity still remains. Lack of additivity leads to inequivalence of statistical ensembles. Before relaxing to thermodynamic equilibrium, isolated systems with LR forces become trapped in out-of-equilibrium quasi-stationary state (qSS), the lifetime of which diverges with the number of particles. Therefore, in thermodynamic limit LR systems will not relax to equilibrium. The qSSs are attained through the process of collisionless relaxation. Density oscillations lead to particle-wave interactions and excitation of parametric resonances. The resonant particles escape from the main cluster to form a tenuous halo. Simultaneously, this cools down the core of the distribution and dampens out the oscillations. When all the oscillations die out the ergodicity is broken and a qSS is born. In this report, we will review a theory which allows us to quantitatively predict the particle distribution in the qSS. The theory is applied to various LR interacting systems, ranging from plasmas to self-gravitating clusters and kinetic spin models.

We study, using both theory and molecular dynamics simulations, the relaxation dynamics of a microcanonical two dimensional self-gravitating system. After a sufficiently large time, a gravitational cluster of N particles relaxes to the Maxwell-Boltzmann distribution. The time to reach the thermodynamic equilibrium, however, scales with the number of particles. In the thermodynamic limit, $N\to\infty$ at fixed total mass, equilibrium state is never reached and the system becomes trapped in a non-ergodic stationary state. An analytical theory is presented which allows us to quantitatively described this final stationary state, without any adjustable parameters.

Theory and simulations are used to study collisionless relaxation of a gravitational $N$-body system. It is shown that when the initial one particle distribution function satisfies the virial condition -- potential energy is minus twice the kinetic energy -- the system quickly relaxes to a metastable state described {\it quantitatively} by the Lynden-Bell distribution with a cutoff. If the initial distribution function does not meet the virial requirement, the system undergoes violent oscillations, resulting in a partial evaporation of mass. The leftover particles phase separate into a core-halo structure. The theory presented allows us to quantitatively predict the amount and the distribution of mass left in the central core, without any adjustable parameters. On a longer time scale $\tau_G \sim N$ collisionless relaxation leads to a gravothermal collapse.