  We use high-resolution cosmological zoom-in simulations from the Feedback in
Realistic Environment (FIRE) project to study the galaxy mass-metallicity
relations (MZR) from z=0-6. These simulations include explicit models of the
multi-phase ISM, star formation, and stellar feedback. The simulations cover
halo masses Mhalo=10^9-10^13 Msun and stellar mass Mstar=10^4-10^11 Msun at z=0
and have been shown to produce many observed galaxy properties from z=0-6. For
the first time, our simulations agree reasonably well with the observed
mass-metallicity relations at z=0-3 for a broad range of galaxy masses. We
predict the evolution of the MZR from z=0-6 as
log(Zgas/Zsun)=12+log(O/H)-9.0=0.35[log(Mstar/Msun)-10]+0.93 exp(-0.43 z)-1.05
and log(Zstar/Zsun)=[Fe/H]-0.2=0.40[log(Mstar/Msun)-10]+0.67 exp(-0.50 z)-1.04,
for gas-phase and stellar metallicity, respectively. Our simulations suggest
that the evolution of MZR is associated with the evolution of stellar/gas mass
fractions at different redshifts, indicating the existence of a universal
metallicity relation between stellar mass, gas mass, and metallicities. In our
simulations, galaxies above Mstar=10^6 Msun are able to retain a large fraction
of their metals inside the halo, because metal-rich winds fail to escape
completely and are recycled into the galaxy. This resolves a long-standing
discrepancy between "sub-grid" wind models (and semi-analytic models) and
observations, where common sub-grid models cannot simultaneously reproduce the
MZR and the stellar mass functions.
