Equation of state of Lennard-Jones spheres One way of accounting for van der Waals interactions in fluids is to use the Lennard-Jones potential [*], which has a repulsive term that dominates at short distances which mimics the hard-core repulsion of particles, and a short-range attraction term (due to van der Waals). LJ is usually written as:
$$
V_{LJ} = 4\epsilon \left[(\frac{\sigma}{r})^{12}-(\frac{\sigma}{r})^{6} \right] \tag{1}
$$
where $\sigma$ denotes the sphere diameter, $r$ is the distance between the spheres, and $\epsilon$ is the depth (minimum) of the attraction potential. 
[*]: This effective potential is derived from the consideration of all 3 key dipolar interactions (between permanent (p) and induced (i) ones), namely: Keesom interaction (p-p), Debye (p-i), and London dispersion interaction (i-i).

Questions:


*

*In the context simple spherical fluids, where the only interaction between them is LJ, has the corresponding phase diagram been studied by theory and simulation previously? In contrast to the hard sphere system which is well studied. [1, 2]. In other words, are there known estimates for the equation of state of such systems? Alternatively,  has the phae diagram been studied for instance using Monte Carlo simulations?

*More generally, be it for hard spheres (pure steric repulsion and no attraction) or LJ spheres, what order parameters are usually of interest for tracking the phase transitions that such systems undergo (e.g., liquid-solid)? 

 A: There's a considerable amount of simulation data (both molecular dynamics and Monte Carlo) for the Lennard-Jones potential. 
You should be aware that there are many different versions of the potential. In simulations it is truncated at a maximum "cutoff" distance; $r_\text{cut}=2.5\sigma$ is a common choice. It may or may not be shifted up by a constant amount, so as to make its value zero at the cutoff distance. Other choices apply a distance-dependent correction, so that the potential and its derivative(s) vanish at the cutoff. For any of these choices, it is possible to apply long-range corrections to simulation results, to get an estimate of the equation of state for the true infinite-range potential. Some features of the phase diagram are very sensitive to these choices. For example, the liquid-vapour coexistence curve is noticeably different: the critical temperature is at $k_BT_c/\epsilon\approx 1.32$ for the full potential, but $k_BT_c/\epsilon\approx 1.086$ for the cutoff and shifted potential with $r_\text{cut}=2.5\sigma$.
A fitted equation of state for $r_\text{cut}=2.5\sigma$, based on a wide range of simulation data, has been presented by Thol, Rutkai, Span, et al. Int J Thermophys, 36, 25 (2015), and for the full potential by Thol, Rutkai, Köster, et al. J Phys Chem Ref Data, 45, 023101 (2016). In each case, those authors have thoughtfully provided a computer program which you can download, to generate thermodynamic values for any given density and temperature in the fluid region.
When it comes to theoretical studies of thermodynamic properties and the phase diagram, there is also a considerable body of work. One approach is to predict the structure by mapping the repulsive part of the Lennard-Jones potential onto an effective hard sphere model (referring back to the theories that you mention in your question), and then treat the attractive part of the potential as a perturbation, to predict thermodynamic properties. This is sometimes called the Weeks-Chandler-Andersen approach, after the authors of the original paper. But really, a comprehensive answer cannot be given here: you need to read a book such as Theory of Simple Liquids by J-P Hansen and I R McDonald (4th edition). This covers many different theoretical approaches.
You also asked about order parameters. For the liquid-vapour transition, the density is the important one. For the liquid-solid transition, one can define parameters based on peaks at particular wave-vectors in the solid-state structure factor (which will be large in the solid, and small in the liquid) or bond-order parameters, termed Steinhardt parameters, characterizing the distribution of nearest-neighbour directions in the crystal. Some more details may be found in this online tutorial.
