Combination of the power spectrum and bispectrum is a powerful way of breaking degeneracies between galaxy bias and cosmological parameters, enabling us to maximize the constraining power from galaxy surveys. Recent cosmological constraints treat the power spectrum and bispectrum on an uneven footing: they include one-loop bias corrections for the power spectrum but not the bispectrum. To bridge this gap, we develop the galaxy bias description up to fourth order in perturbation theory, conveniently expressed through a basis of Galilean invariants that clearly split contributions that are local and nonlocal in the second derivatives of the linear gravitational potential. In addition, we consider relevant contributions from short-range nonlocality (higher-derivative terms), stress-tensor corrections and stochasticity. To sidestep the usual renormalization of bias parameters that complicates predictions beyond leading order, we recast the bias expansion in terms of multipoint propagators, which take a simple form in our split-basis with loop corrections depending only on bias parameters corresponding to nonlocal operators. We show how to take advantage of Galilean invariance to compute the time evolution of bias and present results for the fourth-order parameters for the first time. We also discuss the possibilities of reducing the bias parameter space by using the evolution of bias and exploiting degeneracies between various bias contributions in the large-scale limit. Our baseline model allows to verify these simplifications for any application to large-scale structure data sets.

36 pages, 6 figures; added two figures and various references, some minor changes throughout