Consistent cosmological structure formation on all scales in relativistic extensions of MOND

General relativity manifests very similar equations in different regimes, notably in large scale cosmological perturbation theory, non-linear cosmological structure formation, and in weak field galactic dynamics. The same is not necessarily true in alternative gravity theories, in particular those that possess MONDian behaviour ("relativistic extensions" of MOND). In these theories different regimes are typically studied quite separately, sometimes even with the freedom in the theories chosen differently in different regimes. If we wish to properly and fully test complete cosmologies containing MOND against the $\Lambda$CDM paradigm then we need to understand cosmological structure formation on all scales, and do so in a coherent and consistent manner. We propose a method for doing so and apply it to generalised Einstein-Aether theories as a case study. We derive the equations that govern cosmological structure formation on all scales in these theories and show that the same free function (which may contain both Newtonian and MONDian branches) appears in the cosmological background, linear perturbations, and non-linear cosmological structure formation. We show that MONDian behaviour on galactic scales does not necessarily result in MONDian behaviour on cosmological scales, and for MONDian behaviour to arise cosmologically, there will be no modification to the Friedmann equations governing the evolution of the homogeneous cosmological background. We comment on how existing N-body simulations relate to complete and consistent generalised Einstein-Aether cosmologies. The equations derived in this work allow consistent cosmological N-body simulations to be run in these theories whether or not MONDian behaviour manifests on cosmological scales.