  We study the nature of melting of a two dimensional (2D) Lennard-Jones solid
using large scale Monte Carlo simulation. We use systems of up to 102,400
particles to capture the decay of the correlation functions associated with
translational order (TO) as well as the bond-orientational (BO) order. We study
the role of dislocations and disclinations and their distribution functions. We
computed the temperature dependence of the second moment of the two order
parameters, i.e., of the TO order parameter as well as of the order parameter
associated with BO order and using finite-size scaling we determined the
corresponding two anomalous dimension critical exponents eta. We also computed
the temperature dependent distribution of these two order parameters on the
complex plane which support a two stage melting with a hexatic phase as an
intermediate phase. The analysis of our results leads to a consistent picture
strongly supporting a two stage melting scenario as predicted by the
Kosterlitz, Thouless, Halperin, Nelson, and Young (KTHNY) theory. We find that
the two correlation lengths associated with translational and
bond-orientational order have a distinctly different temperature dependence
each diverging at different temperature and that their finite size scaling
properties are consistent with the KTHNY theory. We also used the temperature
dependence of etas and their theoretical bounds to provide estimates for the
two critical temperatures, which are found to be in agreement with those
determined from the divergent correlation lengths and from the Binder ratio.
