In this work, I investigate the properties of Lyman limit systems (LLSs) using state-of-the-art zoom-in cosmological galaxy formation simulations with on the fly radiative transfer, which includes both the cosmic UV background (UVB) and local stellar sources. I compare the simulation results to observations of the incidence frequency of LLSs and the HI column density distribution function over the redshift range $z=2-5$ and find good agreement. I explore the connection between LLSs and their host halos and find that LLSs reside in halos with a wide range of halo masses with a nearly constant covering fraction within a virial radius. Over the range $z=2-5$, I find that more than half of the LLSs reside in halos with , indicating that absorption line studies of LLSs can probe these low-mass galaxies which H$2$-based star formation models predict to have very little star formation. I study the physical state of individual LLSs and test a simple model (Schaye 2001) which encapsulates many of their properties. I confirm that LLSs have a characteristic absorption length given by the Jeans length and that they are in photoionization equilibrium at low column densities. Finally, I investigate the self-shielding of LLSs to the UVB and explore how the non-sphericity of LLSs affects the photoionization rate at a given $NHI$. I find that at $z\approx 3$, LLSs have an optical depth of unity at a column density of $\sim 1018cm-2$ and that this is the column density which characterizes the onset of self-shielding.