  The next generation of cosmological surveys will operate over unprecedented
scales, and will therefore provide exciting new opportunities for testing
general relativity. The standard method for modelling the structures that these
surveys will observe is to use cosmological perturbation theory for linear
structures on horizon-sized scales, and Newtonian gravity for non-linear
structures on much smaller scales. We propose a two-parameter formalism that
generalizes this approach, thereby allowing interactions between large and
small scales to be studied in a self-consistent and well-defined way. This uses
both post-Newtonian gravity and cosmological perturbation theory, and can be
used to model realistic cosmological scenarios including matter, radiation and
a cosmological constant. We find that the resulting field equations can be
written as a hierarchical set of perturbation equations. At leading-order,
these equations allow us to recover a standard set of Friedmann equations, as
well as a Newton-Poisson equation for the inhomogeneous part of the Newtonian
energy density in an expanding background. For the perturbations in the
large-scale cosmology, however, we find that the field equations are sourced by
both non-linear and mode-mixing terms, due to the existence of small-scale
structures. These extra terms should be expected to give rise to new
gravitational effects, through the mixing of gravitational modes on small and
large scales - effects that are beyond the scope of standard linear
cosmological perturbation theory. We expect our formalism to be useful for
accurately modelling gravitational physics in universes that contain non-linear
structures, and for investigating the effects of non-linear gravity in the era
of ultra-large-scale surveys.
