  When studied at finite temperature, Yang-Mills theories in $3+1$ dimensions
display the presence of confinement/deconfinement phase transitions, which are
known to be of first order -- the $SU(2)$ gauge theory being the exception.
Theoretical as well as phenomenological considerations indicate that it is
essential to establish a precise characterisation of these physical systems in
proximity of such phase transitions. We present and test a new method to study
the critical region of parameter space in non-Abelian quantum field theories on
the lattice, based upon the Logarithmic Linear Relaxation (LLR) algorithm. We
apply this method to the $SU(3)$ Yang Mills lattice gauge theory, and perform
extensive calculations with one fixed choice of lattice size. We identify the
critical temperature, and measure interesting physical quantities near the
transition. Among them, we determine the free energy of the model in the
critical region, exposing for the first time its multi-valued nature with a
numerical calculation from first principles, providing this novel evidence in
support of a first order phase transition. This study sets the stage for future
high precision measurements, by demonstrating the potential of the method.
