We use field-level forward models of galaxy clustering and the EFT likelihood formalism to study, for the first time for self-consistently simulated galaxies, the relations between the linear b1 and second-order bias parameters b2 and bK2. The forward models utilize all of the information available in the galaxy distribution up to a given order in perturbation theory, which allows us to infer these bias parameters with high signal-to-noise, even from relatively small volumes (Lbox = 205 Mpc/h). We consider galaxies from the IllustrisTNG simulations, and our main result is that the b2(b1) and bK2(b1) relations obtained from gravity-only simulations for total mass selected objects are broadly preserved for simulated galaxies selected by stellar mass, star formation rate, color and black hole accretion rate. We also find good agreement between the bias relations of the simulated galaxies and a number of recent estimates for observed galaxy samples. The consistency under different galaxy selection criteria suggests that theoretical priors on these bias relations may be used to improve cosmological constraints based on observed galaxy samples. We do identify some small differences between the bias relations in the hydrodynamical and gravity-only simulations, which we show can be linked to the environmental dependence of the relation between galaxy properties and mass. We also show that the EFT likelihood recovers the value of σ8 to percent-level from various galaxy samples (including splits by color and star formation rate) and after marginalizing over 8 bias parameters. This demonstration using simulated galaxies adds to previous works based on halos as tracers, and strengthens further the potential of forward models to infer cosmology from galaxy data.
