  We use cosmological hydrodynamic simulations with stellar feedback from the
FIRE project to study the physical nature of Lyman limit systems (LLSs) at z<1.
At these low redshifts, LLSs are closely associated with dense gas structures
surrounding galaxies, such as galactic winds, dwarf satellites, and cool
inflows from the intergalactic medium. Our analysis is based on 14 zoom-in
simulations covering the halo mass range M_h~10^9-10^13 Msun at z=0, which we
convolve with the dark matter halo mass function to produce cosmological
statistics. We find that the majority of cosmologically-selected LLSs are
associated with halos in the mass range 10^10 < M_h < 10^12 Msun. The incidence
and HI column density distribution of simulated absorbers with columns 10^16.2
< N_HI < 2x10^20 cm^-2 are consistent with observations. High-velocity outflows
(with radial velocity exceeding the halo circular velocity by a factor >~2)
tend to have higher metallicities ([X/H] ~ -0.5) while very low metallicity
([X/H] < -2) LLSs are typically associated with gas infalling from the
intergalactic medium. However, most LLSs occupy an intermediate region in
metallicity-radial velocity space, for which there is no clear trend between
metallicity and radial kinematics. Metal-enriched inflows arise in the FIRE
simulations as a result of galactic winds that fall back onto galaxies at low
redshift. The overall simulated LLS metallicity distribution has a mean
(standard deviation) [X/H] = -0.9 (0.4) and does not show significant evidence
for bimodality, in contrast to recent observational studies but consistent with
LLSs arising from halos with a broad range of masses and metallicities.
