  According to the separate universe conjecture, spherically symmetric
sub-regions in an isotropic universe behave like mini-universes with their own
cosmological parameters. This is an excellent approximation in both Newtonian
and general relativistic theories. We estimate local expansion rates for a
large number of such regions, and use a scale parameter calculated from the
volume-averaged increments of local scale parameters at each time step in an
otherwise standard cosmological $N$-body simulation. The particle mass,
corresponding to a coarse graining scale, is an adjustable parameter. This mean
field approximation neglects tidal forces and boundary effects, but it is the
first step towards a non-perturbative statistical estimation of the effect of
non-linear evolution of structure on the expansion rate. Using our algorithm, a
simulation with an initial $\Omega_m=1$ Einstein--de~Sitter setting closely
tracks the expansion and structure growth history of the $\Lambda$CDM
cosmology. Due to small but characteristic differences, our model can be
distinguished from the $\Lambda$CDM model by future precision observations.
Moreover, our model can resolve the emerging tension between local Hubble
constant measurements and the Planck best-fitting cosmology. Further
improvements to the simulation are necessary to investigate light propagation
and confirm full consistency with cosmic microwave background observations.
