  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 $b_1$ and second-order bias
parameters $b_2$ and $b_{K^2}$. 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 ($L_{\rm box} = 205{\rm
Mpc}/h$). We consider galaxies from the IllustrisTNG simulations, and our main
result is that the $b_2(b_1)$ and $b_{K^2}(b_1)$ 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 $\sigma_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.
