The KBC void and Hubble tension contradict $\Lambda$CDM on a Gpc scale $-$ Milgromian dynamics as a possible solution

The KBC void is a local underdensity with the observed relative density contrast $\delta \equiv 1 - \rho/\rho_{0} = 0.46 \pm 0.06$ between 40 and 300 Mpc around the Local Group. If mass is conserved in the Universe, such a void could explain the $5.3\sigma$ Hubble tension. However, the MXXL simulation shows that the KBC void causes $6.04\sigma$ tension with standard cosmology ($\Lambda$CDM). Combined with the Hubble tension, $\Lambda$CDM is ruled out at $7.09\sigma$ confidence. Consequently, the density and velocity distribution on Gpc scales suggest a long-range modification to gravity. In this context, we consider a cosmological MOND model supplemented with $11 \, \rm{eV}/c^{2}$ sterile neutrinos. We explain why this $\nu$HDM model has a nearly standard expansion history, primordial abundances of light elements, and cosmic microwave background (CMB) anisotropies. In MOND, structure growth is self-regulated by external fields from surrounding structures. We constrain our model parameters with the KBC void density profile, the local Hubble and deceleration parameters derived jointly from supernovae at redshifts $0.023 - 0.15$, time delays in strong lensing systems, and the Local Group velocity relative to the CMB. Our best-fitting model simultaneously explains these observables at the $1.14\%$ confidence level (${2.53 \sigma}$ tension) if the void is embedded in a time-independent external field of ${0.055 \, a_{_0}}$. Thus, we show for the first time that the KBC void can naturally resolve the Hubble tension in Milgromian dynamics. Given the many successful a priori MOND predictions on galaxy scales that are difficult to reconcile with $\Lambda$CDM, Milgromian dynamics supplemented by $11 \, \rm{eV}/c^{2}$ sterile neutrinos may provide a more holistic explanation for astronomical observations across all scales.


INTRODUCTION
The Cosmological Principle (CP) states that the Universe is homogeneous and isotropic on very large scales. This concept is the foundation of the current Lambda-Cold Dark Matter (ΛCDM) standard model of cosmology (Ostriker & Steinhardt 1995), which assumes that Einstein's General Relativity is valid on all astrophysical scales. Applying it to the non-relativistic outskirts of galaxies yields nearly the same result as Newtonian dynamics − the rotation curve Email: mhaslbauer@astro.uni-bonn.de mhaslbauer@mpifr-bonn.mpg.de (Moritz Haslbauer) ibanik@astro.uni-bonn.de (Indranil Banik) † Alexander von Humboldt Fellow should undergo a Keplerian decline beyond the extent of the luminous matter (de Almeida et al. 2016). The observed flat rotation curves of galaxies (e.g. Babcock 1939;Rubin & Ford 1970;Rogstad & Shostak 1972) demonstrate that Newtonian gravity of the baryons alone is insufficient to hold them together, leading to the concept that each galaxy is surrounded by a CDM halo (Ostriker & Peebles 1973). However, no experiment has ever confirmed the existence of CDM, with stringent upper limits coming from e.g. null detection of γ-rays from DM annihilation in dwarf satellites of the Milky Way (MW; Hoof et al. 2020). In addition to the hypothetical ingredient of CDM, the ΛCDM model also requires a cosmological constant Λ in Einstein's gravitational field equations to explain the anomalous faintness of distant Type Ia supernovae (SNe Ia; Riess et al. 1998;Schmidt et al. 1998;The Supernova Cosmology Project 1999). Λ may be associated to a vacuum energy (dark energy).
This 'concordance' flat ΛCDM model explains the cosmic microwave background (CMB) as relic radiation from the Universe at redshift z ≈ 1100 (e.g. Bennett et al. 2003;Planck Collaboration VI 2020). The temperature fluctuations within the CMB are of the order δT/T ≈ δρ/ρ ≈ 10 −5 (Wright 2004). These are interpreted as tracers of density contrasts in the baryons alone, with the CDM being significantly more clustered by that time due to it not feeling radiation pressure. After recombination, baryons fell into the potential wells of the DM, starting the process of cosmic structure formation via gravitational instability.
Observations have shown that this widely used ΛCDM model faces several challenges, especially on galactic up to Mpc scales (e.g. Kroupa 2012Kroupa , 2015, and references therein). One of the most serious problems is the distribution of dwarf galaxies in the Local Group (LG). The MW is surrounded by a thin co-rotating disc of satellite galaxies (Kroupa et al. 2005), which is part of the vast polar structure (Pawlowski et al. 2012) that also includes ultra-faint galaxies, globular clusters, and gas and stellar streams. Recently, Pawlowski & Kroupa (2020) showed that its kinematic coherence has increased further with Gaia Data Release 2 (Gaia Collaboration 2018). A thin plane of co-rotating satellites is also observed around M31 (Ibata et al. 2013).
It is very difficult to understand these structures if their member satellites are primordial . However, such phase space-correlated structures can arise during an interaction between two disc galaxies, as observed e.g. in the Antennae galaxies (Mirabel et al. 1992). Due to the higher velocity dispersion of the DM, such tidal dwarf galaxies (TDGs) should be free of DM in ΛCDM, as shown with simulations of galaxy interactions (Barnes & Hernquist 1992;Wetzstein et al. 2007) and in cosmological simulations (Ploeckinger et al. 2018;Haslbauer et al. 2019b). This would lead to very low internal velocity dispersions, which are in conflict with observations for satellites of the MW (McGaugh & Wolf 2010) and M31 (McGaugh & Milgrom 2013).
A disc of satellites has also been observed around Centaurus A (Cen A; Müller et al. 2018), suggesting that such structures are ubiquitous and in any case not unique to the LG. Although they may well consist of TDGs, these are quite rare in ΛCDM due to their weak Newtonian selfgravity (Haslbauer et al. 2019a,b). This makes the Cen A satellite plane hard to explain even though we lack internal velocity dispersion measurements for its members (Müller et al. 2018). A review on satellite planes in the local Universe can be found in Pawlowski (2018), who suggested that the TDG hypothesis could work in an alternative gravitational framework where all galaxies are DM-free. We consider this possibility further in Section 1.3. Some evidence in favour of this scenario is the strong correlation between the bulge fractions and the number of satellite galaxies for the MW, M31, M81, Cen A, and M101 (Javanmardi & Kroupa 2020). This is unexpected in standard cosmology (Kroupa 2012(Kroupa , 2015Javanmardi et al. 2019), but may indicate that bulges and satellite galaxies formed simultaneously in galactic interactions.
Although ΛCDM is widely considered a successful theory in explaining large-scale structure, the observed Universe appears to be much more structured and organized than it predicts. In particular, Peebles & Nusser (2010) reported that standard ΛCDM theory is in conflict with the distribution of galaxies within ≈ 8 Mpc of the LG. The local void contains much fewer galaxies than expected (e.g. Tikhonov & Klypin 2009), while massive galaxies are located away from the matter sheets where they ought to reside. These facts suggest a more rapid growth rate of structure (Peebles & Nusser 2010; though see Xie et al. 2014). Karachentsev (2012) studied the matter distribution of the Local Volume in more detail, finding that the average density of matter within ≈ 50 Mpc is only Ω m,loc = 0.08±0.02, much lower than the global cosmic density at the present time (Ω m,0 = 0.315, Planck Collaboration VI 2020). This is consistent with a more recent work which obtained Ω m,loc = 0.09 − 0.14 within a sphere of radius 40 Mpc around the LG (Karachentsev & Telikova 2018). This is striking because the Harrison-Zeldovich spectrum and the current value of σ 8 = 0.811 ± 0.006 (Planck Collaboration VI 2020) predict root mean square (rms) density fluctuations of 23% on this scale. Indeed, recent studies have questioned the assumption of homogeneity and isotropy (e.g. Javanmardi et al. 2015;Kroupa 2015;Javanmardi & Kroupa 2017;Bengaly et al. 2018;Colin et al. 2019;Mészáros 2019;Migkas et al. 2020).
Therefore, observations of the galaxy distribution on large scales can constrain various cosmological models and their different underlying gravitational theories. In this study, we investigate the local matter density and velocity field within 1 Gpc in ΛCDM and in a previously developed Milgromian cosmological model (Angus 2009). This allows us to assess the implications for the CP and Hubble tension.

KBC void
Several observations at different wavelengths have found evidence for a large local underdensity around the LG. The first indication for a deficiency in the galaxy luminosity density was observed in optical samples (e.g. Maddox et al. 1990). Using the ESO Slice Project galaxy survey that covers ≈ 23 deg 2 on the sky, Zucca et al. (1997) found a local underdensity out to a distance of ≈ 140h −1 Mpc in the b J band, where h ≈ 0.7 is the present Hubble constant H 0 in units of 100 km s −1 Mpc −1 .
Galaxy counts in the near-infrared (NIR) revealed that the local Universe is significantly underdense on a scale of 200 − 300h −1 Mpc around the LG (e.g. Huang et al. 1997;Frith et al. 2003;Busswell et al. 2004;Frith et al. 2005Frith et al. , 2006Keenan et al. 2013;Whitbourn & Shanks 2014). NIR photometry accurately traces the stellar mass and is therefore a good proxy for the underlying matter distribution.
A local underdensity is also evident in the X-ray galaxy cluster surveys REFLEX II (Böhringer et al. 2015) and CLASSIX (Böhringer et al. 2020). The latter work found a 15 − 30% (10 − 20%) underdensity in the matter distribution within a radius of ≈ 100 Mpc (140 Mpc).
At the opposite end of the spectrum, Rubart & Schwarz (2013) found that the cosmic radio dipole from the NRAO VLA Sky Survey is ≈ 4× stronger than can be explained purely kinematically given the magnitude of the CMB dipole. Interestingly, the radio dipole points towards Galactic coordinates (245 • , +43 • ) which, given the uncertainty of ≈ 30 • , is consistent with the direction in which the LG moves with respect to (wrt.) the CMB (276 • ± 3 • , +30 • ± 3 • ; Kogut et al. 1993). In a subsequent study, Rubart et al. (2014) showed that the unusually strong radio dipole could be explained by a single void with a size of 11% of the Hubble distance and a density contrast of δ ≡ 1 − ρ/ρ 0 = 1/3, where ρ is the local density and ρ 0 is the cosmic mean.
Moreover, Bengaly et al. (2018) studied the dipole anisotropy of galaxy number counts over the redshift range 0.10 < z < 0.35, revealing a large anisotropy for z < 0.15 that could be the imprint of a large local density fluctuation. Thus, a significant local underdensity is evident across the entire electromagnetic spectrum.
Here, we focus on the study by Keenan, Barger & Cowie (2013), who found clear evidence for a large local underdensity by measuring the K-band galaxy luminosity function at different distances over a large part of the sky (see their figures 9 and 10). They used the 2M++ catalogue (Lavaux & Hudson 2011), which combines photometry from the Two Micron All Sky Survey Extended Source Catalog (2MASS-XSC) with redshifts from the Sloan Digital Sky Survey (SDSS), the Two Micron Redshift Survey (2MRS), and the Six-degree Field Galaxy Redshift Survey (6DFGRS). This sample covers 37 080 deg 2 (90% of the whole sky) and is ≈ 98% complete to a limiting magnitude of K s = 13.36. Using this sample, Keenan et al. (2013) estimated the luminosity density and derived a relative density contrast of δ ≈ 0.5 in the redshift range 0.0025 < z < 0.067 compared to larger redshifts (see the pink down-pointing triangle in their figure 11, and their table 1). In addition, they also probed the density field to a deeper magnitude limit of K s = 14.36, but only in the SDSS and 6DFGRS regions. This yielded a slightly smaller density contrast of δ = 0.46 ± 0.06 between z = 0.01 (≈ 40 Mpc) and z = 0.07 (≈ 300 Mpc; see the light blue dot in their figure 11). In the following, we will show that the Keenan-Barger-Cowie (KBC) void is highly unexpected within the ΛCDM framework by virtue of its sheer size and depth. In order to minimize the tension, we assume for our analysis that δ = 0.46 ± 0.06, and refer to this as the KBC void. Calculating the K-band luminosity density in different regions suggests that it reaches the cosmic mean at a distance of ≈ 500 Mpc.
In the ΛCDM framework, the existence of such a deep and extended void is a puzzle given the expected Harrison-Zeldovich scale-invariant power spectrum, which states that the power P (L) on some length scale L varies as P(L) ∝ L −n s , with n s = 1 (Harrison 1970;Zeldovich 1972). Since the CMB anisotropies require a power of σ 8 = 0.811 ± 0.006 on a scale of 8h −1 Mpc (Planck Collaboration VI 2020), we expect density fluctuations of only ≈ 3.2% between spheres of radius L = 300 Mpc.
Combining measurement errors with cosmic variance, we can estimate that the KBC void would falsify the ΛCDM model by well over 5σ because 0.46 √ 0.06 2 + 0.032 2 = 6.8 .
In Section 2, we provide a much more sophisticated analysis of how likely the KBC void is in standard cosmology. Since the measurement uncertainty of 6% is much larger than the cosmic variance of 3.2%, the latter is not the main source of uncertainty in how far off ΛCDM is from matching the observations − as explicitly calculated in Section 2.2.1. Con-sequently, if we assume that ΛCDM is the correct model, the most likely explanation for the detection of such a deep void would be a measurement error. However, the KBC void is evident over the entire electromagnetic spectrum.
The above prediction of 3.2% rests on two fundamental assumptions − that the CMB reflects baryonic density fluctuations at z = 1100, and that General Relativity is valid on all scales. The existence of the KBC void might indicate that either or both of these assumptions must be relaxed. In this contribution, we focus on modifying gravity because the standard approach leads to problems in galaxies (e.g. Kroupa 2012, 2015, andreferences therein).
A large local void should also have implications for local measurements of cosmological parameters such as the Hubble constant and deceleration parameter. If mass is conserved in the Universe and it was nearly homogeneous initially, a large fractional underdensity would show up in the velocity field. This is because the co-moving radius enclosing a fixed amount of mass must exceed its initial value, and changes in co-moving coordinates imply a peculiar velocity.
Suppose that we are living near the centre of a void whose true density relative to the cosmic mean is This implies that the co-moving radius enclosing a fixed mass must exceed its initial value by a factor α −1/3 . Depending on details of how the void grows, the impact on the locally measured Hubble parameter would be approximately the same. In other words, where H local 0 is the locally measured H 0 , whose background (true) value is H global 0 ≡ a/a at the present time, with a the cosmic scale factor and an overdot indicating a time derivative. The mismatch between these H 0 values would create a redshift space distortion (RSD) effect whereby the physical volume of a survey with known redshift range would be reduced by a factor α compared to the case of no void. In this way, RSD would further reduce the observed α obs by a factor α if it is not accounted for and a constant H 0 is used to convert redshifts to distances (as done in the work of Keenan et al. 2013, see their section 4.7). Thus, we expect that Combining Equations 3 and 4, we get that Given that α obs = 0.54, the measured H local 0 should exceed the background value H global 0 by 0.54 −1/6 , i.e. by 11%. This would raise H 0 from the Planck-based prediction of 67.4 km s −1 Mpc −1 (Planck Collaboration VI 2020) to 74.7 km s −1 Mpc −1 , which is very close to the observed value (Section 1.2). This is unlikely to be a coincidence − it is more parsimoniously explained as a consequence of the observed void under the standard assumption of matter conservation.

Hubble tension
In this context, we consider the Hubble tension, a statistically significant discrepancy between the locally measured cosmic expansion rate and the ΛCDM prediction based on the early universe properties needed to match the CMB power spectrum (e.g. Riess 2019). The local Hubble constant can be determined through the distance ladder technique. Recently, the Supernova H 0 for the Equation of State (SH0ES) team ) calibrated the distance ladder with eclipsing binaries in the Large Magellanic Cloud, masers in NGC 4258, and parallaxes of Galactic Cepheid variables via the Leavitt law. They derived a local Hubble constant of H local 0 = 74.03 ± 1.42 km s −1 Mpc −1 , which results in 4.4σ tension with the Planck-based prediction (H Planck 0 = 67.4 ± 0.5 km s −1 Mpc −1 ; Planck Collaboration VI 2020).
The systematic error of the Cepheid background subtraction is only 0.029 ± 0.037 mag, which is not sufficient to explain the ≈ 0.2 mag Hubble tension . Moreover, calibrating the SN Ia luminosity using instead Mira variables in the galaxy NGC 1559 with periods of 240 − 400 d and using NGC 4258 (the Large Magellanic Cloud) as an anchor, Huang et al. (2020) obtained H local 0 = 72.7 ± 4.6 km s −1 Mpc −1 (H local 0 = 73.9 ± 4.3 km s −1 Mpc −1 ; see also their table 6 and figure 11). Both values are consistent with H local 0 derived from Cepheid variables within the 1σ confidence range, though the Mira-calibrated H 0 is less precise.
It is also possible to go beyond the traditional Cepheid-SN Ia route using Type II SNe as standard candles. These yield a high H local 0 of 75.8 +5.2 −4.9 km s −1 Mpc −1 , which is very consistent with H 0 derived from Type Ia SNe − albeit with larger uncertainties (de Jaeger et al. 2020). Thus, systematic errors in Type Ia SNe data are likely not driving the Hubble tension. Camarena & Marra (2020a) analysed the Pantheon SNe Ia sample without fixing the deceleration parameter (q ≡ −a a/ a 2 ) to the present ΛCDM prediction of q 0 = −0.55. They jointly derived H local 0 = 75.35 ± 1.68 km s −1 Mpc −1 and q 0 = −1.08 ± 0.29 from SNe in the redshift range 0.023 ≤ z ≤ 0.15. This is in 4.54σ tension with ΛCDM. The unexpectedly low q 0 is robust to the choice of data set (table 5 of Camarena & Marra 2020b).
Interestingly, it is highly implausible to get such low q 0 values at the background level. Even in a pure dark energydominated (de Sitter) universe, it is not possible to get q 0 < −1. Thus, first-and second-order effects in the local Hubble diagram seem to provide additional evidence for the KBC void. To quantify this, we compare the standard expansion rate history (H 0 = 67.4 km s −1 Mpc −1 , q 0 = −0.55) with an extrapolation of the Camarena & Marra (2020a) results. Approximating both as quadratic functions of time t with a(t 0 ) ≡ 1 at the present time t 0 , we get that the reconstructed a(t) parabolas coincide 4.2 Gyr ago. This provides strong evidence for a Gpc-scale void independently of the galaxy luminosity density (discussed earlier in Section 1.1).
A method of measuring H 0 independently of the cosmic distance ladder relies on time delays between multiple images of the same source, as occurs in strong gravitational lensing. Jee et al. (2019)  −3.0 km s −1 Mpc −1 from the strong lens system DES J0408−5354, whose deflector lies at an angular diameter distance of D d = 1711 +376 −280 Mpc (z = 0.597). This is broadly consistent with measurements of the H 0 Lenses in COSMO-GRAIL's Wellspring (H0LiCOW; Wong et al. 2020). Using a blinded analysis protocol (see their section 3.6), they obtained H 0 = 73.3 +1.7 −1.8 km s −1 Mpc −1 from six lensed quasar systems in the redshift range z = 0.295 − 0.745. Combining their results with the measurement of Riess et al. (2019) leads to a 5.3σ discrepancy with ΛCDM expectations based on the CMB (Wong et al. 2020). The latter work showed for the first time that the Hubble tension exceeds the 5σ threshold typically used to judge the validity of scientific theories.
Although Kochanek (2020) suggested there might be biases in the strong lensing analysis causing ≈ 10% uncertainties on the inferred H 0 , Pandey et al. (2020) showed that the SNe and strong lensing measurements are consistent and likely have systematics much smaller than the Hubble tension, as also found by Millon et al. (2020). Indeed, the near-perfect agreement between the SNe and lensing determinations despite the blinded protocol of the latter does suggest rather small uncertainties. Moreover, Wong et al. (2020) found that H 0 measured from strong lensing decreases as a function of lens redshift at a significance of 1.9σ. Their measurements converge towards the Planck prediction for more distant lenses (see their figure A1). This again strongly suggests that the Hubble tension is indeed driven by a local environmental effect.
Another technique to determine H 0 uses maser-derived distance and velocity measurements, as done by the Megamaser Cosmology Project (Reid et al. 2009). This method is independent of distance ladders, standard candles, and the CMB. It also faces rather different systematics to techniques that rely on gravitational lensing (Pesce et al. 2020). They used measurements for the six maser galaxies UGC 3789, NGC 6264, NGC 6323, NGC 5765b, CGCG 074-064, and NGC 4258. Except for the well-studied case of NGC 4258 (e.g. Reid et al. 2019), these galaxies are located at distances between 51.5 +4.5 −4.0 Mpc and 132.1 +21 −17 Mpc. The resulting H local 0 = 73.9 ± 3.0 km s −1 Mpc −1 , consistent with Wong et al. (2020) and again larger than predicted by Planck.
So far, we have distinguished between the Planck prediction and H 0 measurements from the local Universe that avoid assumptions about early Universe physics. Baryon acoustic oscillation (BAO) measurements combine the two through a CMB-based prior on the sound horizon at the time of last scattering. The co-moving length of this standard ruler is assumed to remain fixed, allowing its angular size at different epochs to constrain the expansion history (Eisenstein et al. 2005). Such BAO-based H 0 measurements are available from redshift surveys at effective redshifts of z eff = 0.38, 0.51, and 0.61 (Alam et al. 2017), with the range recently extended to z eff ≈ 1.5 ). These yield a Hubble parameter consistent with the Planck prediction.
The combination of clustering and weak lens-ing data, BAO, and light element abundances gives 67.4 +1.1 −1.2 km s −1 Mpc −1 (Dark Energy Survey & South Pole Telescope Collaborations 2018). Estimating H 0 using cosmic chronometers yields a nearly direct measure of the background cosmology. This is also consistent with Planck (Gómez-Valent & Amendola 2018). Assuming spatial flatness of the Universe, Ruan et al. (2019) combined cosmic chronometers with information on HII galaxies to show that the true value of a is much closer to the Planck value than the local value of Riess et al. (2016), with the latter discrepant at ≈ 3σ. Migkas et al. (2020) inferred H 0 from the X-ray luminosity-temperature relation of galaxy clusters, finding that it ranges from 65.20 ± 1.48 to 76.64 ± 1.41 km s −1 Mpc −1 for different sky regions (see their figure 23). This range is similar to that between H global 0 (Planck Collaboration VI 2020) and H local 0 as found using SNe (e.g. Riess et al. 2016Riess et al. , 2019Camarena & Marra 2020a) or strong lensing systems (Wong et al. 2020). The apparent anisotropy of the local velocity field could potentially be caused by our off-centre location within the KBC void, a non spherical void shape, or a combination of both. However, these considerations are beyond the scope of this work.
Remarkably, all these studies reveal that only the lowredshift probes prefer a high value for the Hubble constant, with high-redshift probes yielding similar results to the Planck-based prediction (see e.g. figure 12 in Wong et al. 2020, or figure 1 in Verde et al. 2019). Some recent reviews on the Hubble tension can be found in Verde et al. (2019) and Riess (2019). All these results point to the overall picture that the Hubble tension is driven by a local environmental effect like a void. In particular, the KBC void shows up not only in galaxy counts but also in the velocity field as an unexpected first and second time derivative of the apparent scale factor (as evidenced by the reported anomalies in H 0 and q 0 , respectively). As discussed in Section 1.1, a large local underdensity can potentially resolve the Hubble tension if mass conservation is assumed. Therefore, this would be a natural resolution to the Hubble tension that would minimize adjustments to the ΛCDM model on cosmological scales. In particular, there would be no need to assume a novel expansion rate history driven by yet more undetected sources such as early dark energy (e.g. Karwal & Kamionkowski 2016;Alexander & McDonough 2019;Poulin et al. 2019;Sakstein & Trodden 2020). 1 Instead, the standard ΛCDM expansion rate history could be preserved. In Section 5.3, we discuss some of the objections to this approach.
The works of Enea (2018) and Shanks et al. (2019) constitute attempts to relate the Hubble tension and KBC void on the basis of mass conservation. In a next step, one has to perform more sophisticated dynamical modelling with reasonable initial conditions provided by the CMB. As we will argue, this is not possible with the standard governing equations of ΛCDM (Section 2.2). In particular, Macpherson et al. (2018) explicitly showed that cosmic variance caused by inhomogeneities of the underlying density field cannot resolve the Hubble tension. This is because the expected 1 The work of Hill et al. (2020) argues that early dark energy cannot resolve the Hubble tension due to constraints from other data.
cosmic variance is too low, implying the Hubble tension and KBC void must both be measurement errors. Given the very different ways in which they are measured, this is highly implausible.
Thus, a large void and high H local 0 could well point to a different theory where both are explained by enhancing the long-range strength of gravity, which would promote the growth of structure. In principle, any alternative cosmological model that enhances cosmic variance through faster structure formation could explain the KBC void and Hubble tension, insofar as the model faces the Hubble tension. However, it is important for the model to explain phenomena in addition to those for which the model was explicitly designed, and to address observations on galaxy scales. Therefore, we concentrate on detailed dynamical modelling in the framework of an approach known to satisfy galaxyscale constraints, and to promote the growth of structure on larger scales.

Milgromian dynamics
Milgrom (1983) originally developed Milgromian dynamics (MOND) to explain the flattening of galactic rotation curves without the need of massive CDM haloes. MOND is a classical potential theory of gravity with a Lagrangian formalism (Bekenstein & Milgrom 1984). It explains the dynamical effects usually attributed to CDM by an accelerationdependent modification to Newtonian gravity. In particular, the gravity at radius r from an isolated point mass M becomes where G is the Newtonian gravitational constant, and a 0 is Milgrom's constant. Empirically, a 0 = 1.2 × 10 −10 m s −2 to match galaxy rotation curves (e.g. Begeman et al. 1991;McGaugh 2011). For a more complicated mass distribution, g follows a non-relativistic field equation (Bekenstein & Milgrom 1984). We use a more computer-friendly version known as quasilinear MOND (QUMOND; Milgrom 2010). In this approach, where Φ is the gravitational potential, g N is the Newtonian gravitational field, and r ≡ |r | for any vector r. The function ν g N interpolates between the Newtonian (|∇Φ| a 0 ) and deep-MOND (|∇Φ| a 0 ) regimes. Throughout this project, we apply the widely used 'simple' interpolating function (Famaey & Binney 2005): This closely approximates the empirically determined radial acceleration relation (RAR) between g N obtained from photometry and g ≡ −∇Φ obtained from rotation curves (McGaugh 2016;Lelli et al. 2017 where M b is the baryonic mass, v f is the asymptotic rotation velocity of a disc galaxy, and the exponent ξ = 4. Empirically, a tight relation of this form is evident with ξ ≈ 3 − 4 (e.g. McGaugh et al. 2000;McGaugh 2005;Stark et al. 2009;McGaugh 2011;Torres-Flores et al. 2011;Ponomareva et al. 2018). The more recent investigations put ξ very close to the MOND-predicted value of 4, which is also what we expect empirically based on the RAR. Since v f can be measured independently of distance but M b depends on the adopted distance, the BTFR provides another independent method to obtain H local 0 . Recently, Schombert et al. (2020) calibrated the BTFR with redshift-independent distance measurements from Cepheids and/or the tip magnitude of the red giant branch for 30 galaxies in the Spitzer Photometry and Accurate Rotation Curves catalogue (SPARC; Lelli et al. 2016) and 20 galaxies from Ponomareva et al. (2018). The so-calibrated BTFR was then applied to 95 independent SPARC galaxies for which only the redshift is known. Since the SPARC catalogue contains galaxies up to distances of ≈ 130 Mpc, Schombert et al. (2020) derived H 0 of the very local Universe. They got H local 0 = 75.1 ± 2.3(stat) ± 1.5(sys) km s −1 Mpc −1 (see also their table 5). This is quite consistent with other measurements from the late Universe and significantly exceeds the ΛCDM prediction based on the CMB (Section 1.2). Interestingly, the dominant source of systematic uncertainty is how to correct redshifts of SPARC galaxies for peculiar velocities induced by large-scale structure. This points towards mis-modelled peculiar velocities as a possible cause for the entire Hubble tension.
According to Equation 7, MOND is non-linear in the acceleration, which yields the interesting concept of the external field effect (EFE; Milgrom 1986). In contrast to Newtonian gravity, the non-linearity of Milgrom's law causes the internal gravitational forces within a MONDian subsystem to be affected by the external gravitational field from its environment even without any tides. This breaks the strong equivalence principle. The EFE has likely been observed in the declining rotation curves of some disc galaxies (Haghi et al. 2016) and the internal dynamics of dwarfs. For example, Crater II is a diffuse dwarf satellite galaxy of the MW at a distance of ≈ 120 kpc (Torrealba et al. 2016). Its observed velocity dispersion of 2.7±0.3 km s −1 (Caldwell et al. 2017) is below the isolated MOND prediction of 4 km s −1 (McGaugh 2016). Taking into account the Galactic EFE reduces the MOND prediction to 2.1 +0.9 −0.6 km s −1 , matching the observed value within uncertainties. Similar examples are the ultradiffuse dwarf galaxies Dragonfly 2 (DF2) and DF4, where the MOND predictions agree with observations only if the EFE is included Haghi et al. 2019a). For the more isolated galaxy DF44, the MOND prediction without the EFE is consistent with observations (Bílek et al. 2019;Haghi et al. 2019b).
The EFE is also important within the MW, whose MONDian escape velocity curve is similar to observations (Banik & Zhao 2018a). Since Equation 6 yields a logarith-mically divergent potential, escape from an isolated object is not possible in MOND unless the EFE is taken into account. Recently, Pittordis & Sutherland (2019) showed that MOND without an EFE is completely ruled out by the observed relative velocity distribution of wide binary stars in the Solar neighbourhood at separations of ≈ 10 kAU. Including the EFE leads to nearly Newtonian behaviour, though the predicted 20% difference is likely detectable in a more thorough analysis (Banik & Zhao 2018c) that must include contamination by undetected close companions (Clarke 2020).
In addition to its successes with internal dynamics of galaxies (reviewed in Famaey & McGaugh 2012), MOND may also explain the discs of satellites around the MW and M31 as TDGs born out of a past MW-M31 flyby. A previous close interaction is required in MOND (Zhao et al. 2013) due to the almost radial MW-M31 orbit (van der Marel et al. 2012(van der Marel et al. , 2019. In such an interaction, structures resembling satellite planes can be formed (Bílek et al. 2018). Using restricted N-body models to explore a wide range of flyby geometries, Banik et al. (2018) identified models where the tidal debris around the MW and M31 align with their observed satellite planes and have a similar radial extent. A past MW-M31 interaction would naturally explain the apparent correlation between their satellite planes, and with other structures in the LG (Pawlowski & McGaugh 2014). It may also account for the anomalous kinematics of the NGC 3109 association, which is difficult to understand in ΛCDM (Peebles 2017;Banik & Zhao 2018b).
Interestingly, there is an order of magnitude coincidence between the value of a 0 and the cosmic acceleration rate: where c is the speed of light (Milgrom 1983). This may indicate that MOND is related to a fundamental theory of quantum gravity (e.g . Milgrom 1999;Pazy 2013;Smolin 2017;Verlinde 2017). A bigger clue would come from tighter empirical constraints on the time evolution of a 0 , which at present are still weak (Milgrom 2017). Even so, his work showed that current data are sufficient to rule out the a −3/2 scaling required by the model of Zhao (2008), which additionally would have a very significant impact on the CMB (Sections 3.1.3 and 5.2.3).
Another intriguing coincidence is that the total matter density is very nearly 2π times the baryonic density, i.e. Ω m ≈ 2πΩ b (Milgrom 2020a). This could imply that the effective gravitational constant in a MONDian Friedmann equation is a factor of 2π larger than for a system decoupled from the cosmic expansion. However, we will not follow this interpretation here.
The first relativistic version of MOND was developed by Bekenstein (2004). This was modified slightly by Skordis & Z lośnik (2019) so that gravitational waves propagate at the speed of light, as required for consistency with the near-simultaneous detection of gravitational waves and their electromagnetic counterpart (Virgo & LIGO Collaborations 2017). The theory of Skordis & Z lośnik (2019) allows solutions where the background cosmology follows the standard Friedmann equations to high precision (see their section 4). We discuss this further in Section 3.1, where we explain why the expansion rate history and the power spectrum of the CMB should be nearly the same as in ΛCDM. Thus, MOND would suffer from the Hubble tension in just the same way as ΛCDM if H local 0 = a at the sub-per cent level. Fortunately, this might not be the case − Sanders (1998) showed that due to the long-range modification to gravity, MOND produces much larger and deeper voids than predicted by ΛCDM cosmology. Thus, MOND could be a promising framework to explain both the KBC void and the Hubble tension. We therefore extrapolate Milgrom's law of gravity from sub-kpc to Gpc scales. For the first time, we study the Hubble tension and KBC void in the context of MOND. We emphasize that MOND was originally designed to address discrepancies on galactic scales (Milgrom 1983), so no new assumptions are made specifically to address the latest data on the low-z distance-redshift relation and galaxy counts − apart from the usual assumption that the background follows a standard evolution to high precision (Section 3.1.1), and that MOND applies only to density deviations from the cosmic mean (e.g. Llinares et al. 2008;Angus & Diaferio 2011;Angus et al. 2013;Katz et al. 2013;Candlish 2016). In this context, we aim to provide a unified explanation for both the dynamical discrepancies on galaxy scales and the z < ∼ 0.2 matter density and velocity field given current constraints from the CMB.
The layout of this paper is as follows: In Section 2, we quantify the likelihood of the observed KBC void and how it might relate to the Hubble tension in a ΛCDM context. After introducing a cosmological MOND model in Section 3, we compare it to observations of the local Universe (Section 4). The implications for ΛCDM and MOND cosmologies are discussed in Section 5. We finally conclude in Section 6. Throughout this paper, co-moving distances are marked with the prefix 'c' (e.g. cMpc, cGpc).

ΛCDM FRAMEWORK
In this section, we describe how we use a cosmological ΛCDM simulation to quantify cosmic variance and thereby determine the likelihood of finding ourselves inside the observed KBC void in standard cosmology. We also consider the implications of our results when combined with the Hubble tension.

Cosmic variance in the Millennium XXL simulation
Millennium XXL (MXXL; Angulo et al. 2012) is a standard ΛCDM cosmological simulation that evolves 6720 3 DM particles from z = 63 forwards to z = 0. Though it only considers DM, baryonic physics should have a negligible role on the 300 Mpc scale we consider. The simulation box has a length of 3h −1 cGpc, resulting in a volume that is 216× larger than that of the Millennium simulation (Springel et al. 2005). The mass of a particle is 8.456 × 10 9 M and its Plummer-equivalent softening length is 13.7 kpc. The MXXL simulation assumes a flat ΛCDM cosmology consistent with WMAP-7 results, i.e. the present matter density parameter is Ω m,0 = 0.25, that of dark energy is Ω Λ,0 = 0.75, σ 8 = 0.9, H 0 = 73 km s −1 Mpc −1 , and the power spectrum is assumed to be of the Harrison-Zeldovich form (n s = 1). The baryonic mass of each subhalo is obtained by applying the semianalytic galaxy formation code l-galaxies (Springel et al. 2005) to the MXXL data (see also section 2.2 in Angulo et al. 2014).
We use MXXL to calculate the relative density contrast given by the stellar mass distribution in subhaloes with stellar mass M * > 10 10 h −1 M at z = 0. For this purpose, we consider 10 6 vantage points distributed on a Cartesian grid with a spacing of 30h −1 Mpc in each direction. To maximize the accuracy of our results, we use the nearest subhalo as our final choice for the vantage point. Our adopted minimum mass avoids an excessive computational cost, but still leaves enough subhaloes to accurately determine the expected cosmic variance. Using only stellar masses makes our results more comparable to observations in the NIR.
We need to allow for the incomplete sky coverage of Keenan et al. (2013). Following their section 2.5, we adopt a sky area of 37 080 deg 2 , which in dimensionless units is We assume the incompleteness is caused by observational difficulties at low Galactic latitudes. Thus, we define a mock Galactic spin axis by randomly generating a unit vector n i drawn from an isotropic distribution. We can then define an angle θ j based on the direction towards another subhalo at position r j relative to our vantage point.
The subscript i refers to the vantage point, while j refers to another subhalo observed from there. We mimic incomplete sky coverage by requiring that Since most of the sky is surveyed, cos θ obs = 0.10. The observed density contrast is calculated for galaxies in the redshift range 0.01 < z < 0.07 (table 1 in Keenan et al. 2013). Therefore, we further require selected subhaloes to satisfy where r min = 40 Mpc and r max = 300 Mpc. The relative density contrast around vantage point i is then The sum is taken over all subhaloes with M * > 10 10 h −1 M that satisfy Equations 13 and 15. These conditions restrict us to a volume V. The cosmic mean density ρ 0 is found by relaxing the position-related conditions and dividing the much larger sum by the whole simulation volume.

Comparison with observations
We now compare our so-obtained list of δ i with the observed local matter distribution. By combining our results with prior analytic work in ΛCDM, we also assess the implications for the Hubble tension and conduct a joint analysis. ) of spheres with a 300 Mpc radius less an inner 40 Mpc hole in the ΛCDM MXXL simulation, calculated at redshift z = 0 (Section 2.1). The red solid curve shows the observed density contrast of δ obs = 0.46±0.06 with Gaussian errors (see also figure 11 and table 1 in Keenan et al. 2013). Theδ values closely follow a Gaussian distribution with a dispersion of σ ΛCDM = 0.048 (black curve). A more detailed Gaussianity test is performed in Appendix A. Both curves are normalised to the same area.

KBC void
As discussed in Section 1.1, Keenan et al. (2013) discovered a large local underdensity with an apparent density contrast of δ obs = 0.46 ± 0.06 around the LG assuming a fixed distance-redshift relation with H 0 = 70 km s −1 Mpc −1 (see their section 4.7). To compare their reported δ obs with ΛCDM expectations, we need to account for the fact that any underdensity δ would also affect the local Hubble parameter by where e.g. Marra et al. (2013) showed that for δ 1 in ΛCDM, with the bias factor b = 1 (see also Section 5.3.1). As a result, the volume within a fixed redshift would be reduced below that assumed in Keenan et al. (2013) by a fraction The apparent underdensityδ i uncorrected for RSD would then be For the small underdensities expected in ΛCDM (see below), this approximately implies In other words, the apparent (RSD-uncorrected) underdensity would be 1.5× larger than the actual value. Figure 1 shows the distribution ofδ i in the standard ΛCDM MXXL simulation. This yields true rms density fluctuations of 3.2%, so observations uncorrected for RSD should exhibit fluctuations of 4.8%. To a very good approximation, these should be normally distributed, as demonstrated in Appendix A. Since 46/ √ 6 2 + 4.8 2 ≈ 6.0, we expect the discrepancy to be at the ≈ 6σ level.
Comparing the density contrast predicted by standard cosmology with the observed KBC void reveals a very significant discrepancy (Figure 1). This is usually quantified by finding the likelihood P of observing a more severe discrepancy, which we find for each vantage point and then average: Here, N = 10 6 is the number of vantage points, δ obs = 0.46 is the observed underdensity, and σ obs = 0.06 is its uncertainty. The function f χ →P gives the likelihood that a 1D Gaussian is more than χ standard deviations away from its mean. We use the inverse function f P →χ to convert the so-obtained Pvalue into a more easily understood form, as will usually be done throughout this article. In this way, we find that the KBC void is in 6.04σ tension with ΛCDM cosmology if it is accurately represented by the MXXL simulation on a 300 Mpc scale.

Implications for the Hubble tension
In any matter-conserving cosmological model, we expect an underdensity to be associated with some change in the local expansion rate (Equation 5). Figure 2 illustrates the manner in which this occurs for ΛCDM. In principle, the KBC void can boost the global Hubble constant to its local value observed by the SH0ES and H0LiCOW teams (H local 0 = 73.8 ± 1.1 km s −1 Mpc −1 , Riess et al. 2019;Wong et al. 2020). In fact, the straight line drawn on Figure 2 should curve to the right for large δ because as δ → 1, we expect that H local 0 /H global 0 → ∞ due to mass conservation (Equation 5, see also figure 1 of Marra et al. 2013). Thus, the expected relation between H local 0 andδ would pass rather close to the observations (red point). However, a 10σ density fluctuation would be necessary to reduce the Hubble tension to the 2σ level. Moreover, even a 5σ underdensity in ΛCDM is still not enough to get within 5σ of the local observations. This suggests that combining the KBC void and Hubble tension leads to a discrepancy with ΛCDM that slightly exceeds 5 √ 2σ = 7.07σ. We next perform a more detailed joint analysis.

Combined implications for ΛCDM
As discussed in Section 1.1, the locally measured H 0 is discrepant at the 5.3σ level with the Planck-based ΛCDM prediction (Wong et al. 2020) if we neglect the small expected impact of cosmic variance (Wojtak et al. 2014). In the previous section, we have shown that the KBC void is in 6.04σ tension with ΛCDM ( Figure 2). Therefore, both the KBC void and Hubble tension are difficult to explain within the = 67.4 ± 0.5 km s −1 Mpc −1 (Planck Collaboration VI 2020) and a local density equal to the cosmic mean (δ = 0). The red data point is the local Hubble constant combined from the SH0ES and H0LiCOW projects (H local 0 = 73.8 ± 1.1 km s −1 Mpc −1 ; Riess et al. 2019;Wong et al. 2020) and the locally observed δ obs = 0.46 ± 0.06 (Keenan et al. 2013). The grey contour lines show the indicated confidence levels assuming the measurements are independent. The blue points show the expected cosmic variance in ΛCDM corrected for RSD (Equations 18 and 22) at the indicated confidence level. Notice that a 5σ fluctuation is not enough to get within 5σ of the local observations. ΛCDM framework − we can explain both simultaneously, but this would require a 10σ density fluctuation ( Figure 2). In this context, the most plausible explanation is that both are caused by measurement errors. If so, we would have to assume two independent > 5σ errors, an unlikely scenario. The combined tension would correspond to χ 2 = 5.30 2 +6.04 2 for 2 degrees of freedom. This results in a probability of P = exp − χ 2 /2 = 9.4 × 10 −15 , which is equivalent to 7.75σ for one variable.
Measurements of the local density and velocity fields rely on rather different techniques, justifying our assumption of independence. For instance, a miscalibration of SNe magnitudes would affect H local 0 but not δ obs as the latter is a relative density contrast between different redshift bins. Thus, it is extremely unlikely that both phenomena are caused purely by measurement errors. Moreover, the KBC void is evident at different wavelengths as well as independently on smaller (< 50 Mpc) scales (Karachentsev 2012), while several independent teams have measured a higher local expansion rate than the Planck-based ΛCDM prediction (Sections 1.1 and 1.2, respectively).
A more rigorous way to estimate the combined tension is to average the P-values across different vantage points considering their individual δ i , how this would perturb the local expansion rate, and how the resulting RSD would lead to an enhanced apparentδ i . The average P-value is thus is the apparent local Hubble constant. Here, H local 0 = 73.8 km s −1 Mpc −1 and σ H 0 = 1.2 km s −1 Mpc −1 , with the latter including an allowance for the 0.5 km s −1 Mpc −1 uncertainty from Planck Collaboration VI (2020). This procedure reveals that the KBC void and Hubble tension falsify the ΛCDM framework at 7.09σ, in agreement with our earlier estimate. Our calculation of the cosmic variance in ΛCDM is derived from the stellar masses of subhaloes with M * > 10 10 h −1 M , which should be more than sufficient to accurately trace the matter distribution on a 300 Mpc scale. Moreover, our results are consistent with expectations from the Harrison-Zeldovich spectrum (Harrison 1970;Zeldovich 1972) and its early Universe normalisation required to match the CMB (Planck Collaboration VI 2020). In a ΛCDM context, this is parametrized using σ 8 , which implies rms fluctuations of 3.2% on a 300 Mpc scale at the present epoch. This agrees with our much more rigorous estimate using MXXL (Section 2.2.1).
Therefore, the KBC void is not a consequence of random measurement errors or density fluctuations expected in standard cosmology. Structure formation mainly depends on the underlying gravitational law, strongly suggesting that the observed KBC void cannot be explained by treating baryonic physics differently on galaxy scales.
Although cosmic variance in a standard context is insufficient to explain the KBC void and H 0 from low-redshift probes (e.g. Macpherson et al. 2018), Figure 2 indicates that a large local void appears to be a promising explanation for these local observations. Consequently, we next consider a long-range modification to gravity which should enhance cosmic variance while accurately explaining observations on galactic scales with a fixed acceleration threshold (Famaey & McGaugh 2012). Section 5.3 discusses some commonly used arguments for why the KBC void cannot solve the Hubble tension.

MOND FRAMEWORK
As shown in the previous section, the cosmic variance expected within the ΛCDM framework is insufficient to explain the KBC void and Hubble tension. Thus, we aim to investigate structure formation and the velocity field in MOND (Milgrom 1983). In this section, we first introduce a conservative MOND cosmology that has the same expansion rate history and overall matter content as ΛCDM, but with CDM replaced by hot dark matter (HDM) to account for light element abundances, galaxy clusters, and the CMB without much affecting galaxies (Angus 2009). We then explain how we parametrize the initial void density profile and evolve it forwards to the present time (Section 3.2). Finally, we describe how predictions for local observables are extracted from our models (Section 3.3).

The νHDM cosmological model
Any viable cosmological model has to explain the angular power spectrum of the CMB and the primordial abundances of light elements. Angus (2009) provided a promising cosmological model that seeks to address the shortcomings of MOND on galaxy cluster and larger scales using an extra sterile neutrino species with a mass of m ν s = 11 eV/c 2 . Thermally produced neutrinos of this mass would have the same relic abundance as CDM particles in standard cosmology, but would behave as HDM in the sense of not clustering on galaxy scales. 2 The composition of the universe as a whole would be similar to ΛCDM − baryons would still comprise ≈ 5% of the present critical density of the universe, sterile neutrinos would replace the ≈ 25% contribution of CDM, and dark energy would yield the remaining ≈ 70% (i.e. Ω m,0 = Ω b,0 + Ω ν s ,0 ≈ 0.3 and Ω Λ,0 ≈ 0.7). We refer to this model as the νHDM paradigm, where ν stands for both the interpolating function in QUMOND (Equation 8) and sterile neutrinos, maximizing the chance that it is physically meaningful. The observed expansion history of the Universe seems broadly consistent with ΛCDM cosmology (e.g. Joudaki et al. 2018). As shown by Angus (2009), νHDM yields the same expansion history as ΛCDM due to the same overall matter content and the same Friedmann equations at the background level ). This issue is discussed further in Section 3.1.1.
Although the existence of sterile neutrinos is not experimentally confirmed yet, they are theoretically consistent with standard particle physics (Merle 2017). Observationally, the νHDM model is motivated mainly by galaxy clusters, where the dynamical discrepancy cannot be explained in MOND without DM (Sanders 2003). Furthermore, DM is necessary to address the offset between X-ray and lensing peaks in the Bullet Cluster (Clowe et al. 2006), since MOND acting on the baryons alone is unable to fully replace the role played by CDM in standard cosmology (Angus et al. 2007). We emphasize that these observations do not uniquely require CDM since they are on a much larger spatial scale than the hypothesized CDM haloes of individual galaxies (Ostriker & Peebles 1973).
In this context, Angus et al. (2010) analysed 30 of the most virialized galaxy groups and clusters in the νHDM paradigm. They found that the required HDM density in all cases reaches the so-called Tremaine-Gunn limit (Tremaine & Gunn 1979) at the centre for sterile neutrinos with m ν s = 11 eV/c 2 . This is a strong indication that the DM density in galaxy cluster cores is limited by quantum degeneracy pressure (the Pauli Exclusion Principle). Note that MOND fits to galaxy rotation curves are hardly affected by sterile neutrinos with m ν s < ∼ 100 eV/c 2 , even if their number density reaches the Tremaine-Gunn limit (section 4.4 of Angus et al. 2010). As a result, νHDM is likely to explain the internal dynamics of both galaxies and galaxy clusters. Introducing sterile neutrinos is thus well consistent with astronomical observations and almost consistent with the standard model of particle physics (unlike CDM particles), but they nevertheless require experimental verification.
In the following, we address the background evolution of a (t) in the νHDM framework, allowing us to address the primordial abundances of light elements and the CMB. We also consider the implications for large-scale structure, where substantial differences are expected from ΛCDM. The theoretical uncertainties of the here applied MOND approach are summarized in Section 5.2.3, which focuses on how density perturbations should be treated in MOND.

Background cosmology
The background evolution a (t) requires a relativistic theory that yields the appropriate MOND limit in galaxies. In this contribution, we make certain assumptions about the parent relativistic theory that gives rise to MOND. These assumptions are based on prior work, in particular with the tensorvector-scalar (TeVeS) theory that was the first covariant framework with an appropriate MOND limit (Bekenstein 2004). His section 7 indicates that the background evolution should be very similar to General Relativity at all epochs for the same matter-energy content.
The background evolution and perturbations in TeVeS were addressed in detailed calculations done by Skordis (2006). To avoid detectable departures from the standard expansion history during the nucleosynthesis era, the free dimensionless parameter µ 0 must be rather large . 3 In particular, if we allow the extra energy density contributed by the scalar field to comprise a fraction X of the critical density during the radiation-dominated era, then the contribution in the matter and Λ-dominated eras would be X/9. Primordial light element abundances then imply that the standard Friedmann equation would differ from the TeVeS cosmology at only the sub-per cent level (see their figure 1). The very small contribution of the scalar field density was also demonstrated in figure 2 of Dodelson & Liguori (2006). Therefore, we will assume that the background cosmology is identical to that of ΛCDM. Since the CMB is also expected to have similar properties in both frameworks (Section 3.1.3), they both lead to the Hubble tension in a similar manner provided that a = H local 0 , i.e. if cosmic variance in the local measurements is much smaller than the Hubble tension. Our main argument is that this assumption is valid in ΛCDM but need not be in MOND.
While the original version of TeVeS is inconsistent with gravitational waves travelling at c, a slightly modified version does have this property, even in the presence of perturbations (Skordis & Z lośnik 2019). The above-mentioned results should carry over to the updated version of TeVeS, though this should be carefully demonstrated in future work. The preliminary results of Skordis & Z losnik (2020) are an important step in this direction.
Throughout this work, we assume dark energy not to be an artefact of an observer in an underdense region seeing an apparently accelerating expansion due to the developing inhomogeneities (Buchert 2000). However, we emphasize that proper time-averaging of global properties of the universe would be required to further study the present model (Wiltshire 2007).

Big Bang nucleosynthesis
Big Bang nucleosynthesis (BBN) occurred at a temperature of kT ≈ 1 MeV, where k is the Boltzmann constant. A review on BBN can be found e.g. in Cyburt et al. (2016). In the νHDM framework, Skordis et al. (2006) showed that it is possible to have essentially no departure from the standard expansion history during the radiation-dominated era. However, the model would still have an effect on BBN because at kT ≈ 1 MeV, sterile neutrinos with m ν s ≈ 11 eV/c 2 would be relativistic. Their weaker interactions would cause them to decouple earlier, so they would add an extra 7/8 to g * , the number of effective relativistic degrees of freedom. Since the Hubble parameter H ≡ a/a scales as H ∝ √ g * and standard physics predicts g * = 10.75, this would increase H by only 4%, causing a slight impact on the primordial abundances of light elements. As shown in equation 13 of Cyburt et al. (2016), any increase in H raises the primordial He-4 mass fraction Y p because free neutrons have less time to decay. Their detailed calculations have shown that this dependence can be fitted with a power law of the form In standard cosmology, the effective neutrino number is N ν = 3.046, which slightly exceeds 3 because neutrinos decouple only slightly before electron-positron annihilation at kT = 511 keV. Thus, an extra sterile neutrino species would increase Y p by a factor of (4.046/3.046) 0.163 = 1.047, implying the standard value of Y p = 0.247 would rise to 0.259. This is only a small effect, so observations of the primordial He abundance in ancient gas clouds currently do not set a strong constraint on the existence of an extra sterile neutrino. For instance, measurements of the He abundance of a gas cloud at z = 1.724 backlit by a quasar yield Y = 0.250 +0.033 −0.025 (Cooke & Fumagalli 2018). Using a sample of H ii regions, Aver et al.
Measurements of the primordial abundances of D and Li-7 are less sensitive to N ν (Cyburt et al. 2002). However, primordial D abundances are relatively well known.  obtained N ν = 3.41 ± 0.45 based on (D/H) p derived from a metal-poor damped Ly α system. Therefore, both D and He measurements allow an extra sterile neutrino, which was actually favoured by the earlier analysis of Steigman (2012). We do not consider the more problematic case of Li-7, though see Howk et al. (2012) for a gas phase measurement in the Small Magellanic Cloud that seems to resolve the lithium problem.
These considerations only hold for sterile neutrinos in thermal equilibrium during the nucleosynthesis era. However, if sterile neutrinos decoupled much earlier, their number density could be lower depending on whether any other particle subsequently became non-relativistic. If so, ∆N ν would be lower, reducing the impact on g * and on BBN.
This scenario would require a higher sterile neutrino mass to recover the standard value of Ω m .

Radiation-dominated era and the CMB
After BBN, the next major constraint on any cosmological model comes from the CMB. This occurred shortly after the epoch of matter-radiation equality at z eq = 3411±48 (Planck Collaboration VI 2020). 4 This corresponds to a photon temperature of kT ≈ 0.80 eV, which is much less than the mass of the here considered sterile neutrinos. Consequently, they would behave just like non-relativistic CDM, causing z eq to be the same as in the ΛCDM model.
The CMB was emitted at z CMB ≈ 1100, corresponding to kT ≈ 0.26 eV. At this time, matter dominated the energy budget of the universe. Since the background cosmology of the νHDM model is the same as for ΛCDM and the plasma physics is unchanged, the sound horizon at recombination would still have the standard value of 147.09 ± 0.26 cMpc (Planck Collaboration VI 2020). This is directly related to the angular scale of the first acoustic peak in the CMB, which should thus be unaffected in our model. 11 eV/c 2 sterile neutrinos would be non-relativistic at the time of last scattering. Since both T and the peculiar velocity v pec should decline ∝ 1/a, we expect the sterile neutrinos to typically have This implies a free-streaming length of L fs ≈ 3.5 cMpc, which is much shorter than the horizon scale. Since the first acoustic peak of the CMB occurs at a multipole moment of ≈ 200 (Jaffe et al. 2001), free-streaming becomes important only for > ∼ 200/ 0.024 √ 3 = 4900, beyond the range accessible by Planck Collaboration VI (2020). This is consistent with section 6.4.3 of Planck Collaboration XIII (2016), which explicitly states that any particles with m > 10 eV/c 2 "are so massive that their effect on the CMB spectra is identical to that of CDM." The νHDM paradigm does more than simply replace CDM with HDM. Because of the Milgromian force law, the paradigms differ with regards to the evolution of sub-horizon perturbations. In the following, we estimate the gravitational field from inhomogeneities around t CMB , the time of recombination.
The peculiar velocities are of order v pec ≈ cδ and were built up over a duration of t CMB = 380 kyr. Assuming rms density fluctuations of δ CMB = 10 −5 as observed in the baryons, we can obtain a lower bound on the peculiar acceleration g CMB sourced by inhomogeneities.
This already exceeds Milgrom's constant a 0 . However, the gravity must have been significantly stronger to compensate for resistance from radiation pressure. In order to estimate the density fluctuations in the HDM component at t CMB , we consider the value of σ 8 = 0.811 ± 0.006 on a scale of 8h −1 ≈ 12 cMpc that is required to fit the CMB anisotropies (Planck Collaboration VI 2020). For a scale-invariant power spectrum, the density fluctuations on the 147 cMpc scale of the first acoustic peak in the CMB are 12σ 8 /147 ≈ 0.065 at the present epoch, as can also be seen by scaling our results of Section 2.2 for fluctuations on a 300 cMpc scale. 5 Since ΛCDM predicts that δ ∝ a in the matter-dominated era and neglecting the effect of dark energy, we would expect density fluctuations of δ CMB ≈ 5.9 × 10 −5 at t CMB . Taking into account that structure formation slowed down when the Universe became dark energy-dominated at z ≤ 0.7 and was slower around the time of recombination due to the still significant amount of radiation, we estimate that Thus, the typical gravitational field at recombination was implying that MOND had only a minor impact at that time.
In the matter-dominated era (a a eq ), the density perturbations grow ∝ a after their mode enters the horizon. Therefore, the Harrison-Zeldovich power spectrum predicts that the power of the density perturbations scales inversely with their length L (Harrison 1970;Zeldovich 1972) Since the mass enclosed by the mode is M ∝ L 3 , the mass perturbation must scale as Therefore, the perturbation's Newtonian gravity is independent of L, i.e.
The Harrison-Zeldovich power spectrum breaks down for length-scales that enter the horizon before a eq . Since no modes would be able to grow during the radiationdominated era, these short-wavelength modes would have much less power than predicted by a 1/L scaling relation. Thus, g N would be smaller. However, in MOND, these shortwavelength modes would be embedded in the EFE generated especially by long-range modes (Section 1.3). This would severely limit the MOND boost to the internal gravity of shorter modes, since their total g N depends on both their internal gravity and any external field. For this reason, we expect that modes of any L were unaffected by MOND around the epoch of recombination. We next consider how this picture changes with time. Since Newtonian density perturbations are expected to grow as δ ∝ a in the matter-dominated era, the mass perturbation should also scale as For linear (δ 1) perturbations whose co-moving size 5 The here used MXXL simulation is calibrated to the CMB data gathered by WMAP-7 (Angulo et al. 2012).
hardly changes, the Newtonian gravity should scale as Our previous estimation showed that the gravitational field sourced by inhomogeneities is g a 0 at t CMB (Equation 32). We now see that even larger gravitational fields are expected at earlier times, further justifying our assumption that MOND would have little effect then. 6 We can combine Equations 32 and 37 to deduce that MOND does not play a significant role in structure formation until z < ∼ z MOND = 50. This underpins the commonly used assumption that MOND does not play a role in the very early universe, but would promote the formation of the first galaxies (Sanders 1998).
The high accelerations around the time of recombination strongly suggest that the MOND gravity law would not by itself affect the acoustic oscillations in the CMB. This issue was investigated further by Skordis et al. (2006), who considered a covariant formulation of MOND. Their figure 2 confirms our conclusion that the modification to gravity has by itself only a very small effect for plausible choices of the model parameters consistent with BBN. However, their use of three ordinary neutrino species with a much lower mass of 2 eV/c 2 led to significant free streaming effects that are totally inconsistent with the latest observations (Planck Collaboration XXVII 2014). If instead a single 11 eV/c 2 sterile neutrino is used, a very good fit can be obtained to the CMB power spectrum for the reasons just discussed (figure 1 of Angus 2009). Note also that with a standard a (t), the angular diameter distance to the CMB would be the same as in ΛCDM, placing the acoustic peaks at the correct angular scales. Indeed, figure 1 of Angus & Diaferio (2011) shows that the CMB power spectra in the νHDM and ΛCDM models agree quite closely, so both paradigms are consistent with observations taken by WMAP-7, the Atacama Cosmology Telescope (ACT), and the Arcminute Cosmology Bolometer Array Receiver up to = 2500.

Evolution of perturbations and large-scale structure
Even if the CMB power spectrum is correct in our framework, the observed CMB is also influenced by foreground structures. Section 5.3.3 discusses the gravitational redshift of the entire last scattering surface due to the rather high MOND potential of the KBC void. Foreground lensing of the CMB by large scale structures and the integrated Sachs-Wolfe (ISW) effect would also be stronger in MOND. There are some observational hints that these effects are stronger than expected in ΛCDM (Section 5.3.1). These tensions could be eased in a theory where structure formation is more efficient. However, it is possible that νHDM overcorrects the problem and produces too much foreground lensing and/or a Sachs-Wolfe effect in disagreement with observations. These issues are beyond the scope of our work, but should be addressed before the νHDM framework can be considered to fully account for all observed aspects of the CMB. This would almost certainly require numerical simulations of structure formation. In addition, photon propagation through such a simulation would need to be handled with care, taking account of inhomogeneities and their time evolution (e.g. Wiltshire 2007). Nusser (2002) considered the growth of density perturbations in a Milgromian framework. Their section 2 introduced the basic principle used in all subsequent MOND cosmological simulations (Llinares et al. 2008;Angus & Diaferio 2011;Angus et al. 2013;Katz et al. 2013;Candlish 2016). These simulations make the ansatz that a MONDified Poisson equation (usually Equation 7) is applied only to the density perturbations about the mean background value, as evident e.g. in equation 2 of Candlish (2016). 7 This 'Jeans swindle' (Binney & Tremaine 1987) approach to MOND was justified using the earlier work of Sanders (2001), who showed its validity in a non-relativistic Lagrangian formulation of MOND (see his section 2). The approach is certainly valid for systems such as galaxies that are much denser than the cosmic mean. The use of non-relativistic gravitational equations should be sufficient when dealing with structures such as the KBC void that are much smaller than the cosmic horizon, since gravity travel time effects would not be too significant. Falco et al. (2013) showed that the Jeans swindle is formally correct in Newtonian gravity − including the background would simply add on the force required to maintain the time-dependent Hubble flow velocity. However, it still needs to be rigorously demonstrated that the swindle remains mathematically valid in a MONDian model with a non-linear gravity law. Therefore, although this ansatz is commonly used by the MOND community, it is one of the strongest assumptions in the here presented cosmological model.
One of the few works that does not make this assumption is Sanders (2001), whose model is a non-relativistic two-field Lagrangian-based theory of MOND. The coupling between these two fields is described by an adjustable parameter β in his modified Poisson equation 8. Setting β = 0 is equivalent to applying the Jeans swindle approach. However, if β 0, there exists a coupling between the peculiar acceleration sourced by inhomogeneities and the zeroth-order Hubble flow acceleration g Hubble (Equation 38). Sanders (2001) adopted β = 3.5 for his main analysis. As discussed in the cosmology section of Sanders & McGaugh (2002), g Hubble essentially contributes an extra source of gravity to the total entering the ν calculation in Equation 7, limiting the MOND boost to gravity. We call this the 'Hubble field effect' (HFE), since it is similar to but distinct from the usual EFE in MOND − both make the behaviour more Newtonian. In Section 5.2.3, we address theoretical uncertainties arising from the HFE, which is neglected in our main analysis. A non-zero HFE would substantially affect large-scale structures especially at scales > ∼ 100 cMpc, which could be used to constrain it in future studies (Section 5.2.3). However, we argue there that even with a strong HFE, cosmic variance would still be enhanced 3× compared to ΛCDM expectations on a 300 Mpc scale under conservative assumptions, enough to reproduce the KBC void.
Nusser (2002) built on the model of Sanders (2001) but assumed instead that β = 0 because he could not find any physical justification for coupling both fields, i.e. for the HFE. This uncoupled (Jeans swindle) approach is generally the one adopted in MOND cosmological simulations (e.g. Llinares et al. 2008;Angus et al. 2013;Katz et al. 2013;Candlish 2016). In particular, Angus et al. (2013) used it in a cosmological N-body simulation designed to address the formation of large-scale structure in MOND supplemented by sterile neutrinos. Although their work was novel and very advanced for its time, it faces some conceptual and numerical problems. In particular, they concluded that their model with 11 eV/c 2 sterile neutrinos significantly underestimates the number of low-mass galaxy clusters and slightly overestimates the number of very massive clusters (see e.g. their figure 4). This inconsistency between the model and observational data could arise for several reasons. Their conclusion is based on a simulation with a box size of 256h −1 cMpc and a particle resolution of only ≈ 3.78 × 10 10 M . The underproduction of low-mass galaxy clusters could be explained by the low particle resolution and therewith by an absence of low-mass particles needed to form such systems. In addition, they do not use a grid with adaptive mesh refinement (AMR), which causes that the potential wells especially of the smaller clusters may not be resolved properly, making them difficult to form. Therefore, it would be highly valuable to revisit their cosmological simulations with an AMR grid code such as phantom of ramses (Lüghausen et al. 2015), which adapts the potential solver of the widely used ramses algorithm (Teyssier 2002).
In general, small simulation boxes lack large-scale modes. Since the EFE is mainly sourced by very massive objects, a too small simulation box would potentially underestimate the EFE on MONDian subsystems. Thus, the internal gravitational field would be too strong, which could also explain the efficient formation of massive galaxy clusters in Angus et al. (2013).
As already discussed at the beginning of this section, Angus et al. (2010) demonstrated that the required neutrino density in 30 virialized galaxy groups and clusters reaches the Tremaine-Gunn limit at the centre, which supports the νHDM model. However, the neutrino degeneracy pressure in the cores of galaxy clusters has not been included in the simulations of Angus et al. (2013). If one would account for this effect, it would be more difficult to form massive galaxy clusters because gravity is resisted by neutrino degeneracy pressure.
Finally, Angus et al. (2013) compared their simulated halo mass functions with cluster mass functions derived from observations at z ≤ 0.3 (Reiprich & Böhringer 2002) and z ≤ 0.1 (Rines et al. 2008). As we have seen in Section 1.1, the KBC void has a similar extent. It is evident in X-ray galaxy cluster surveys (e.g. Böhringer et al. 2015Böhringer et al. , 2020. Therefore, local observations are biased against high-mass clusters, e.g. the massive merging galaxy cluster El Gordo (ACT-CL J0102-4915, Marriage et al. 2011) with a mass of 3 × 10 15 M (Jee et al. 2014) at z = 0.87 (Menanteau et al. 2012) would almost certainly not be evident in local observations from within a deep void. Thus, local observa-tions do not provide a representative cluster mass function of the whole Universe, so cannot be compared with the entire simulated halo population.
Consequently, the Angus et al. (2013) cosmological model has never been tested in full detail on large scales. An object similar to El Gordo was identified in the νHDM simulation of Katz et al. (2013), so initial results seem promising. It would be highly valuable to revisit their analysis in more physically and numerically advanced large-scale simulations. This is because the νHDM framework provides a viable explanation for BBN and the CMB, but also works on galaxy cluster scales while recovering the successes of MOND in galaxies. At present, there is no N-body or hydrodynamical simulation with a large enough box size to study the KBC void in a MONDian framework. Therefore, we develop a semi-analytic simulation for this purpose. In the following, we introduce the governing equations and parameters of the here discussed νHDM cosmological model.

Governing equations
We develop a simplified simulation in which the trajectories of particles are integrated up to the present time from z = 9, which corresponds to ≈ 0.5 Gyr after the Big Bang (Equation 45). As derived from General Relativity in section 2.2 of Banik & Zhao (2016), the particle's trajectory is described by the background cosmological acceleration term and any additional gravity sourced by inhomogeneities: where r is the particle's position relative to the void centre, g void is the local gravitational acceleration sourced only by density deviations from the cosmic mean, g Hubble is the acceleration in a homogeneously expanding spacetime, and i subscripts denote initial values when a = 0.1. At that time, particles are assumed to be on the Hubble flow. However, the initial matter distribution is assumed to be inhomogeneous. A spherically symmetric underdensity causes a Newtonian gravitational force of where ∆M is the mass deficit within radius r, ρ 0 is the present cosmic mean density of matter, and M enc is the enclosed mass. Since we assume mass conservation and no shell crossing, M enc remains constant for an individual particle. In the case of no void, g N = 0 since ∆M = 0. The exact set-up of the initial void profile is described in Section 3.2.1 and Appendix B.
Applying the Jeans swindle approach to MOND (Section 3.1.4), the gravitational force g is calculated with the 'simple' interpolation function (Equation 8) between the Newtonian and deep-MOND regimes (Famaey & Binney 2005). The EFE is included by quadrature summing g N and the Newtonian-equivalent external field g N,ext ): The EFE and its impact on the void will be described in more detail in Sections 3.2.2 and 3.3.6, respectively. Milgrom's constant a 0 = 1.2×10 −10 m s −2 is taken to be constant over cosmic time. Substantially higher values in the past may conflict with the CMB (Section 3.1.3) and high-redshift rotation curves (Milgrom 2017). Solving Equation 38 requires knowledge of the background cosmology. As argued in Section 3.1.1, assuming this follows a standard Friedmann equation should be accurate at the sub-per cent level. We therefore apply the second Friedmann equation and assume a standard flat background cosmology (Ω m,0 + Ω Λ,0 = 1), yielding where ρ m and ρ Λ are the cosmic mean densities of matter and dark energy, respectively. We assume that ρ m ∝ a −3 while ρ Λ = const. The parameters Ω m,0 and Ω Λ,0 are the present-day matter and the dark energy densities in units of the critical density ρ c = 3H 2 0 /(8πG). We set Ω m,0 = 0.315, Ω Λ,0 = 0.685, and choose a global Hubble constant of H 0 = 67.4 km s −1 Mpc −1 , consistently with the latest Planck data (Planck Collaboration VI 2020). Imposing the boundary conditions a = 0 when t = 0 and a = H 0 at a = 1, we get that

Initial void profile
The implemented void in the fiducial simulation run is initialized with a Maxwell-Boltzmann radial density profile. This is motivated by the observed Local Volume, where the density increases inwards for distances < ∼ 40 Mpc (see e.g. figure 3 in Karachentsev & Telikova 2018). The enclosed mass within co-moving radius r com from the void centre is thus given by The dimensionless radius x ≡ r com /r void , while α void is the initial void strength and r void is the parameter determining its co-moving size at z = 9. The first term in Equation 46 is the mass within a sphere of co-moving radius r com if the density were equal to the cosmic mean, with the void arising from the mass deficit imposed by the second term.
We run different simulations with α void ranging from 10 −5 to 10 −2 and r void ranging from (50 − 1030) cMpc. The parameter range of the initial void strength is motivated by the expected density fluctuations at z = 9 based on CMB data. In addition, we also run simulations in which the void is modelled with a Gaussian or an exponential initial density profile (Appendices B and C).

External field history
As stated in Section 1.3, the EFE is a consequence of the non-linearity of Milgrom's law of gravity (Milgrom 1986). Thus, we allow for the possibility that the void as a whole is embedded in an EF from even larger scales. We follow the usual approach of assuming the EF is sourced by a distant point-like object. This allows us to obtain the present-day Newtonian-equivalent external field using the simple interpolation function (Famaey & Binney 2005): whereg ext is the external field in units of a 0 . The evolution of the EFE over cosmic time is unknown due to the lack of a fully self-consistent MONDian framework. Since the EFE depends on the environment in which the MONDian system is embedded and thus on the formation of structure, we assume that the external field has a power-law dependence on the cosmic scale factor: where t 0 = 13.8 Gyr is the present time, and n EFE is a free parameter ranging from −2 to +2 in steps of 0.5 for different models. For our fiducial simulation run, we adopt a timeindependent external field (n EFE = 0). The results for different external field histories are discussed in Section 5.2.2. Table 1 summarizes the fixed and free parameters of our models.

Extracting mock observables
Our cosmological MOND models are constrained by the observed density contrast of the KBC void (Keenan et al. 2013), the local Hubble constant and deceleration parameter derived jointly from SNe data (Camarena & Marra 2020a), the Hubble constant from strong lensing (Wong et al. 2020;Shajib et al. 2020), and the motion of the LG wrt. the CMB (Kogut et al. 1993). In the following, we explain how we obtain the corresponding simulated quantities. Our approach involves comparing the void models described in Section 3.2 with a control simulation of a void-free standard cosmology. The control trajectories have a fixed comoving radius: Since the lookback time can be derived from SNe luminosities or angular diameter distances in a standard background cosmology, we fix this variable between the void and control models, allowing us to analyse the difference in other variables. The main advantage of this approach is that in the absence of a local void, our calculated late-time cosmological parameters (e.g. H 0 and q 0 ) would revert to their values in standard cosmology. Local observations imply that we are located close to the void centre (Keenan et al. 2013;Karachentsev & Telikova 2018). Therefore, as a simplification we assume in our analysis that we are at the void centre (Sections 3.3.2−3.3.5), except when calculating the likelihood of the observed LG peculiar velocity (Section 3.3.6). It is beyond the scope of our work to analyse the Hubble diagram and density field that might be seen by a substantially off-centre observer.

Apparent scale factor
The main quantity we extract is the redshift experienced by a photon as it travels from a particle to the void centre. This is given by where λ obs and λ emit are the wavelengths of the light as measured by the observer and at the source of emission, respectively, v int is the peculiar velocity of the particle relative to the void centre, and g void is the gravity in the radially outwards direction. The factor of a −1 arises from expansion of the universe while light from the particle is travelling towards us. This is the only factor that needs to be considered even without the void. The term marked 'Doppler' is the special relativistic Doppler effect, while the exponential factor (marked 'GR') is the gravitational redshift that arises because photons must climb up the void potential to reach its centre. As discussed in Section 3.3.5, relativistic lensing in MOND should yield similar results to General Relativity for the same g.
To limit the complexity of our algorithm and because we are dealing with a void at low z, we approximate the GR contribution by assuming the final density profile of the void is also applicable at earlier times. This leads to a time-independent gravitational field g void (r). We use this to calculate the integral in Equation 52 out to the co-moving distance where our past lightcone intersects the particle's trajectory (Section 3.3.3).
Since the observed SNe and lensing Hubble diagrams reported by observers are not corrected for the large peculiar velocities we expect in our model, the apparent scale factor is simply We compare the behaviour of this a app with the corresponding values in our control simulations, which are governed by Equation 51. Since we run a finite number of trajectories for each model, we interpolate between them to ensure the comparison is done at fixed lookback time. Void parameters α void Initial void strength at z = 9 (10 −5 , 10 −2 ) r void Initial void size at z = 9 (50 cMpc, 1030 cMpc)

Density contrast and redshift space distortion
In our models, the fractional underdensity inside a shell between radii r min,now and r max,now at the present time is Here, r min,initial and r max,initial are the initial co-moving distances of particles which are currently at r min,now and r max,now , respectively, and x min ≡ r min,initial /r void . Similar procedures are used to calculate x max and I max . The first term represents the initial density contrast, while the second accounts for expansion of the co-moving volume enclosed by the two shells.
As discussed in Section 1.1, the analysis of Keenan et al. (2013) used a distance-redshift relation based on the assumption of no void (see their section 4.7). Therefore, we apply an RSD correction to the observed relative density contrast in order to estimate the true value: Here, δ obs is the observed relative density contrast between the distances r void,in and r void,out at the present time. However, observations uncorrected for RSD are reported as if the known redshift range of the survey covers the distance range r control,in − r control,out , which are the corresponding distances to the same z in a void-free universe. The number of galaxies counted by the observers thus corresponds to a different δ than what they report, which is the RSD effect. Note that its magnitude will depend on the void model, so it is not possible to know the true density contrast in a model-independent way. This is because it is not possible to convert redshifts to distances without a dynamical model of the void. As a result, the uncertainty σ obs,corr is also model-dependent.
We can compare the so-corrected observed δ to the model prediction (Equation 54). This leads to a χ 2 contribution of which is calculated for the relative density contrasts in the redshift range 0.01 < z < 0.07 and between distances of 600 Mpc and 800 Mpc at the present time. According to Keenan et al. (2013), we estimate that δ obs,in = 0.46 ± 0.06 in the inner part of the void, while δ obs,out = 0.0 ± 0.1 in its outer part (see their table 1 and figure 11).

Lightcone analysis
To determine exactly when we would observe a test particle, we need to determine the intersection between its trajectory and our past lightcone. This occurs when the co-moving distance travelled by a light ray emitted from a particle equals the time-dependent co-moving distance to the particle. In other words, where t LC is the cosmic time when our past lightcone intersects a particle's trajectory. This is obtained by solving Equation 59 using the Newton-Raphson algorithm. We can then calculate relevant quantities at that time, which is used in our analyses related to the Hubble diagram (Section 3.3.4). However, when comparing the simulated v pec with the observed LG peculiar velocity (Section 3.3.6), we need to extract v pec at the present epoch since the measurement relates to the LG motion today. To limit the complexity of our analysis, we also use the present positions of particles when determining the density field of the void (Section 3.3.2). This should be valid if the void has not appreciably changed in the time needed for light to cross it, which is reasonable for a void much smaller than the Hubble distance c/H 0 = 4.4 Gpc.

Hubble constant and deceleration parameter from SNe
We constrain our models with the results of Camarena & Marra (2020a), who derived the local Hubble constant and deceleration parameter jointly from Pantheon SNe in the redshift range 0.023 ≤ z ≤ 0.15. As discussed earlier, we first find the difference in the apparent scale factor between our void model and a control void-free model.
where a app is the apparent scale factor (Equation 53), and a control is the scale factor at the same cosmic time in the control model of a void-free standard cosmology. Expanding Equation 60 as a Taylor series in the vicinity of the present time t 0 , we get that Dividing the above equation by a t 0 ≡ 1 and using the definitions of the Hubble parameter (H ≡ a/a) and deceleration parameter (q ≡ −a a/ a 2 ), we obtain that where ∆H 0 and ∆ q 0 H 2 0 are the boosts to these parameters due to the void, and the 0 subscripts denote present-day quantities. We find these by fitting ∆a (t) using a parabola forced to pass through ∆a = 0 at t = t 0 . The local Hubble and deceleration parameters are thus Because of a historical accident where it was assumed that the expansion of the Universe should decelerate, q 0 was defined as the present deceleration parameter −a a/ a 2 . It was subsequently shown that the Universe accelerates, implying q 0 < 0 (Riess et al. 1998;Schmidt et al. 1998;The Supernova Cosmology Project 1999). In order to minimize confusion from unnecessary use of − signs, we introduce from now on the acceleration parameter: The ΛCDM theory with the parameters obtained from Planck Collaboration VI (2020) predicts q 0 = Ω Λ,0 − 1 2 Ω m,0 = 0.53. In the absence of a void, H local 0 and q local 0 become identical to the Planck values since we use those for the background cosmology (Equation 45). The combined χ 2 contribution from H 0 and q 0 is Observationally, q local 0 = 1.08 ± 0.29 and H local 0 = 75.35 ± 1.68 km s −1 Mpc −1 , with a mutual correlation coefficient of C = 0.515 (Camarena & Marra 2020a). Their section 4 mentions that their posterior inference is very close to Gaussian, justifying our χ 2 approach.

Hubble constant from strong lensing
Empirically, it has been shown that light deflection in strong lenses works similar to General Relativity for the same nonrelativistic g (Collett et al. 2018). In the H0LiCOW lenses, g is constrained using the positions and time delays between images and also with velocity dispersion data. Their analysis should remain valid even in a MOND context, since the results of Collett et al. (2018) can be reproduced in relativistic versions of MOND (Milgrom 2013). The latter work showed that this approach works well empirically even in the deep-MOND regime, which can only be probed using weak lensing. This is because strong lensing always occurs in the Newtonian regime due to MOND's cosmological coincidence (Equation 10), as explained in Sanders (1999). Hence, strong lensing is little affected by MOND. None the less, we discuss in Section 5.2.1 how H 0 measurements from strong lensing impact our analysis, and consider the effect of excluding these measurements. In our main analysis, we constrain our MOND models with H 0 measured from seven strong-lens systems with the deflector at redshift z d . Our data set is derived from Shajib et al. (2020) and Wong et al. (2020). In our models, the Hubble constant at redshift z d is estimated as where ∆a d is the difference in a app between the void and control models, and t d is the cosmic age at redshift z d . The χ 2 contribution from all seven lenses is where H obs 0,lensing,i and σ lensing obs,i are the derived Hubble constant and corresponding uncertainty for lens system i as reported by Shajib et al. (2020) or Wong et al. (2020), which we summarize in Table 2.

Local Group peculiar velocity
An important constraint on our model is the observed motion of the LG relative to the surface of last scattering. The observed CMB dipole indicates that the LG moves with a peculiar velocity of v LG = 627 ± 22 km s −1 towards Galactic coordinates (l, b) = (276 • ± 3 • , 30 • ± 3 • ) (Kogut et al. 1993).
To calculate the expected peculiar velocity in different parts of the void, we first need to consider the motion of the void as a whole. The void peculiar velocity v void arises from the time-integrated EFE. Using the approach stated in section 2.2 of Banik et al. (2018), we get that where g ext is the external field, and the integrating factor a accounts for Hubble drag. The total velocity of a particle wrt. the CMB is where v int is the 'internal' velocity of the particle relative to the void centre, and θ is the angle between v void and the voidcentric position of the particle. A schematic representation of this situation is depicted in Figure 3. For numerical purposes, the simulated void is divided into cells. The volume of cell i is where ∆r and ∆θ are the radial and angular bin size, respectively. V i is determined by the change in r 3 and cos θ across the cell. In order to quantify how the observed v LG = 627 km s −1 affects the relative probability of a model, we define f motion as the proportion of cells which satisfy where is a numerical parameter whose choice should have no bearing on our final results. To get a good balance between reducing numerical noise and increasing the accuracy, we choose = 50 km s −1 . The resulting error should be of order (50/630) 2 , which is acceptable given other uncertainties.
On the other hand, 50 km s −1 is much larger than the change in v tot between adjacent cells. We obtain similar results if = 30 km s −1 is used instead. Using this discretized scheme, v int θ v void Figure 3. Schematic of the KBC void, which as a whole moves with velocity v void due to the time-integrated EFE accounting for Hubble drag (Equation 71, see also section 2.2 in Banik et al. 2018). The total peculiar velocity of a particle wrt. the CMB, v tot , is calculated by combining v void with the internal velocity v int of a particle relative to the void centre (Equation 72). The inner circle illustrates the region in which v tot ≤ v LG = 627 km s −1 . We estimate its volume by adding the volumes of the red cells.
we get that where r rms void is the rms size of the void, and n is a dimensionless factor of order unity that sets our prior expectation for how close we are to the void centre (we must be within a distance of nr rms void ). Since observations suggest that we are located quite close to the centre (e.g. Keenan et al. 2013;Karachentsev & Telikova 2018), we adopt n = 0.5 for our probability calculations. We estimate the void size as with the void profile δ (r) found using Equation 54. In practice, we cut off the integrals at a very large distance much beyond the possible extent of the void. Since δ → 0 at large r, this is sufficient to accurately estimate the limiting values of both integrals. In Section 4.2, we apply a less sophisticated probability calculation where we assume that v tot follows a Gaussian distribution. The extent to which the observed v LG is an outlier to the simulated v tot distribution is given by the proportion of the void volume with v tot ≤ v LG . This allows an easier comparison with the other constraints. Table 3 summarizes the here presented observational constraints, which are used to test our cosmological MOND model in the following section.

RESULTS OF MOND SIMULATIONS
In this section, we perform a detailed parameter study of our cosmological MOND models. This includes an estimation of the tension between our best-fitting model and observations of the local Universe. We focus on a void initialized with a Maxwell-Boltzmann profile (Section 3.2.1). Results for Gaussian and exponential starting profiles are presented in Appendix C. We first quantify the relative probabilities of different models (Section 4.1), and then check how well our best-fitting model agrees with observations (Section 4.2).

Relative probabilities of different models
The observational constraints can mostly be assumed to have a Gaussian distribution, allowing a standard χ 2 -based analysis. This is due to the central limit theorem and the fact that e.g. many SNe are used in the study of Camarena & Marra (2020a). However, the expected distribution of v tot (Equation 72) is based on just one void, so we cannot assume Gaussianity. This constraint requires a more careful treatment, as explained in Section 3.3.6.
Combining the different constraints, we get that the joint probability of each model is We use i to label different observational constraints, each of which has uncertainty σ obs,i . The only model-dependent uncertainties are the density contrasts of the inner and outer parts of the KBC void, a consequence of the applied RSD correction (Equation 57). Figure 4 shows the marginalized posterior distributions of the model parameters and parameter pairs based on 10 6 MOND models. The assumed external field strength has a significant impact on individual models. On the one hand, increasing the EFE typically makes the MONDian subsystem more Newtonian, suppressing the growth of structures. This results in less pronounced voids at the present time, and consequently a local Hubble constant and acceleration parameter closer to the Planck predictions. On the other hand, some EFE is required because otherwise structure formation would be too efficient, causing a very high local Hubble constant and acceleration parameter. These considerations restrict g ext to the range (0.024 − 0.076) a 0 at 2σ confidence, with the most likely value being 0.04 a 0 .
In contrast, the initial void size and strength are not strongly constrained − our analysis merely yields 2σ limits of α void = 10 −5 − 2.91 × 10 −4 and r void = (173.4 − 818.6) cMpc. This is because the local Hubble constant and acceleration parameter are estimated only for SNe in the redshift range 0.023 ≤ z ≤ 0.15, which does not constrain the outer part of the void. Although some constraints are available at higher z from lensed quasars, the uncertainties of H 0 measured in this way are relatively large, allowing for a wide range of possible model parameters ( Table 2).
The best-fitting model is that for which the joint probability (Equation 78) becomes maximal. We mark this as a red dot in Figure 4 and consider it our fiducial model. It has an initial void strength of α void = 3.76 × 10 −5 at z = 9, an initial void size of r void = 228.2 cMpc, and an external field strength of g ext = 0.055 a 0 , causing the void as a whole to move with v void = 1586 km s −1 . We analyse this particular model in more detail in the subsequent section.
The marginalized posterior distributions for MOND models with Gaussian and exponential initial profiles are shown in Appendix C. Those models still assume a timeindependent EFE. In Section 5.2.2, we present and discuss an analysis demonstrating that allowing time-dependence of the EFE reveals no strong preference for a time-varying EFE, though some variation is expected on theoretical grounds.

The fiducial model
In the following, we discuss the results of our best-fitting (fiducial) model.

Density profile
We begin by studying the density contrast of the fiducial model at different times. This is plotted in the lefthand panel of Figure 5. The void starts with an initial size of r void = 228.2 cMpc and a very small initial strength of α void = 3.76 × 10 −5 at z = 9. Equation 76 implies that r rms void = r void √ 3 = 395.2 cMpc at that time. At present, the void has grown to a size of r rms void = 528.7 Mpc and has a density contrast of δ in = 0.172 in the redshift range 0.01 < z < 0.07. Correcting the corresponding observed density contrast from Keenan et al. (2013) by the model-dependent RSD correction factor f in model = 1.38 yields δ in obs,corr = 0.254 ± 0.083. This agrees with the simulated value at the 0.99σ level. The calculated density contrast between 600 and 800 Mpc is δ out = 0.050, which also compares favourably with the RSDcorrected observed density contrast δ out obs,corr = −0.052 ± 0.105 ( f out model = 1.05). The tension in this case is only 0.97σ.  The red dashed, black solid, and black dashed contours mark the 1σ, 2σ, and 3σ confidence levels, respectively. For 1D posteriors, these are shown using the horizontal black lines. The 1σ and 2σ lines are at almost the same level for r void . The red dot or vertical line marks the best-fitting model with an external field strength of g ext = 0.055 a 0 , and an initial void size and strength of r void = 228.2 cMpc and α void = 3.76 × 10 −5 , respectively, at z = 9. This model is analysed in more detail in Section 4.2.
The long-range modification to gravity in MOND causes structure formation to be much more efficient than in ΛCDM cosmology (e.g. Sanders 1999;Famaey & McGaugh 2012, and references therein). This can be seen in the right-hand panel of Figure 5, which shows the density contrast of a 300 cMpc sphere for the best-fitting MOND model (the blue solid line) and two Newtonian models (the red lines, with shaded grey region between them) over cosmic time. The solid red and blue lines correspond to the same initial conditions, but end up with very different δ at the present time.
As expected, Newtonian models with different initial δ show a similar pattern of structure growth since they are all in the linear regime. The density contrast scales as δ ∝ ∼ a 3.8 and δ ∝ ∼ a 0.8 over the interval 0.3 ≤ a ≤ 0.7 in our bestfitting MONDian model and equivalent Newtonian model, respectively. The ΛCDM scaling is slightly < 1 because dark energy slows down the growth of structure at late times. The very rapid structure growth in our MOND model can be reduced by applying a higher EFE in the past. In the case of a time-dependent EFE with n EFE = −1 in Equation 50, the growth rate reduces to δ ∝ ∼ a 3.3 . The initial void strength must then be ≈ 13× larger (α void = 4.98×10 −4 ) to compensate for the higher EFE in the past (the blue dot-dashed line in the right-hand panel of Figure 5). Models with a timedependent EFE are discussed further in Section 5.2.2.

Hubble diagram
Our fiducial model yields H model 0 = 76.15 km s −1 Mpc −1 and q model 0 = 1.07. This is consistent with the observations of Camarena & Marra (2020a) at the 84.20% confidence level (only 0.20σ tension). The combined inference on both parameters is shown in Figure 6, which demonstrates that the best-fitting models with Maxwell-Boltzmann, Gaussian, and exponential initial profiles are all consistent with these observations within 1σ. Thus, we show for the first time that the Hubble tension can be resolved in MOND. Note that the Planck parameters (the green dot) are in ≈ 4.39σ tension with these local observations. Time delays from strong gravitational lenses also provide an important constraint on our model. We use Figure 7 to show H lensing 0,model in dependence of redshift, allowing a comparison with measurements from seven lens systems (Section 3.3.5). Interestingly, our model systematically underestimates H 0 especially at low redshifts, causing a 2.05σ tension with the observations of Wong et al. (2020) and Shajib et al. (2020). We expect that this discrepancy is partly caused by void motion due to the EFE, though there is also some internal inconsistency between the void profile of Keenan et al. (2013)   background H global 0 = 67.4 km s −1 Mpc −1 , since it is difficult to imagine a local void having substantial effects at z = 0.5. We discuss these issues in more detail in Section 5.2. Although it is expected that strong lensing in MOND occurs similarly to standard cosmology (Section 3.3.5), we redo our analysis without constraints from lensing-based H 0 measurements in Section 5.2.1.

LG peculiar velocity
Our model yields the total peculiar velocity wrt. the CMB in different parts of the void, as mapped in Figure 8. The entire void moves in the direction indicated by the arrow, which arises from the EFE (Section 3.2.2). Interestingly, the model allows for very high total and internal peculiar velocities, especially towards the void edge. We can also get partial or total cancellation between internal motions within the void and that of the void as a whole, creating a rather large region in which v tot ≤ v LG = 627 km s −1 (Kogut et al. 1993). This region is at a distance of ≈ (150 − 270) Mpc from the void centre, implying that the LG must be slightly offcentre. Applying Equation 75 to find the fraction this region represents of the whole void, we estimate that the observed v LG represents a 2.34σ outlier to the simulated v tot distribution, which causes therewith the highest tension amongst the here used observational constraints.

Overall agreement with observations
Finally, we quantify the combined tension of our fiducial MOND model with local observations. As discussed earlier, most observables can be treated using a standard χ 2 approach, but additional care is needed for v LG . Thus, we quantify the likelihood of different χ 2 , v tot combinations LG is probably near the right end of this curve because the observed radio dipole (Section 1.1) indicates that we are currently moving away from the void centre in the CMB frame. We show only half of the velocity map because v tot is axisymmetric about v void . The here shown total peculiar velocities are v tot < ∼ 0.01c, justifying the use of nonrelativistic equations for the void gravitational field (Section 3.2).
according to our fiducial model, with χ 2 found using Equation 78. We can then quantify the extent to which the actually observed combination is unlikely.
In our model universe, the joint probability that the observables can be summarized by some χ 2 , v tot combination is P (Observations|Best-fitting model) = P χ 2 · P motion (v tot ) , (79) where P( χ 2 ) is the probability density function for a χ 2 distribution with 8 degrees of freedom (i.e. 11 observational constraints and three model parameters). P motion (v tot ) is estimated from our simulation by splitting the volume into cells (Section 3.3.6) and assigning the volume in each cell to different bins in v tot , thereby building up a discretized picture of its distribution. This procedure does not assume that v tot follows a Gaussian. If our fiducial model is correct, χ 2 must arise solely from measurement errors, while the observed v LG reflects our position within the void. This causes that χ 2 and v tot have independent distributions, allowing them to be multiplied. We neglect the 22 km s −1 uncertainty in v LG (Kogut et al. 1993) because this is much smaller than the ≈ 4000 km s −1 range in v tot allowed by our model (Figure 8).
We use Figure 9 to show the joint χ 2 , v tot distribution based on our fiducial model. This explains the local observations at the 1.14% confidence level (2.53σ tension). Individual observational constraints are summarized and compared with observations in Table 4.
The χ 2 contributions from different constraints are visualized in Figure 10 as a pie chart, which also shows the Table 4. Comparison of individual local Universe observables with our fiducial MOND model (Maxwell-Boltzmann profile, g ext = 0.055 a 0 , r void = 228.2 cMpc, α void = 3.76 × 10 −5 , v void = 1586 km s −1 , r rms void = 528.7 Mpc, n EFE = 0). The last row shows the probability of a higher χ 2 given the number of degrees of freedom. We express this as the equivalent number of standard deviations for a 1D Gaussian (inverting Equation 24). For v LG , we show the proportion of the void volume where v tot ≤ v LG . Results for other void profiles are shown in Appendix C.
Parameter number of data points each constraint represents, and the corresponding level of tension. To facilitate a comparison with the other constraints, we use our previous estimate that v LG is a 2.34σ outlier to the simulated v tot distribution (Section 4.2.3). Therefore, we assign a χ 2 contribution of 2.34 2 to this constraint. The best-fitting MOND models with a Gaussian and an exponential void profile agree with the local observations at the 0.45% (2.84σ) and 0.34% (2.93σ) confidence level, respectively (see also Appendix C). Thus, our best-fitting void model with a constant EFE cannot be rejected regardless of the initial density profile. The implications of our results are discussed in Section 5.2, which also looks at the overall picture of the νHDM model and the theoretical uncertainties of the here applied MOND approach (Section 5.2.3).

DISCUSSION
We discuss what our results in Sections 2 and 4 imply for ΛCDM and MOND cosmologies. This is followed by a consideration of commonly proposed arguments claiming that a  Table 4). The bracketed numbers show the number of degrees of freedom and corresponding level of tension for each constraint. The value for the motion of the LG is estimated based on the fraction of the void volume for which v tot ≤ v LG = 627 km s −1 .
local underdensity cannot solve the Hubble tension (Section 5.3). Keenan et al. (2013) measured the K-band luminosity density as a function of redshift and found evidence for an underdensity around the LG with a radial extent of ≈ 300 Mpc (see their figs. 9 and 10). They used the 2M++ catalogue (Lavaux & Hudson 2011), which covers ≈ 90% of the sky based on photometric data from the 2MASS-XSC catalogue and redshift data from the 2MRS, 6DFGRS, and SDSS catalogues. From the K s < 13.36 luminosity density, Keenan et al. (2013) estimated a relative density contrast of δ ≡ 1 − ρ/ρ 0 ≈ 0.5 in the redshift range 0.0025 < z < 0.067 (pink down-pointing triangle in their figure 11). Probing the luminosity density slightly deeper (K s < 14.36) but only in the SDSS and 6DFGRS regions, they derived a slightly lower density contrast of δ = 0.46 ± 0.06 in the redshift range 0.01 < z < 0.07 (the light blue dot in figure 11 of Keenan et al. 2013). We used this value for the inner density contrast of the KBC void in order to minimize tension with the ΛCDM framework.

Assessing the tension for ΛCDM
Using the MXXL simulation (Angulo et al. 2012), we calculated the density contrast for spheres with an outer radius of 300 Mpc and an inner hole of radius 40 Mpc around 10 6 vantage points at z = 0. We also took into account the sky coverage of the 2M++ survey by generating at each vantage point a random observing direction from which 90% of the mock sky is covered (Section 2.1). The so-selected density fluctuations have an rms amplitude of 3.2%. This is consistent with the prediction of the Harrison-Zeldovich spectrum (Harrison 1970;Zeldovich 1972) in combination with the early universe normalisation of σ 8 = 0.811 ± 0.006 (Planck Collaboration VI 2020). Since Keenan et al. (2013) used a fixed distance-redshift relation (see their section 4.7), we applied an RSD correction to the simulated density fluctuations. The rms fluctuation then became 4.8%, with the individual values closely following a Gaussian of this width (Appendix A). Thus, the observational uncertainty of 6% is larger than the expected cosmic variance.
Based on our analysis, the observed KBC void is in 6.04σ tension with standard ΛCDM cosmology and cannot be explained with cosmic variance (Section 2.2.1). This contrasts with Sahlén et al. (2016), who concluded that supervoids such as the KBC void are consistent with standard theory. However, for this they used a top-hat galaxy density radius of 210h −1 Mpc and a DM density contrast of δ = 0.15 − 0.2. Their assumed δ describes a much less pronounced void than the observed density contrast derived from the 2M++ survey (Keenan et al. 2013). 8 Moreover, even a 15% true underdensity on a 300 Mpc scale is very difficult to reconcile with ΛCDM due to the expected variance being only 3.2% (Section 2.2.1). While it may be possible for such large voids to exist somewhere in the Universe, it would be unlikely for us to live inside one − unless they are more common.

Hubble tension
Several studies have already discussed a potential connection between the local void and the Hubble tension (e.g. Keenan et al. 2013;Enea 2018;Shanks et al. 2019;Kenworthy et al. 2019). Indeed, if mass conservation is assumed, a large underdensity in the local number density of galaxies should also show up in the velocity field. However, given the expected cosmic variance in ΛCDM, Figure 2 indicates that a ≈ 10σ density fluctuation would be necessary to explain the locally observed expansion rate within its 2σ confidence region (H local 0 = 73.8±1.1 km s −1 Mpc −1 , Riess et al. 2019;Wong et al. 2020). The maximum plausible 5σ density fluctuation is still not enough to explain H local 0 at the 5σ level. These findings are broadly consistent with Wu & Huterer (2017), who concluded that a void with δ = 0.8 and a radius of 120h −1 Mpc could resolve the Hubble tension. Such a void would be in ≈ 20σ tension with ΛCDM (Kenworthy et al. 2019). While this by itself does not constitute an argument against such a large local underdensity, Wu & Huterer (2017) stated that observations disfavour it. Section 5.3 of our work contains a more detailed discussion of claimed problems with a local void solution to the Hubble tension.
Combining the mutually consistent SH0ES and H0LiCOW results, Wong et al. (2020) showed that the Hubble tension has now reached the 5.3σ level. Thus, both the 8 In ΛCDM, the RSD effect implies δ ≈ 2 3 δ obs (Equation 22). Thus, δ obs = 0.46 does not correspond to a true underdensity of δ = 0.2.

Mass conservation
Modified gravity necessary z=1100 z=0 Figure 11. The KBC void in context. This large local underdensity and the Hubble tension can both be reconciled if mass is conserved in the Universe. However, such a large void cannot form out of the initial conditions of the CMB if the initial scale-invariant power spectrum (Harrison 1970;Zeldovich 1972) is preserved. In particular, the existence of the KBC void within the ΛCDM framework is ruled out at 6.04σ, as demonstrated in Section 2.2. Thus, modified gravity is required to explain low redshift observables (Section 4.2) if the initial conditions at z = 1100 are indeed set by the CMB (Section 5.3.8). Importantly, any solution must simultaneously solve both the Hubble and underdensity tensions in order to conserve mass.
KBC void and the Hubble tension falsify ΛCDM at > 5σ significance. The most likely explanation in the context of standard theory is that both are caused purely by measurement errors, since there is less cosmic variance than observational uncertainty in both the KBC data and the H 0 measurement. This is especially true for the latter − Equation 18 shows that density fluctuations of 3.2% would impact H local 0 by only 0.4 km s −1 Mpc −1 (see also Wojtak et al. 2014). Since galaxy counting and measurements of the local Hubble constant face rather different observational challenges, the measurement errors would be independent, yielding a combined tension of 7.75σ. Using a more rigorous estimation that allows for cosmic variance results in a 7.09σ falsification of the ΛCDM paradigm (Section 2.2.3).
Importantly, any cosmological model which solves the Hubble tension without addressing the void (or vice versa) would violate the assumption of mass conservation in the Universe. Thus, early dark energy models (e.g. Hill et al. 2020;Khoraminezhad et al. 2020) which simply increase H 0 at the background level by ≈ 10% would overestimate the local H 0 by about this much once the observed KBC void is taken into account (Section 1.1). Perhaps the most important implication of the KBC void is that the Planck value of H 0 is probably correct at the background level, since any attempt to substantially change it would likely cause the void-corrected local value to disagree with local observations if mass is conserved in the Universe. Therefore, the KBC void is a plausible explanation for the Hubble tension if we can preserve a Planck background cosmology but enhance the cosmic variance. A schematic that considers the KBC void in a broader context is presented in Figure 11. Starting from the initial conditions of the CMB at z = 1100 and assuming the scale-invariant Harrison-Zeldovich power spectrum to be valid at that time, we have shown that the existence of a KBC-like void at present is virtually impossible in a standard context. This indicates that a scale-invariant power spectrum is violated today − order unity fluctuations on a 10 Mpc scale (e.g. Mantz et al. 2015) do not give way to 3% fluctuations on a 300 Mpc scale. This is a very strong hint that the gravitational inverse square law has to break down. Thus, the spatial distribution of matter on both an 8 Mpc scale (Peebles & Nusser 2010) and on an ≈ 1 Gpc scale (this work) suggest a long-range enhancement to gravity.

Assessing the tension for MOND
MOND (Milgrom 1983) is a low-acceleration modification to gravity originally designed to explain galaxy rotation curves without CDM. It has enjoyed a great deal of predictive success in this regard (Section 1.3). Therefore, extrapolating MOND from kpc to Gpc scales could be a promising way to address large-scale challenges for standard cosmology such as massive high-redshift galaxy clusters (e.g. El Gordo, Katz et al. 2013) and supervoids.
In this context, we study the possible origin of the KBC void from small initial density fluctuations, and its impact on the local Hubble constant. Unfortunately, we do not presently have a large enough cosmological N-body or hydrodynamical MOND simulation to quantify the likelihood of the KBC void, as done for the ΛCDM framework with the MXXL simulation (Section 2).
We therefore used the νHDM framework (Section 3.1) to develop a semi-analytical MOND simulation in which the expansion history is a standard flat background cosmology consistent with the latest Planck data (Planck Collaboration VI 2020) − which should be a good approximation also in MOND (Section 3.1.1). We applied MOND only to the density deviations from the cosmic mean (Section 3.1.4). In Section 5.2.3, we discuss the possibility of a non-trivial coupling between density perturbations and the background.
Our main MOND models assume an initial Maxwell-Boltzmann density profile motivated by the radial density distribution of the Local Volume. Karachentsev & Telikova (2018) showed that the matter density within a sphere of r = 40 Mpc (r = 135 Mpc) around the LG is only Ω m,loc = 0.09 − 0.14 (Ω m,loc = 0.05 − 0.16). This is ≈ 2 − 3× lower than the cosmic mean density measured by Planck Collaboration VI (2020), confirming the existence of a large local underdensity (Section 1.1). Karachentsev & Telikova (2018) also showed that the density increases inwards for heliocentric distances < ∼ 40 Mpc (see their figure 3), justifying our choice of a Maxwell-Boltzmann void profile (Section 3.2.1). In addition, we also run void models initialized with a Gaussian and an exponential profile (Appendix B). In all cases, the void profiles are parametrized by an initial void size and strength at z = 9. The initial void strengths range from α void = 10 −5 − 10 −2 , with the lower limit based on the observed density fluctuations in the CMB. By the time that z = 9, we expect significantly larger perturbations. Our upper limit on α void is sufficient to capture the range of models preferred by our analysis (Section 5.2.2).
The EFE is strongly constrained in our models because it affects the formation of cosmic structure and thus internal velocities within the void, in addition to the void's motion as a whole (Section 3.3.6). Models with a very small EFE create extremely deep and extended voids, which disagrees with the density contrast of the KBC void − especially for its outer region. This also results in a much larger local Hubble constant than observed. Increasing the EFE leads to v void v LG , so the observed v LG can only be explained by nearly complete cancellation with a large v int . However, a strong EFE makes the MONDian system more Newtonian and suppresses therewith the growth of structure. Consequently, models with a very high EFE produce very shallow voids and a local Hubble constant similar to its global value, causing that v int is not large enough to cancel v void .
Our analysis for the Maxwell-Boltzmann profile restricts g ext to the range (0.030 − 0.053) a 0 at the 1σ confidence level. Models with a Gaussian and an exponential void profile prefer a slightly larger EFE, i.e. (0.054 − 0.094) a 0 ( Figure C1) and (0.054 − 0.092) a 0 ( Figure C2) at the 1σ confidence level, respectively. This is because the Maxwell-Boltzmann profile reduces δ near the void centre and therewith slows down the internal peculiar velocities of individual particles within the void. Thus, a lower EFE is required to achieve v tot ≤ v LG = 627 km s −1 over a large part of the void.
Our analysis rules out models without an EFE, which is in any case a logical consequence of the non-linearity inherent to Milgrom's law of gravity (Milgrom 1986). Observationally, MOND without the EFE is strongly disfavoured by the velocity distribution of wide binary stars in the Solar neighbourhood (Pittordis & Sutherland 2019). The EFE is also necessary to explain the internal velocity dispersions of dwarf galaxies (Section 1.3). We discuss the timedependence of the EFE in more detail in Section 5.2.2.
In contrast to the EFE, Figure 4 indicates that the initial void size and strength are not strongly constrained by observations. Thus, other initial void parameters could in principle also yield reasonable results at the present time. In particular, our analysis of Maxwell-Boltzmann voids yields 1σ confidence intervals on r void and α void of (173.9 − 636.9) cMpc and (1.07 − 8.12) × 10 −5 , respectively. Models with Gaussian or exponential initial profiles allow for larger voids, but with a similar void strength (Appendix C). There are two main reasons why both void parameters are only weakly constrained by local observations. First of all, H 0 and q 0 are derived from data in the redshift range 0.023 ≤ z ≤ 0.15 and constrain therewith only the inner and not the outer part of the void. Secondly, the uncertainties of H 0 measured from strong lens systems are relatively large, which allows for a wide range of possible void behaviours in the outskirts.
We found that our best-fitting Maxwell-Boltzmann MOND model requires an EFE of g ext = 0.055 a 0 , an initial void size of r void = 228.2 cMpc, and an initial void strength of α void = 3.76 × 10 −5 at z = 9. The EFE causes a bulk flow of v void = 1586 km s −1 at z = 0. Our fiducial model explains the local observations listed in Table 3 at the 1.14% confidence level (2.53σ tension). Figure 4 shows that models with somewhat different initial conditions also provide reasonable results.
The rms density fluctuation in the total matter field at the CMB (z = 1100) is δ rms ≈ 10 −4 (Section 3.1.3), which implies δ rms ≈ 10 −2 at z = 9 for a ΛCDM cosmology. Thus, α void of our best-fitting model is much lower than the expected cosmic variance in ΛCDM. This could make KBC-like voids very common in the universe, potentially conflicting with the observed foreground lensing of the CMB (Section 5.3.1). This problem could be alleviated if the EFE was stronger in the past (Section 5.2.2), or if the peculiar accelerations and the Hubble flow are coupled (Section 5.2.3) − both would slow down the growth of structure. It would be highly interesting to quantify the existence of KBC-like voids in a large cosmological MONDian N-body simulation, especially if it accounts for the HFE in some way.
As already shown in several previous studies (e.g. Sanders 1998;Famaey & McGaugh 2012, and references therein), we affirm that structure formation is much more efficient in MOND compared to the Newtonian case (righthand panel of Figure 5). Applying an RSD correction based on the best-fitting model, the observed underdensity is 25.4±8.3% (−5.2±10.5%) in the inner (outer) part of the KBC void. The enhanced growth of structure allows our fiducial model to match these constraints at the 0.99σ (0.97σ) confidence level. In contrast, the KBC void rules out the ΛCDM framework at 6.04σ (Section 2.2.1).
In MOND, the long-range modification to gravity causes a very shallow decrease of the density contrast with distance, causing our model to systematically underestimate the density at the outer part of the KBC void as derived from the K-band luminosity data of Keenan et al. (2013). However, observational uncertainties on the density contrast there are still relatively large. Future surveys would be necessary to more precisely measure the density profile beyond ≈ 400 Mpc. This may provide an important test of our model because the radial density profile should be sensitive to the underlying growth rate.
Interestingly, Angus & Diaferio (2011) found some evidence for large voids with a diameter of 250h −1 Mpc in their 512h −1 cMpc N-body cosmological MOND simulation with massive neutrinos. Although both large voids and massive galaxy clusters are expected in a MOND cosmology (e.g. Sanders 1998), it is not fully clear if those were formed artificially due to the low particle resolution . Their simulations also assume no coupling between peculiar accelerations and the Hubble flow (i.e. β = 0 in equation 8 of Sanders 2001). As discussed further in Section 5.2.3, a coupling to the Hubble flow would suppress the formation of massive voids and clusters on scales > ∼ 100 Mpc. Therefore, it would be very valuable to revisit their cosmological simulations with a higher particle resolution and an AMR grid code such as phantom of ramses (Teyssier 2002;Lüghausen et al. 2015). Such MOND simulations would require very large box sizes in order to include large-scale modes and the resulting EFE on smaller regions (Section 3.1.4). For very long modes, light travel time effects could be important such that a relativistic code is required. This could be based on the model of Skordis & Z lośnik (2019).
A unique characteristic of our void model is the prediction of very high total peculiar velocities, especially towards the void edge in the direction parallel to g ext . In the bestfitting model, the void as a whole moves with a peculiar velocity of v void = 1586 km s −1 due to the EFE from source(s) beyond the void (i.e. at z > ∼ 0.15). Thus, our model predicts a sphere centred on the LG should have a large bulk flow of ≈ 1000 km s −1 in a similar direction to g ext . This is qualitatively similar to the results of previous νHDM simulations (Katz et al. 2013). Interestingly, some evidence for a large bulk flow has been found (Kashlinsky et al. 2008(Kashlinsky et al. , 2011. We discuss this further in Section 5.3.2. Partial cancellation between the void's motion and internal motions within it leads to a region ≈ (150 − 270) Mpc from the void centre in which v tot ≤ v LG = 627 km s −1 . The fraction that this volume represents of the whole void corresponds to a 2.34σ event, implying that the LG is statistically not at a special position in the void. Note that the LG motion causes the highest tension amongst our constraints (Table 4 and Figure 10).
Our fiducial model gives an apparent expansion history very close to ΛCDM (Figure 12 Camarena & Marra (2020a) at the 84.20% confidence level (0.20σ tension). The best-fitting models with a Gaussian and an exponential void profile agree at the 0.83σ and 0.89σ level, respectively (Figure 6). Thus, we showed for the first time that the KBC void can arise in MOND and solve therewith the Hubble tension.
The locally observed acceleration parameter q 0 = 1.08 ± 0.29 (Camarena & Marra 2020a) disagrees with the ΛCDM expectation of q 0,ΛCDM = 0.53 (Planck Collaboration VI 2020) at the 1.9σ level. In combination with the H 0 discrepancy between these studies, this would falsify ΛCDM at 4.54σ confidence (see also Figure 6). Interestingly, q 0 > 1 is not possible for a standard background cosmology. The locally observed high Hubble constant and acceleration parameter provide compelling evidence that the Hubble tension is caused by a local effect like the KBC void. This addresses the concern of Kenworthy et al. (2019) that the KBC void is not evident in the SNe distance-redshift relation (Section 5.3.6) − both the first and second derivatives of the distance-redshift relation very much point to a local void. Observationally, a discrepancy could also appear as a third order effect in the jerk parameter j ≡ a 2 a/ a 3 , but given the already large uncertainty of q 0 , it would be difficult to measure j 0 precisely.
As discussed in Section 3.3.5, strong lensing does not occur in the MOND regime and so should be similar to in ΛCDM cosmology (Sanders 1999). Thus, our main analysis includes H 0 constraints from seven strongly lensed quasars as obtained by Shajib et al. (2020) and Wong et al. (2020). The latter work applied a blinded analysis (described in their section 3.6) and found that H 0 decreases as a function of lens redshift at 1.9σ significance (see their appendix A). H 0 becomes consistent with Planck expectations at z > ∼ 0.5, well beyond the void. This is again a very strong indication that the Hubble tension is driven by a local environmental effect such as the KBC void. A decrease of the inferred H 0 with redshift is also apparent in our MOND model (Figure 7) and is a generic consequence of any local resolution to the Hubble tension. The redshift dependence of H model 0,lensing depends mainly on the density profile of the void. For our fiducial model, the combined tension with all seven lensing-based H 0 measurements is χ 2 = 14.66, which represents 2.05σ tension for 7 degrees of freedom. In the case of a Gaussian and an exponential void profile, this would reduce to 1.76σ and 1.83σ, respectively, because the best-fitting models have a larger void (Appendix C). This discrepancy with the lensing data is mainly caused by a systematic underestimation of H 0 by our models, especially for the two lowest redshift lenses RXJ1131−1231 at z = 0.295 and PG 1115+080 at z = 0.311 (Section 5.2.1).
The high values of H 0 at z > ∼ 0.4 are also conspicuous because according to Keenan et al. (2013) the density should have reached the cosmic mean already at z ≈ 0.2 (see their figure 11). Consequently, we would expect that H 0 obtained from lenses located at z > ∼ 0.4 must be very similar to the Planck prediction. The H 0 values from the four highest redshift lenses of Wong et al. (2020) give a median (mean) of 71.3 km s −1 Mpc −1 (70.8 km s −1 Mpc −1 ), which is 3.9 km s −1 Mpc −1 (3.4 km s −1 Mpc −1 ) higher than the Planck prediction. This systematic offset can be reduced if the background H 0 is underestimated due to errors in the Planck measurements caused by intergalactic dust (Yershov et al. 2020; see also Vavryčuk 2018Vavryčuk , 2019. We discuss this issue further in Section 5.3.8. It is also possible that there is some systematic offset in H 0 measurements using strong lensing time delays (Kochanek 2020). Since we assume mass conservation in our models, this discrepancy cannot be fully resolved in our analysis − the algorithm searches for a compromise between the KBC and lensing data. It seems that there is some internal inconsistency between them. A strong test of our model would be to infer H 0 very accurately from nearby and high-redshift lens systems. The model predicts that H 0 should be almost identical to the Planck prediction at z > ∼ 0.9. However, a substantial anomaly is expected for a lens at z = 0.1. A measurement here would nicely complement the SNe results of Camarena & Marra (2020a), which go out to z = 0.15.

Excluding strong lensing time delays
Although strong lensing in MOND should be similar to the ΛCDM framework (Sanders 1999), the H 0 measurements from strong lensing could be affected by the EFE on the void. This requires a better understanding of its origin. We have assumed that the EFE in our simulations affects the void as a whole, implying that it must be sourced by something beyond the void, i.e. at z > ∼ 0.15. This approach would be valid for deriving H 0 and q 0 from SNe data in the range 0.023 ≤ z ≤ 0.15 since the EFE would similarly move everything in this region. However, this may not be true for even the lowest redshift lens as its z = 0.295 (Table 2). If the EFE is sourced by something at lower z, it would move the void − but not the lens.
In this context, we consider in more detail the geometry of the void and the sky positions of the lenses. The observed CMB dipole shows that the LG moves with a velocity of 627 km s −1 wrt. the CMB towards Galactic coordinates (276 • ± 3 • , 30 • ± 3 • ), which roughly matches the direction of the radio dipole (Section 1.1). Thus, the LG motion wrt. the CMB frame is probably directed away from the void centre (see also Figure 8). Interestingly, the two lowest redshift lens systems (RXJ1131−1231 and PG 1115+080) are located at Galactic coordinates (274.4 • , +45.9 • ) and (249.9 • , +60.6 • ), respectively, which roughly coincides with the directions of the CMB and radio dipoles. Thus, both lenses are also most likely located on the opposite side to the void centre. Assuming these lenses are not affected by the EFE on the void, this would cause an extra redshift. Consequently, a larger Hubble anomaly would be expected than calculated thus far.
Relative to the CMB dipole, the direction to each lens subtends an angle θ, where cos θ = 0.96 and cos θ = 0.82 for RXJ1131−1231 and PG 1115+080, respectively. The impact on the measured H 0 from these two low-z lenses can be estimated as where a is the scale factor at which the lens is observed, and ∆t is the corresponding lookback time. Thus, assuming the lenses are unaffected by the EFE, H 0 derived from the lens systems RXJ1131−1231 and PG 1115+080 would be overestimated by 1.10 and 0.88 km s −1 Mpc −1 , respectively. This would only slightly reduce the tension with our bestfitting models. Therefore, the sharp rise in H 0 values for the two lowest redshift lenses cannot be fully accounted for with the EFE and is still hard to explain. Additional lenses are needed to confirm this feature.
We address the possible impact of the EFE on the lensing-based H 0 measurements by redoing the analysis for our Maxwell-Boltzmann model without the constraints from strong lensing time delays. In this case, the bestfitting model has g ext = 0.050 a 0 causing v void = 1442 km s −1 , r void = 208.4 cMpc, and α void = 2.15 × 10 −5 . This model yields H model 0 = 76.47 km s −1 Mpc −1 and q model 0 = 1.21 (0.26σ combined tension), δ in = 0.167 (δ in obs,corr = 0.258 ± 0.082; 1.11σ), and δ out = 0.037 (δ out obs,corr = −0.038 ± 0.104; 0.73σ). The fraction of the void with v tot ≤ v LG represents a 2.25σ tension. The model explains all these local observations at the 1.96σ (5.0%) confidence level. Thus, excluding the lensing data allows for a somewhat better overall fit, but has little effect on the preferred model parameters.

Structure formation and external field history in MOND
Since the EFE acting on a MONDian subsystem depends on surrounding structure and therewith on the scale factor, we made in Section 3.2.2 the ansatz g N,ext (t) = g N,ext t 0 a n EFE (t).
So far, we have restricted attention to the case n EFE = 0. Letting n EFE vary in the range (−2, 2), Figure 13 shows its marginalized posterior based on 9 × 10 6 different models.
In the case of a Maxwell-Boltzmann profile, the analysis yields n EFE > −0.62 at the 1σ confidence level. Gaussian and exponential initial profiles allow for (−1.60, +0.43) and (−1.59, +0.52) at the 1σ level, respectively. Thus, only the Maxwell-Boltzmann profile prefers a weaker EFE in the past, while the other profiles prefer the opposite. A timeindependent EFE (n EFE = 0) − assumed for all our models thus far − lies within the 1σ range for all three considered void profiles. This justifies our assumption of a timeindependent EFE in our main analysis. The marginalized posterior distribution for n EFE and the initial α void are shown in Figure 14. As expected, a stronger EFE in the past (n EFE < 0) requires a deeper initial void. In particular, values of α void up to ≈ 10 −3 are now allowed, contrary to the case where we fix n EFE = 0 (Figure 4). This would be closer to the expected density fluctuations in ΛCDM when a = 0.1, since the ≈ 10 −5 density fluctuations in the CMB should have grown ≈ 100×. This is only mildly disfavoured by our analysis.
We can estimate the external field history based on the parameters of our fiducial model and the gravitational acceleration at recombination: where t CMB = 380 kyr, and t 0 = 13.8 Gyr. As explained in Section 3.1.3, g N t CMB ≈ 21 a 0 , causing MOND to have only a very small effect before recombination. Since our fiducial model prefers an external field of 0.055 a 0 , we use this as our estimate of g t 0 , which via Equation 8 yields g N t 0 = 0.0029 a 0 . Comparing the expected large-scale Newtonian gravitational fields at these times gives n EFE = −1.27 (the red vertical line in Figure 13). This is consistent with the marginalized posterior distribution for the Gaussian and exponential profile models at the 1σ level and for the Maxwell-Boltzmann models at the 2σ level. The Gaussian and exponential models prefer n EFE ≈ −0.5. Equation 83 implies n EFE = −1.27, which is a faster decline than the n EFE = −1 suggested by Equation 37 for linear perturbations in the matter-dominated era (Section 3.1.3). This could be due to the effect of dark energy, which is not taken into account in Equation 37 as it only dominates at z < ∼ 0.8. Since dark energy slows down the growth of structure, n EFE would be shifted to more negative values − if the fractional density perturbations are frozen in co-moving coordinates, we would get n EFE = −2.
Though it is beyond the scope of our semi-analytic study, we mention briefly that as structure grows in a MON-Dian universe, it imposes an external field on surrounding structures, thereby hampering their growth. This leads to structure formation in different regions becoming mutually correlated. In particular, since MOND gravity declines as 1/r whereas Newtonian tides scale as 1/r 3 (there being no EFE), any structure in a Milgromian universe is affected by much more distant structures compared to the Newtonian case. This makes it difficult to conduct a MOND simulation with sufficiently large volume to satisfy the CP.

Theoretical uncertainties in the MOND approach & outlook for further studies
At present, it is not known if MOND is related to a fundamental (quantum) theory (i.e. 'FUNDAMOND', Milgrom 2020a,b). As a result, we do not have a completely secure understanding of cosmology and structure formation in MOND. In fact, it is likely that the implications on these scales are not uniquely derivable from the RAR in disc galaxies. Although the relativistic MOND theory of Skordis & Z lośnik (2019) seems quite promising, its consequences for cosmology are not yet established. Therefore, the here applied cosmological model required us to make some ansatzes, whose uncertainties will be summarized and discussed in the following (see also Section 3.1). Motivated by previous theoretical studies, we assumed that the background cosmology in a Milgromian universe obeys the same Friedmann equation as in ΛCDM (Section 3.1.1). While this is not necessarily true in MOND, the observations of Joudaki et al. (2018) suggest that this works well empirically. Moreover, constraints from BBN imply rather small deviations from the standard a (t) during the radiation-dominated era. In a MOND context, this forces the expansion history to obey the standard Friedmann equation to sub-per cent precision in the matter and Λdominated eras . Moreover, the source term for the Friedmann equation remains the same if CDM is replaced with the same density in sterile neutrinos with m ν s = 11 eV/c 2 since both would be non-relativistic up to  Figure 14. Marginalized posterior distribution of the indicated model parameters based on 9 × 10 6 MOND models for a Maxwell-Boltzmann (left), Gaussian (middle), and exponential (right) initial profile. The red dashed, black solid, and black dashed lines mark the 1σ, 2σ, and 3σ confidence levels, respectively. A stronger EFE in the past (n EFE < 0) requires a stronger initial void strength at z = 9. The red dots mark the best-fitting models: g ext = 0.030 a 0 , r void = 218.3 cMpc, α void = 7.56 × 10 −4 , n EFE = −2 (Maxwell-Boltzmann profile, left-hand panel); g ext = 0.065 a 0 , r void = 1030.0 cMpc, α void = 1.07 × 10 −4 , n EFE = −0.5 (Gaussian profile, middle panel); and g ext = 0.070 a 0 , r void = 1030.0 cMpc, α void = 1.75 × 10 −4 , n EFE = −0.5 (exponential profile, right-hand panel).
very high z z eq . Consequently, there is very good reason to suppose that the background a (t) is very nearly the same as in ΛCDM.
Sterile neutrinos would only slightly affect primordial nucleosynthesis and plasma physics prior to recombination (Sections 3.1.2 and 3.1.3, respectively, see also figure 1 in Angus 2009). Given also a very nearly standard expansion history and high peculiar gravitational accelerations at that time, we expect the νHDM framework to yield the same CMB power spectrum as ΛCDM. The angular diameter distance to the CMB would also be the same in both frameworks, causing both to suffer from the Hubble tension if H local 0 is little affected by cosmic variance. Our main argument is that this last assumption holds in ΛCDM but not MOND.
To simulate structure formation in MOND, we made the usual ansatz that MOND should be applied only to the density deviations from the cosmic mean (e.g. Llinares et al. 2008;Angus & Diaferio 2011;Angus et al. 2013;Katz et al. 2013;Candlish 2016). We justified this in Section 3.1.4 based on the fact that Sanders (2001) showed in his section 2 that this approach (elaborated further in Sanders & McGaugh 2002) is valid in a non-relativistic Lagrangian that has the MOND behaviour. This so-called Jeans Swindle (Binney & Tremaine 1987) is one of the strongest assumptions of current cosmological MOND models. Falco et al. (2013) has formally shown that it can be justified in an expanding General Relativistic universe, but this needs to be mathematically demonstrated for a MONDian framework in which the Poisson equation is non-linear. Despite this uncertainty, it seems inevitable that structure formation would be significantly faster in MOND compared to ΛCDM on a 100 Mpc scale. This is because 100 Mpc is much larger than the free streaming length of both sterile neutrinos and CDM, so the only major difference between the ΛCDM and νHDM frameworks is a different gravity law. Since the accelerations are much smaller than a 0 (Section 1.3), we expect any MOND theory to yield a significant enhancement to the gravity generated by density perturbations. As a result, we argue that MOND models naturally possess the ability to explain the KBC void and Hubble tension.
We now discuss whether this conclusion remains valid  Figure 15. Accelerations along the trajectory of a particle ending up 300 Mpc from the void centre in our best-fitting model. We show the time evolution of g void (the blue solid line), g ext = 0.055 a 0 (the red horizontal dashed line), and |g Hubble | for a standard background cosmology (the green-dotted line). The thick horizontal grey line refers to a 0 . Note that g Hubble changes sign when a = 0.61 (z = 0.63) − it is positive at later times and negative earlier on. At present, g void ≈ 0.1 a 0 and g Hubble ≈ 0.2 a 0 , so the latter would limit the MOND boost to gravity at the void edge in case of a strong HFE.
if the Jeans swindle approach is not applicable because of the HFE, a coupling between the background cosmology and structures within it. As mentioned in Section 3.1.4, Sanders (2001) developed a non-relativistic two-field Lagrangianbased theory of MOND that couples the Hubble flow and the peculiar accelerations from inhomogeneities. This coupling is described by the adjustable parameter β in his equation 8, and is elaborated further in Sanders & McGaugh (2002). If β = 0, the coupling between these fields vanishes, which is equivalent to the above-mentioned Jeans swindle. In the case β 0, the background cosmology would remain intact, but the Hubble flow acceleration (g Hubble in Equation 38) would appear as an additional source of gravity that suppresses the ν factor in Equation 7.
As argued in Section 2.2.2, the existence of the KBC void and the Hubble tension can be simultaneously reconciled in the ΛCDM framework due to mass conservation, but only for a 10σ density fluctuation (Figure 2). We demonstrated that in the absence of any HFE, our best-fitting model can explain both the KBC void and the Hubble tension because of an enhanced cosmic variance compared to ΛCDM (Section 4.2). A coupling to g Hubble would reduce the cosmic variance in our MOND model, and therewith also the frequency of KBC-like voids. To estimate the possible impact, we use Figure 15 to plot the various accelerations entering Equation 38 over time for the same test particle presently 300 Mpc from the void centre. MOND should boost structure formation somewhat at all epochs with a > ∼ 0.2 since |g Hubble | < a 0 and the other acceleration terms are even smaller. However, in order to obtain a lower limit on the cosmic variance expected in the MOND framework, we assume the most conservative scenario in which structures grow only as fast as in the ΛCDM framework δ ∝ ∼ a 0.8 at all epochs when g Hubble is dominant − it is after all unclear exactly what would happen then. As shown in the right-hand panel of Figure 5, density fluctuations in our best-fitting model grow as δ ∝ ∼ a 3.0 for the period 0.5 ≤ a ≤ 0.8. During this period, our previous approach should be valid because |g Hubble | < g void (Figure 15). Thus, the 3.2% standard ΛCDM cosmic variance is increased by a factor of at least (0.8/0.5) 3.0−0.8 ≈ 2.8, implying a 9.0% cosmic variance. Since our best-fitting MOND model yields a present underdensity of δ in = 0.172 on a 300 Mpc scale 9 , it requires an ≈ 1.9σ fluctuation. In a MONDian model with n EFE = −1 (Equation 50), density fluctuations grow slower as δ ∝ ∼ a 2.8 for 0.5 ≤ a ≤ 0.8, implying an 8.2% cosmic variance and thus a 2.1σ fluctuation. The latter scenario may be more realistic in light of our discussion in Section 5.2.2. Since the present value of δ is essentially fixed by the observations of Keenan et al. (2013), the slower structure growth in such a model implies a deeper void in the past, increasing g void relative to our fiducial model. This would increase the timespan during which g Hubble is sub-dominant, thereby allowing MOND to enhance structure growth to a greater extent. We conclude that even in the case of a strong HFE and making very conservative assumptions, our MOND approach still succeeds in explaining the KBC void, which by mass conservation also resolves the Hubble tension. This is mainly because g Hubble was much smaller a few Gyr ago − it is ∝ a (Equation 38), which crossed zero at a lookback time of ≈ 6 Gyr ( Figure 15). r was also smaller in the past if we consider the same co-moving scale. While g Hubble can undoubtedly suppress the growth of structure, it was sub-dominant for an extended period − during which it should be appropriate to neglect the HFE. Indeed, a much more rapid growth of perturbations around the time a = 0 is also evident in figure 1 of Sanders (2001), who considered MOND with a strong HFE. Thus, the presence of dark energy can actually promote the growth of structure in a MOND context by reducing |g Hubble |. Observations of structure growth on a > ∼ 100 Mpc scale at a ≈ 0.61 might reveal evidence for this growth spurt, which would occur at a redshift beyond the extent of the KBC void. Sanders & McGaugh (2002) showed that since g Hubble scales directly with the size of a system, perturbations on smaller length-scales would be shielded from the HFE, allowing them to grow much faster in case of a strong coupling (see their figure 12). For instance, g Hubble ≈ 7 × 10 −3 a 0 for a particle at a scale of 10 Mpc, so g Hubble would be very sub-dominant compared to typical external fields of a few percent of a 0 (e.g. Famaey et al. 2007). This justifies the Jeans swindle approach used in numerical simulations if the goal is to address problems on this scale.
The initial α void required by our fiducial model is lower than the expected rms density fluctuations at z = 9 for a ΛCDM cosmology, which we can estimate by scaling the present value of 0.032 by 0.1 0.8 to obtain ≈ 0.005. Our results in Figure 14 show that this remains true even with a stronger EFE in the past (n EFE < 0), implying that KBC-like voids would be quite common in MOND. Such a void could explain the Cold Spot in the CMB, which is often interpreted as a huge underdensity (Nadathur et al. 2014). However, such structures are quite rare, suggesting that large voids are not very frequent in the Universe (Section 5.2). Moreover, a high frequency of KBC-like voids could cause too much foreground lensing of the CMB − though this is far from clear (Section 5.3.1).
In principle, the implications of MOND on large scales depend on an adjustable parameter analogous to β in equation 8 of Sanders (2001), which can be used to alter the frequency of KBC-like voids or massive galaxy clusters such as El Gordo at high redshift. A strong HFE would reduce the frequency of such structures. We apply Occam's Razor and assume β = 0 since there is no compelling observational or theoretical evidence for β 0. In particular, it is not yet clear if covariant theories (such as that proposed recently by Skordis & Z lośnik 2019) have any flexibility in the coupling between peculiar and Hubble flow accelerations, at least when we impose other constraints, e.g. that gravitational waves travel at c. Since g Hubble acts to suppress the void gravity and increases with distance, the outer density profile of the KBC void could empirically constrain the coupling. Our best-fitting model implies that at present, g ext g void a 0 at 300 Mpc, implying that g void ∝ ∼ r −1 in the outer part of the void. Since g Hubble > g void , a strong background coupling would make the system more Newtonian, causing therewith a steeper decline of gravity with distance. However, current measurements of the KBC void's outer density profile are not sufficiently precise to strongly constrain the HFE.
If a particular cosmological MOND model predicts that structure formation on 300 cMpc scales is very similar to standard cosmology, such a model would not be able to account for the KBC void − and would have to be rejected for similar reasons to ΛCDM (Section 2). Any viable covariant formalism of MOND has to describe the density and velocity field of the local Universe. Interestingly, the here applied approach can reproduce the KBC void, which solves therewith the Hubble tension due to mass conservation. Therefore, our model can serve to guide further theoretical development of MOND in a cosmological context.
The g Hubble term would be much larger in the CMB era. For the sound horizon scale at recombination (147.09 ± 0.26 cMpc, Planck Collaboration VI 2020), |g Hubble | ≈ 35000 a 0 . Thus, even a slight background-perturbation coupling would completely eliminate any MOND effects, making the pre-CMB universe purely Newtonian. This is also true in the absence of any such coupling, since the gravity from inhomogeneities is ≈ 21 a 0 (Section 3.1.3). Thus, our conclusions regarding the CMB are not affected by a possible HFE. This is also true for the model of Zhao (2008), which implies that the MOND acceleration scale a † is timedependent, with At the time of recombination, a † = 36500 a 0 , which exceeds |g Hubble | at that epoch. Thus, even in the presence of a very strong coupling of perturbations to g Hubble , his model implies significant MOND effects in the CMB, making it very difficult to fit its power spectrum. Furthermore, Milgrom (2017) showed that the model of Zhao (2008) is in tension with the rotation curves of galaxies at high redshift.
In this contribution, we assumed that Milgrom's constant a 0 is constant in space and cosmic time. While the former is expected in a fundamental theory, a time-varying a 0 is in principle quite possible. Although this is observationally not supported at the moment (e.g. Milgrom 2017), a 0 ≈ cH 0 /(2π), which could be a hint that MOND is fundamentally related to cosmology. If this relation is true, a 0 would decrease over cosmic time, implying that the early universe was more MONDian than assumed in our models. Thus, postulating that a 0 ∝ H would cause strong MOND effects in the CMB, arguing against the idea (Section 3.1.3). On the other hand, a much lower a 0 in the past would significantly raise the MOND timing argument mass of the LG, which for a constant a 0 ends up rather similar to its baryonic mass (Banik et al. 2018). Since HDM should not significantly cluster on such a small scale to avoid disrupting MOND fits to galaxy rotation curves, it appears that very strong time evolution of a 0 in either sense is ruled out empirically if not theoretically.
MOND as currently understood cannot explain the CMB power spectrum and the dynamics of galaxy clusters without an extra matter component. Therefore, we follow Angus (2009) in postulating the existence of HDM. This is another strong assumption of our model − but not more hypothetical than the existence of CDM. If anything, sterile neutrinos have been described as "almost part of the standard model" of particle physics (Merle 2017), while CDM particles are generally thought to require supersymmetry. In future, it will be very important to directly search for sterile neutrinos in terrestrial experiments.
In addition, precise measurements of the CMB power spectrum at > 4900 could put strong constraints on our model. This is because small shifts to the CMB power spectrum may arise when applying MOND supplemented with sterile neutrinos rather than Newtonian gravity with CDM. We have assumed that the effects are either not detectable or can be compensated through small adjustments to the cosmological parameters.
Furthermore, we modelled the gravitational field of the void using non-relativistic equations. This should be quite accurate because the total peculiar velocities are v tot < ∼ 0.01c within ≈ 250 Mpc of the void centre ( Figure 8). Moreover, the void is much smaller than the cosmic horizon. Gravity travel time effects should thus not be too significant if gravitational waves travel at c, as occurs in the model of Skordis & Z lośnik (2019).
The exact density profile of the KBC void is not fully known, so further assumptions are required when modelling it. Motivated by observations of an increasing density as one goes inwards for distances < ∼ 40 Mpc (Karachentsev & Telikova 2018), we assumed an initial Maxwell-Boltzmann profile (Section 4). We demonstrated the robustness of our results by also implementing Gaussian and exponential void profiles (Appendix C). The best-fitting models with those profiles yield a slightly larger overall tension compared to our main analysis, with the best parameters shifting to a stronger EF and larger initial void with comparable depth (Table C1). Thus, other void profiles could yield reasonable results with adjusted EF and void parameters. This issue could be constrained with better knowledge regarding the exact density profile of the KBC void.
Although the here presented MOND approach suffers from theoretical uncertainties especially with regards to the HFE, the encouraging results of our best-fitting models (Section 4.2 and Appendix C) suggest that our assumptions are reasonable. Furthermore, our models allow a wide range of possible void parameters (Figs. 4,C1,and C2), so adjusting these could in principle compensate theoretical uncertainties. In particular, a stronger HFE would require an initially deeper void − though we argued that the required depth would be reasonable even under conservative assumptions.
Once it is clear which covariant MONDian framework should be applied, the role of the HFE (if any) would become apparent. It would then be valuable to statistically quantify the existence of the KBC void within a large numerical MONDian cosmological simulation, enabling a comparison with our analysis of its likelihood in ΛCDM (Section 2). Such a simulation would also deliver a better understanding of void profiles at low redshift, and on how the growth of structure is regulated by the EFE from surrounding structures (Section 5.2.2).

Claimed problems for a local void solution to the Hubble tension
In the following, we address some commonly used arguments for why a void model cannot resolve the Hubble tension.

Other anomalies in large-scale structure
If the growth of structure is much more rapid than predicted by standard theory, large underdensities such as the KBC void should also exist at higher redshift. Large voids are not evident in the galaxy two-point correlation function, suggesting that large-scale structure seems to be consistent with the ΛCDM paradigm. However, it must be borne in mind that the underlying matter density field is not measured directly − it is estimated from the distribution of galaxies. At large distances, only the brightest galaxies can be observed, so one has to assume the so-called bias factor: which relates the galaxy density contrast δ galaxy to that of the underlying matter distribution. This bias factor is typ-ically chosen to match the ΛCDM expectation for cosmic variance at the relevant scale. Thus, an accurate modelindependent estimation of the density contrast can only be achieved for low-mass galaxies observed in the NIR, for which b ≈ 1 on 100 Mpc scales in any cosmological model. This makes it very difficult to perform a similarly detailed analysis to Keenan et al. (2013) for z > ∼ 0.5. Even the 2MASS survey they used only covers 57 − 75% of the luminosity function (see their figure 9). This fraction would be much lower for higher redshift galaxy samples.
We are also faced with the problem that galaxy positions are in general unknown − in the distant Universe, only redshifts are available. Even the redshifts are often not measured directly but are estimated photometrically. This leads to a significant smearing effect along the line of sight, making it rather difficult to identify distant supervoids (DES Collaboration 2019). The situation is reminiscent of the LG satellite planes − due to distance uncertainties, it is difficult to know if the satellites of a distant galaxy are distributed anisotropically. In both cases, if similar anomalies had not been reported at larger distances, this would not tell us whether such anomalies exist. 10 However, supervoids identified in the Dark Energy Survey do seem to show an enhanced ISW effect (DES Collaboration 2019). Their stacked analysis of 87 supervoids found that the effect has an amplitude of 5.2 ± 1.6 times the conventional expectation when combined with the earlier results of Kovács (2018).
Another possibly related anomaly is that the lensing amplitude implied by the CMB power spectrum is stronger than predicted (Di Valentino et al. 2020a). They suggested that this problem could instead be an indication that the Universe has a positive curvature. This would have serious implications for our entire understanding of the Universe and require a completely different cosmological model. For instance, a closed universe would imply a very low H 0 of 54 +3.3 −4.0 km s −1 Mpc −1 , which is completely inconsistent with local measurements (see their figure 7). However, the enhanced lensing amplitude − evident also in Di Valentino et al. (2020b) − could be the imprint of unexpectedly large density fluctuations caused by more rapid growth of structure. In this scenario, a supervoid would be a more likely explanation for the CMB Cold Spot (Nadathur et al. 2014). Indeed, their suggested void profile has a central underdensity of 0.25 and characteristic size of 280 Mpc, rather similar to the KBC void (see their equation 1). They concluded that such a void is highly unlikely in ΛCDM, as also implied by our results in Figure 1. It would be very valuable to empirically determine the actual frequency of such voids.

High peculiar velocities
Our fiducial model implies the existence of void regions with total peculiar velocity v tot ≤ v LG = 627 km s −1 relative to the surface of last scattering. Such low velocities are unlikely but allowed at the 2.34σ level.
Interestingly, our model predicts v tot of up to ≈ 4000 km s −1 for objects ≈ 250 Mpc from the void centre in the direction of its motion (Figure 8). Such high peculiar velocities are potentially detectable with the kinematic Sunyaev-Zeldovich (kSZ) effect of galaxy clusters (Sunyaev & Zeldovich 1980). However, this is difficult to disentangle from the thermal Sunyaev-Zeldovich (tSZ) effect because our model predicts similar peculiar velocities to the internal velocity dispersions of galaxy clusters. A large local underdensity also reduces the number of clusters available for kSZ studies, increasing the uncertainty. Hoscheit & Barger (2018) concluded that the KBC void is consistent with the linear kSZ effect (see their figure 6). Similar results were obtained by Ding et al. (2020). Some evidence for a bulk flow of ≈ 1000 km s −1 has been found (Kashlinsky et al. 2008(Kashlinsky et al. , 2011. This is broadly consistent with the expected motion of the whole void due to the EFE (v void = 1586 km s −1 ), though the bulk flow of a smaller region will depend on our exact location within the void and the survey volume. Bulk flows of ≈ 1000 km s −1 are not possible on a 100 Mpc scale in a ΛCDM universe, but would be expected in MOND (Katz et al. 2013).

Gravitational redshifting of the CMB monopole
A large local underdensity like the KBC void should also affect the mean temperature (monopole) of the CMB. This is because the height of the potential at our location causes a gravitational redshift. Using Equation 52, the general relativistic redshift for a photon travelling uphill from distance r to the centre becomes In the best-fitting MOND model, we obtain z GR = 8.4 × 10 −3 for the most distant test particle, which is 700× larger than the 1σ rms fluctuations of z GR = 1.2 × 10 −5 assumed in the study of Yoo et al. (2019). Their figure 2 shows that the impact of such a gravitational redshift on the inferred cosmological parameters is very small, even with an extra factor of 700. Moreover, we expect that the actual gravitational redshift at our position in the void should be much smaller. This is because we are not exactly at the centre of the void, and thus not at the highest part of its gravitational potential hill (Figure 8). Redshifting from the void's gravity would also be partially counteracted by the EFE, which is required in order to explain the rather slow motion of the LG wrt. the CMB (Section 4.2.3). None the less, it is possible that gravitational redshifting of the CMB would change the best-fitting HDM and baryon fractions by a few times their official uncertainties. Since these are nowadays rather small (Planck Collaboration VI 2020), we conclude that this effect has only a small impact on the CMB power spectrum, which moreover could probably be compensated through slight adjustments to the cosmological parameters. Of particular relevance for the Hubble tension is that gravitational redshifting of the CMB has a negligible impact on the precisely measured angular scale of the first acoustic peak (figure 2 of Yoo et al. 2019).

Assumption of Newtonian gravity for the void dynamics
Applying Newtonian gravity to the dynamics of any void would lead to sharp gradients in the predicted density and velocity profiles due to the steep inverse square law. Kenworthy et al. (2019) found no evidence for such a sharp edge in the SN luminosity-distance relation, which would − according to them − rule out the existence of a large local void with δ > 0.2 at the 4σ − 5σ confidence level. Also, Hoscheit & Barger (2018) applied the large scale void radial profile of Keenan et al. (2013) to show that H 0 is 1.27 ± 0.59 km s −1 Mpc −1 higher in the redshift range 0.0233 < z < 0.07 compared to 0.07 < z < 0.15. This only modestly reduces the Hubble tension, e.g. with the SNe data of Riess et al. (2016). However, as shown in Section 2.2, an inverse square law is too weak to produce a deep and extended underdensity like the KBC void. Therefore, the assumption of Newtonian gravity for the void dynamics is not sustainable. In MOND, the long-range modification to gravity would cause a much more gradual return from the void-induced peculiar velocities to the background cosmology, as demonstrated in Figure 12. Therefore, sharp features in the density profile and Hubble diagram are not expected in a MONDian model. This holds especially for H 0 derived from SNe data because in order to constrain the cosmological model, one has to consider many individual SNe. Consequently, the inferred H 0 only gradually declines towards the Planck prediction as SNe beyond the void are included in the analysis (as e.g. done by Colgáin 2019).

Restrictive upper limit on the void size
The present void size can be treated as a model parameter independently of the applied gravity theory. Adopting a very low upper limit on the allowed void size would unavoidably cause sharp features in the density and velocity profiles in any framework. Moreover, the Hubble tension cannot be resolved by a small void unless we postulate that it is extremely deep. This issue affected the analysis of Wu & Huterer (2017), who assumed a void size of 180 Mpc. They noticed that since the SNe data go out much further, it is difficult for such a small void to resolve the Hubble tension. None the less, they did not consider a larger void, opting instead for a very large density contrast of δ = 0.8. This led to poor agreement with direct measurements of the density field. However, a larger and shallower void would have provided much better agreement with observations, as shown in this work.
Using the high-resolution ΛCDM N-body cosmological simulation called Millennium-II, Xie et al. (2014) obtained that ≈ 14% of LG-like systems are located in a region that resembles the observed local void. Thus, they concluded that "the emptiness of the Local Void is indeed a success of the standard ΛCDM theory." However, by "Local Void", they meant a sphere of radius ≈ 8 Mpc, which is much smaller than the KBC void. Thus, their work cannot be used as an argument that the local void observed by Keenan et al. (2013) is consistent with ΛCDM cosmology, as is done in section 5 of Sahlén et al. (2016).  Figure 16. Difference in the redshift between our best-fitting model and a standard void-free cosmology (latter subtracted), shown as a function of lookback time. The first two dashed vertical lines show the range 0.023 ≤ z ≤ 0.15 covered by Camarena & Marra (2020a). The dotted vertical line marks z = 0.5. For comparison, the red points and green squares show fractional distance uncertainties using SNe (table 6 of Riess et al. 2018) and BAO (table 8 of Alam et al. 2017), respectively. Since fractional redshift errors are generally rather small, this gives an estimate of the uncertainty in the inferred Hubble constant, which has similar sensitivity to redshift and distance.

Fixing the acceleration parameter
The acceleration parameter q 0 describes the second time derivative of the scale factor (Equation 65). It is therefore a measure of the void's gravity. As discussed in Section 5.3.4, Kenworthy et al. (2019) concluded that the KBC void is not evident in the SN luminosity-distance relation. In addition to assuming Newtonian gravity, they fixed the acceleration parameter to the ΛCDM prediction of q 0 = 0.55. In general, q 0 would have a higher value if there is a large local void. To allow for this possibility, q 0 must be treated as a free parameter when using the apparent expansion rate history to constrain the properties of a local void.
Fortunately, Camarena & Marra (2020a) address this shortcoming by deriving q 0 and H local 0 jointly from SNe data without a restrictive choice of prior. Their analysis yields q 0 = 1.08±0.29, much higher than in the Planck cosmology. A high q 0 is also evident when using BAO data or treating the SNe Ia absolute magnitude as a free parameter (Camarena & Marra 2020b). This is a strong hint for the existence of a local void independently of the galaxy luminosity density (e.g. Keenan et al. 2013). Indeed, Colgáin (2019) inferred a local underdensity at z < ∼ 0.15 using SNe data alone. Both the local Hubble constant and acceleration parameter can be explained in our best-fitting MOND models ( Figure 6).

Effect of the void at high redshift
The apparent expansion rate history in our fiducial MOND model is very similar to the Planck cosmology ( Figure 12). None the less, the fractional difference in z between void and void-free cosmologies (i.e. ∆z/z) reaches the 12% level and is sufficient to solve the Hubble tension ( Figure 16). Observations at higher z could distinguish our void model from other possible solutions, e.g. miscalibrated SNe, early dark energy, etc. This is because a local void predicts that the inferred H 0 decreases with the redshift of the data set used, and asymptotically approaches the Planck prediction ( Figure 7).
Unfortunately, Figure 16 shows that high-redshift SNe currently do not pose strong constraints on our model. This is because the simulated ∆z/z decreases with redshift, while the uncertainty of binned SNe distance measurements to fixed z increases for z > ∼ 0.2. It crosses the simulated ∆z/z curve at z ≈ 0.38, corresponding to a lookback time of ≈ 4.2 Gyr. At z ≈ 0.5 (the dotted vertical line in Figure 16), the predicted ∆z/z ≈ 3%, but the observational uncertainty is much larger (see also Cuceu et al. 2019;Macaulay et al. 2019). As a result, even the 9% Hubble tension could not be reliably detected in SNe at these redshifts, which is reasonable given that the uncertainty of H 0 derived from SNe at 0.023 ≤ z ≤ 0.15 is already ≈ 2% (Camarena & Marra 2020a). Moreover, there are much fewer observed SNe at high redshifts, which could increase systematic errors.
In contrast to high-redshift SNe, the current BAO precision (Alam et al. 2017) lies slightly below the predicted ∆z/z. However, the uncertainties are still too large to distinguish our model from a void-free Planck cosmology at high significance. We note that BAO-based H 0 measurements (Alam et al. 2017;Zhang et al. 2019) are very close to the Planck prediction (Section 1.2), which is consistent with our void model. The small excess it predicts can only be confirmed or ruled out with more precise observations.
In conclusion, it is currently difficult to distinguish void and void-free models with data only at z > ∼ 0.5. Data at lower z are more useful in this regard. In particular, the redshift range 0.023 ≤ z ≤ 0.15 covered by Camarena & Marra (2020a) brackets the peak of the simulated ∆z/z curve and poses therewith a strong test of our model. In this redshift range, our void model differs from ΛCDM by ∆z/z = 7−12%, which is quite consistent with the 9.5% difference between local and early universe measures of H 0 . Note that the position of the peak in ∆z/z depends on the underlying void profile − it would occur at the void centre for a Gaussian or an exponential profile.

CMB contamination by intergalactic dust
Finally, we consider the possibility that the CMB is contaminated by intergalactic dust, which in turn would affect the H global 0 required by Planck Collaboration VI (2020). Some distant foreground emission can increase H global 0 (Yershov et al. 2020), which would slightly reduce the mild tension between our model and the strong lensing data (Figure 7). This is because the z > 0.4 lenses all give H 0 systematically above the Planck prediction by a similar extent. But changing the Planck H 0 would not explain the high inferred H 0 from the two lowest z lenses, which would continue to hint at a local void.
It is also possible that the CMB is more substantially affected by dust. In contrast to the ΛCDM model in which the CMB is explained as relic radiation from the early Universe at z ≈ 1100 (e.g. Bennett et al. 2003; Planck Collaboration VI 2020), an alternative model is that the entire CMB is thermal radiation from intergalactic dust particles heated up by starlight (Vavryčuk 2018). Assuming the observed intergalactic dust is in thermal equilibrium with the radiation field from galaxies, the model of Vavryčuk (2018) implies a dust temperature of T D = 2.776 K. This is only slightly higher than the measured T CMB = 2.72548±0.00057 K (Fixsen 2009). The exact value of T D depends on the amount of intergalactic dust and the intergalactic opacity ratios, which are both poorly known observationally. In future, it would be important to study which dust parameters are necessary to match the observed T CMB . In other words, it would be important to quantify the uncertainty on T D , which was not explicitly addressed by Vavryčuk (2018). Furthermore, it needs to be demonstrated that the model can yield the observed perfect black body spectrum within rather small uncertainties (Planck Collaboration XXVII 2014), and also yield nearly Gaussian temperature fluctuations (Planck Collaboration XXIII 2014). It is therefore possible that intergalactic dust significantly affects the CMB, but a detailed consideration of such a scenario is beyond the scope of our work. Interestingly, Vavryčuk (2019) showed in a subsequent study that the anomalous dimming of SNe Ia can in principle be explained by light extinction due to intergalactic dust.
In addition to heating by starlight, dust grains would also be heated by the primordial CMB, especially at high z. This may have caused rethermalization of the CMB by dust from the first stars at z ≈ 15. In this scenario, the angular scale evident in BAO measurements corresponds to a different co-moving length than the sound horizon at the time of last scattering. However, agreement can be recovered if we assume a non-standard background cosmology where a ∝ t (Melia 2020). In fact, it is not possible that a ∝ t without such a late rethermalization of the primordial CMB (Fujii 2020). While the late-time expansion history is indeed approximately of this form (Figure 12), the model does not yet explain the nature of the acoustic oscillations in the CMB power spectrum. Moreover, BBN would be modified to a very substantial extent, making it difficult to explain the observed light element abundances (Lewis et al. 2016).

CONCLUSIONS
Cosmic structure − and therewith the distribution of galaxies − provide strong constraints on the underlying cosmological model. In this context, we used the framework of the standard ΛCDM theory and MOND (Milgrom 1983) to investigate the KBC void, a large underdensity with a relative density contrast of δ ≡ 1− ρ/ρ 0 = 0.46±0.06 between z = 0.01 and z = 0.07 (Keenan et al. 2013). A large local underdensity is evident throughout the whole electromagnetic spectrum (Section 1.1).
Using the MXXL simulation (Angulo et al. 2012), we showed that the KBC void is in 6.04σ tension with standard cosmology (Section 2.2.1). In principle, if mass conservation is assumed, such an immense void should also show up in the velocity field, and would approximately solve the Hubble tension (Equation 5). This tension nowadays exceeds the 5σ threshold based on numerous independent techniques (Section 1.2). However, we demonstrated that a 10σ density fluctuation would be necessary to solve the Hubble tension respectively. The model yields a local Hubble constant of H model 0 = 76.15 km s −1 Mpc −1 and an acceleration parameter of q model 0 = 1.07 in the redshift range 0.023 ≤ z ≤ 0.15, consistent with the observations of Camarena & Marra (2020a) at the 84.20% confidence level (0.20σ tension, Section 4.2). Similar results are obtained for models initialized with a Gaussian or an exponential profile (Appendix B). Thus, we have shown for the first time that the Hubble tension can be solved in MOND. Several other tensions are also simultaneously resolved, notably the KBC void and that in q 0 (see also Camarena & Marra 2020b).
While all our best-fitting models generally imply larger peculiar velocities than the observed v LG of only 627 km s −1 , the possibility that v tot ≤ v LG cannot be excluded at the 99% confidence level (Section 4.2 and Appendix C). Thus, we do not require the LG to be at a special position within the KBC void in a statistically significant sense. Our results indicate that we should be 150 − 270 Mpc from the void centre in roughly the opposite direction to the external field on the void ( Figure 8).
As we go beyond the void, all our models predict that the inferred H 0 decreases with redshift. Indeed, observations of strongly lensed quasars taken by the H0LICOW team have shown that H 0 decreases with the lens redshift at a significance level of 1.9σ (Wong et al. 2020). However, our best-fitting model systematically underestimates the lensing-inferred H 0 , which may be related to the EFE sourced by a massive object beyond the void at z > ∼ 0.15 but closer than the lenses at z > 0.3 (Section 5.2.1). It could also be a sign of systematic errors (Kochanek 2020), but the sharp rise for the two nearest lenses is suggestive of a void-induced effect.
Taking into account all observational constraints (Section 3.3), our fiducial MOND model explains these local observations at the 1.14% confidence level, representing 2.53σ tension (Section 4.2). The best-fitting MOND models with a Gaussian and an exponential void profile are consistent with observations at 0.45% (2.84σ) and 0.34% (2.93σ), respectively (Appendix C).
Although strong lensing does not occur in the MOND regime (Sanders 1999) and works similarly to standard cosmology (Section 3.3.5), we also redo our analysis without the H 0 constraints from this method. Our best-fitting model is then consistent with observations at the 5.0% (1.96σ) confidence level, with only small changes to the best-fitting parameters (Section 5.2.1).
Our analysis strongly disfavours models without an EFE, consistent with results from wide binaries (Pittordis & Sutherland 2019). Furthermore, we showed that allowing time variation of the EFE has only a minor impact on our results because a constant EFE is well within uncertainties (Section 5.2.2). The main effect of allowing a stronger EFE in the past is to raise the required void strength at z = 9, with values up to ≈ 10 −3 becoming allowed at 1σ (Figure 14). This is more in line with the expected cosmic variance at that epoch.
We also discussed structure formation and the implications for the KBC void in MOND if peculiar accelerations are coupled to the Hubble flow acceleration g Hubble , as proposed by Sanders (2001). Such a coupling (or HFE) would effectively add g Hubble as an extra source of gravity when calculating the MOND boost to gravity, making the behaviour more nearly Newtonian (Sections 3.1.4 and 5.2.3). However, even a strong HFE implies a significant enhancement to gravity and the formation of voids compared to the Newtonian case. This is because g Hubble ≈ 0.2 a 0 on a 300 Mpc scale, and completely vanished 6 Gyr ago ( Figure 15). As a result, we conservatively estimated that even with a strong HFE, the cosmic variance in MOND would still be at least 2.8× that in standard ΛCDM on a 300 Mpc scale (≈ 9.0% instead of 3.2%). This would mean that whereas ΛCDM needs a 10σ density fluctuation to simultaneously explain the KBC void and Hubble tension (Figure 2), a MOND cosmology would only need an ≈ 2σ fluctuation (Section 5.2.3). Thus, MOND can successfully describe the density and velocity field on a Gpc scale under a wide range of plausible theoretical assumptions on how density perturbations couple to the background cosmology. In principle, the strength of the coupling introduces additional degrees of freedom that could be used to match the observed frequency of KBC-like voids, the observed lensing of the CMB, and the ISW effect. However, it is not clear if a covariant version of MOND has this flexibility when other constraints are imposed, e.g. that gravitational waves should travel at c. These theoretical uncertainties should be addressed in future work.
While the MONDian framework provides a reasonable fit to the locally observed density and velocity field, we emphasize that other alternative cosmologies might do so as well. Our results suggest that a successful model should have an expansion history similar to ΛCDM, but yield significantly more cosmic variance on a 300 Mpc scale. Additionally, the model must also accurately describe the dynamics of galaxies in order to provide a holistic explanation of the observed Universe. In this regard, a modification to gravity at length-scales beyond e.g. 10 Mpc would not be sufficient as it would face the same issues as ΛCDM on galaxy scales.
There are still considerable theoretical uncertainties in the here developed cosmological MOND simulation (Sections 3.1 and 5.2.3) because we lack an understanding of the fundamental theory behind MOND (i.e. FUNDA-MOND, Milgrom 2020a,b). Nevertheless, a promising relativistic MOND version was recently developed in which gravitational waves travel at the speed of light (Skordis & Z lośnik 2019). Its implications for cosmology should be explored, though a rather large box size would be required to reach the scale at which the CP holds in a Milgromian universe. This is because in MOND the EFE suppresses the growth of structure, causing structure formation in different regions to become correlated (Section 5.2.2). Without such simulations and/or further analytic work, we cannot draw any strong conclusions on the expected time evolution of the EFE. We nonetheless expect our results to hold because a wide range of possible EFE histories yield reasonable results, and because other void parameters such as its initial size and strength could be adjusted to optimize the fit (Figure 14). Any viable cosmological model has to explain both the local and global Universe. The KBC void is virtually impossible within the ΛCDM framework (Section 2.2). Consequently, the ΛCDM model faces serious challenges on Gpc scales, as shown in this contribution − the KBC void and Hubble tension falsify the ΛCDM paradigm at the 7.09σ level, and point towards much more rapid growth of structure than predicted by standard cosmology. Moreover, Di Valentino et al. (2020a) reported "a possible crisis for cosmology" based on the Planck power spectra, while Di Valentino et al. (2020c) concluded that the ΛCDM paradigm has to be replaced. These large-scale issues should be addressed together with the severe problems faced by ΛCDM on galactic scales (e.g. the satellite planes and the RAR, see also Kroupa 2015, and references therein).
Previous studies have shown that MOND is successful on several astrophysical scales ranging from the equilibrium dynamics of galaxies (Famaey & McGaugh 2012) and their formation out of gas clouds (Wittenburg et al. 2020), to the equilibrium dynamics of virialized galaxy clusters , and the formation of extreme clusters like El Gordo (e.g. Katz et al. 2013). The cluster-scale successes require the assumption of sterile neutrinos as HDM, which allows MOND to produce a standard expansion history and have very little effect on BBN and the high-acceleration CMB (Section 3.1). Consequently, there exist only very few (if any) scales at which the ΛCDM framework provides a unique explanation for the observations. Rather, observations of the local and global Universe strongly suggest that we should replace ΛCDM with the νHDM framework, which relies on MOND and sterile neutrinos.
The encouraging results we obtained using this approach should be put on a more secure theoretical footing using a covariant framework such as that of Skordis & Z lośnik (2019). In particular, it is important to rigorously demonstrate that the background cosmology behaves like in ΛCDM at the sub-per cent level. A covariant framework would also clarify if there is any coupling between the Hubble flow acceleration and that sourced by inhomgeneities. If there is and if its strength is adjustable, the value could be found empirically using numerical simulations of large-scale structure. Calculating photon propagation through the resulting time-varying inhomogenous gravitational field would then allow comparison with the observed lensing of the CMB by intervening structures, and the resulting ISW effect (Buchert 2000;Wiltshire 2007). Although these both appear to be underestimated in the ΛCDM framework (Section 5.3.1), they may be overestimated in νHDM. In this context, it is worth mentioning that the CMB Cold Spot could be caused by a KBC-like void (Nadathur et al. 2014). The expected frequency of such voids should be quantified using numerical simulations, which would also account for more complicated effects such as non-sphericity of the void. This may lead to predictions for angular dependence of the apparent expansion rate, which could be contrasted with observations (e.g. those of Migkas et al. 2020).
We conclude that unlike ΛCDM as presently understood, MOND supplemented by HDM appears to be a promising way to explain observations across all astrophysical scales. In particular, we expect this νHDM model to yield an almost standard expansion history but with enhanced cosmic variance on a 300 Mpc scale, allowing it to explain the observed KBC void and therewith the Hubble tension. This scenario has to be investigated in an open-minded manner in future studies.

DATA AVAILABILITY
The data underlying this article are available in the article.
indicating that the MXXL density fluctuations closely follow a normal distribution.

APPENDIX B: KBC VOID MASS PROFILES
In addition to our fiducial MOND simulation based on a Maxwell-Boltzmann void density profile (Section 3.2.1), we also model the void with a Gaussian and an exponential profile. The enclosed mass of the void within co-moving radius r com for a Gaussian profile is As before, x ≡ r com /r void , α void is the initial void strength, and r void is the initial co-moving void size at z = 9. The corresponding result for an exponential profile is In both cases, α void is the initial underdensity at the void centre. The results of using these void profiles are presented and compared with local observations in Appendix C.

APPENDIX C: RESULTS FOR DIFFERENT VOID PROFILES
The marginalized posterior distribution of the model parameters based on 10 6 MOND models for a Gaussian and an exponential initial void profile are shown in Figures C1 and C2, respectively. All these models assume a time-independent EFE (i.e. n EFE = 0 in Equation 50). As with the Maxwell-Boltzmann profile, models with a very weak or a very strong EFE are ruled out, but the initial void parameters are only weakly constrained by local observations. In particular, models with a Gaussian and an exponential profile restrict g ext to the range (0.045 − 0.127) a 0 and (0.045 − 0.117) a 0 at the 3σ level, respectively.
The best-fitting model for a Gaussian void profile has an external field strength of g ext = 0.070 a 0 , an initial void size of r void = 1030.0 cMpc (the upper limit of the allowed parameter range), and an initial void strength of α void = 3.76 × 10 −5 . This model is in 2.84σ (0.45%) tension with local observations (Section 3.3).
The results for both models are listed and compared with observations in Table C1. A time-dependent EFE and its implications for structure formation are studied in Section 5.2.2 for all three considered profiles. This paper has been typeset from a T E X/L A T E X file prepared by the authors.   ). The red dashed, black solid, and black dashed lines mark the 1σ, 2σ, and 3σ confidence levels, respectively. For 1D posteriors, these are shown using horizontal black lines. The red dot or vertical line marks the best-fitting model with an external field strength of g ext = 0.070 a 0 , an initial void size of r void = 1030 cMpc (the upper limit of the allowed parameter range), and an initial void strength of α void = 3.76 × 10 −5 at z = 9. ). The red dashed, black solid, and black dashed lines mark the 1σ, 2σ, and 3σ confidence levels, respectively. For 1D posteriors, these are shown using horizontal black lines. The red dot or vertical line marks the best-fitting model with an external field strength of g ext = 0.080 a 0 , an initial void size of r void = 1030 cMpc (the upper limit of the allowed range), and an initial void strength of α void = 7.56 × 10 −5 at z = 9.