  On very large scales, density fluctuations in the Universe are small,
suggesting a perturbative model for large-scale clustering of galaxies (or
other dark matter tracers), in which the galaxy density is written as a Taylor
series in the local mass density, delta, with the unknown coefficients in the
series treated as free "bias" parameters. We extend this model to include
dependence of the galaxy density on the local values of nabla_i nabla_j phi and
nabla_i v_j, where phi is the potential and v is the peculiar velocity. We show
that only two new free parameters are needed to model the power spectrum and
bispectrum up to 4th order in the initial density perturbations, once symmetry
considerations and equivalences between possible terms are accounted for. One
of the new parameters is a bias multiplying s_ij s_ji, where s_ij=[nabla_i
nabla_j \nabla^-2 - 1/3 delta^K_ij] delta. The other multiplies s_ij t_ji,
where t_ij=[nabla_i nabla_j nabla^-2 - 1/3 delta^K_ij](theta-delta), with
theta=-(a H dlnD/dlna)^-1 nabla_i v_i. (There are other, observationally
equivalent, ways to write the two terms, e.g., using theta-delta instead of
s_ij s_ji.) We show how short-range (non-gravitational) non-locality can be
included through a controlled series of higher derivative terms, starting with
R^2 nabla^2 delta, where R is the scale of non-locality (this term will be a
small correction as long as k^2 R^2 is small, where k is the observed
wavenumber). We suggest that there will be much more information in future huge
redshift surveys in the range of scales where beyond-linear perturbation theory
is both necessary and sufficient than in the fully linear regime.
