A Bayesian approach to the cosmic dipole in radio galaxy surveys: Joint analysis of NVSS & RACS

We examine the sky distribution of radio galaxies in the NRAO VLA Sky Survey (NVSS) and the Rapid ASKAP Continuum Survey (RACS). Analyses of these samples have reported tension between their inferred dipoles and the kinematic dipole of the Cosmic Microwave Background (CMB). This represents a challenge to the traditional assumption that the Universe is homogeneous and isotropic on large scales: the cosmological principle. We find that NVSS and RACS contain local radio sources which give a non-negligible contribution to the overall dipole signal. These need to be adequately accounted for since the aim is to probe the composition of the Universe at large scales. By appropriately considering these sources, the inferred dipole amplitude in either sample is reduced. Nonetheless, we find support for a dipole aligning with that of the CMB but larger in amplitude, especially in the joint analysis. However, the ‘clustering dipole’ – the contribution of local sources to the net inferred dipole – appears to align with the direction of the CMB dipole, and its magnitude increases as deeper nearby sources are considered up to a comoving distance of ≈ 130 Mpc ( ℎ = 0 . 7). The significance of this observation in the context of the cosmological principle is unclear and prompts further inquiry.


INTRODUCTION
A critical assumption in modern cosmology -the Lambda cold dark matter (ΛCDM) paradigm -is that the Universe is homogeneous and isotropic on large scales.Such a postulate about the symmetry of space was initially made by Einstein, but it was later termed the 'cosmological principle' (CP) by Milne (Milne 1935;Peebles 2020).The CP is implicitly incorporated into the Friedmann-Lemaître-Robertson-Walker (FLRW) spacetime metric, greatly simplifying its form.Now, the CP can be supported by the smoothness of the Cosmic Microwave Background (CMB), although one must additionally couple this empirical observation with the metaphysical assumption that all observers in the Universe perceive the same CMB.
Small-scale CMB temperature fluctuations, generated from matter-photon interactions before they decoupled, are of order Δ/ ≈ 10 −5 .Imprinted on this smooth temperature map is a dipole of order Δ/ ≈ 10 −3 (the 'kinematic dipole'), which is normally interpreted as arising from the Earth's peculiar motion towards (, ) = (264.• 021, 48.• 253) in Galactic coordinates at a speed of  CMB = 369.82± 0.11 km s −1 (Planck Collaboration et al. 2020).But since the small-scale anisotropies are thought to be the progenitors of structures in the late universe under ΛCDM, the dipolesubtracted frame must be one in which the cosmological principle holds.In this frame, distributions of matter in the Universe should not have a preferred direction of higher or lower density.However, where the Earth's peculiar motion has not been factored out, distributions of ★ E-mail: oliver.oayda@sydney.edu.audistant matter should exhibit a dipole (the 'matter dipole') aligning in direction and magnitude with the CMB temperature dipole.If this is not the case, this suggests against the applicability of the CP to observational data, and hence to our prevailing understanding of the Universe more generally.
In recent studies, especially those using large catalogues of radio galaxies covering significant fractions of the sky, the matter dipole and CMB dipole appear to align in direction but not magnitude, with the matter dipole being in excess (see e.g. the reviews of Peebles 2022; Kumar Aluri et al. 2023).In this work, we turn to two radio catalogues that have previously provided evidence for this magnitude discrepancy, namely the Rapid Australian SKA Pathfinder Continuum Survey (RACS; McConnell et al. 2020) at 887.5 MHz (RACS-low; Hale et al. 2021) and the National Radio Astronomy Observatory Very Large Array Sky Survey (NVSS; Condon et al. 1998).We examine both catalogues with a joint Bayesian methodology, probing the anisotropic distribution of their radio galaxies and comparing it to the CMB dipole.This paper is structured as follows.In Section 2, we give an overview of the state of the literature relating to tests of the CP and the matter dipole studies.In Section 3, we present the raw NVSS and RACS-low samples and our approach to preparing them for analysis, and in Section 4, we explain our method of statistical analysis.We give our results in Section 5 and discuss them in Section 6, in which we conclude with our salient finding.

BACKGROUND
The CP forms the backbone of our cosmological framework, hence its empirical verification is extremely important.If the assumption of isotropy and homogeneity is correct, then the thermal dipole of the CMB is explained by the Earth's peculiar motion, which imprints a Doppler modulation on an otherwise isotropic CMB.An independent way of verifying this claim is to search for a dipolar modulation in surveys of distant sources such as radio galaxies and quasars.The Earth's motion is anticipated to induce this dipole, and it should be consistent with the CMB dipole in direction and magnitude if the cosmological principle holds.In this vein, Ellis & Baldwin (1984) proposed a number count dipole test of the CP.They assumed the following: (i) An observer is moving with a velocity much less than  with respect to a homogeneous and isotropic distribution of distant sources.
(ii) The apparent flux density  can be described by a cumulative power law, where  (> ) ∝  −  .
(iii) Within the observer's passband, the spectral energy distribution for sources has a power law dependence on frequency   ∝  −  for spectral index  and frequency  .
Under these assumptions, relativistic aberration and Doppler boosting will modulate the distribution of sources in the observer's frame, inducing a dipole anisotropy with amplitude Thus, the modulation in number density observed for a patch of sky containing  sources and pointing in some direction n will be Δ (2) Ellis & Baldwin (1984) supposed that  (10 5 ) sources would be needed to observe this modulation.It is worth mentioning that some of the assumptions made by the original authors have been questioned.For instance, Dalang & Bonvin (2022) (see also Guandalin et al. 2023) have claimed that the flux distribution exponent  as well as the spectral index  are redshift-dependent, which impacts the resulting dipole measure (although von Hausegger 2024 has suggested the dipole measure is robust to these effects).One other outstanding issue is the assumption that the inferred dipole is solely contributed to by the Earth's peculiar motion, which is a central theme in this work.In general, the inhomogeneous distribution of sources in the Earth's local neighbourhood can contribute to the overall dipole signal.This is an undesired effect, since the aim of the Ellis & Baldwin (1984) test is to probe the distant background of sources where the CP is expected to hold -that is, where the Universe is expected to average out to homogeneity and isotropy.For example, Tiwari & Nusser (2016) suggest that the bulk of sources analysed need to be at a high redshift ( ≈ 1) to avoid local source contamination.We will return to this point below.
The Ellis & Baldwin (1984) test has to date been widely performed with a number of different samples.One significant first result was that of Blake & Wall (2002b), in which the authors examined NVSS and reported broad agreement (within uncertainties) of the NVSS radio galaxy dipole and the kinematic dipole.However, this stands as an isolated result; many subsequent studies found a discrepancy between the matter dipole and the CMB dipole.We outline a handful of the key results here.Singal (2011) reported that the NVSS dipole is consistent in direction with the CMB dipole, but about 4 times larger than expected in magnitude.The trend of an unexpectedly large NVSS amplitude continued in subsequent works (see e.g.Gibelyou & Huterer 2012;Rubart & Schwarz 2013;Colin et al. 2017;Bengaly et al. 2018;Siewert et al. 2021;Singal 2023a,b;Wagenveld et al. 2023).A number of these studies used additional radio catalogues beyond NVSS (such as RACS for the later works), usually arriving at similar conclusions.The weight of these studies has lead to the prevailing sentiment that the matter dipole is consistent with the CMB dipole in direction but anomalously larger in magnitude (Abdalla et al. 2022;Kumar Aluri et al. 2023).
More recently, Secrest et al. (2021) expanded the scope of the matter dipole studies from radio sources to quasars.Analysing the CatWISE2020 catalogue (Marocco et al. 2021), the authors reported the presence of an anomalously large matter dipole at a 4.9 statistical significance.This finding was followed up in Secrest et al. (2022) with a joint analysis of CatWISE2020 and NVSS, confirming an excessive dipole amplitude at a 5.1 joint statistical significance.However, Abghari et al. (2024) has recently questioned the CatWISE2020 result, pointing to peculiarities regarding the ecliptic bias in the sample, as well as coupling between the dipole mode and higher order multipoles -which impacts the dipole estimator used in Secrest et al. (2021).
We note that, in general, the radio galaxy and quasar studies rely on frequentist statistical techniques, using estimators for calculating the dipole.Estimators are sometimes characterised by an inherent bias.As an example, Siewert et al. (2021) illustrate that the linear estimator used in earlier studies has a directional bias depending on which part of the sky is masked.In simulations, this has been shown to draw the inferred dipole away from the direction of the true dipole.In contrast to frequentist methods, the tools of Bayesian statistics have not been used as extensively.Bayesian techniques were implemented in Dam et al. (2023)'s study of CatWISE2020, which supported the findings of Secrest et al. (2021). However, Mittal et al. (2024) found that an independent catalogue -the Quaia quasar sample (Storey-Fisher et al. 2023) -was consistent with the kinematic dipole, although it was noted that the sample suffers from contamination near the Galactic plane which needs to be accounted for.
Apart from the matter dipole studies, many other methods to check for isotropy have been proposed.Kumar Aluri et al. (2023) gives an extensive discussion on these tests, but some of them involve probing Type Ia SNe (see e.g.Horstmann et al. 2022;Singal 2022;Sorrenti et al. 2022) and the FLRW metric by testing spatial curvature (see e.g.Zhou & Li 2020).Recently Oayda & Lewis (2023) proposed a method of calculating the dipole in the time dilation of sources with intrinsic time-scales.
The final point to consider is the effect of local clustering, which we alluded to earlier.Some studies make the assumption that, on account of most sources being at high redshifts, its effects will be negligible (see e.g.Wagenveld et al. 2023 for a recent contribution analysing NVSS and RACS-low).Other studies have anticipated the degree to which clustering will affect the overall dipole signal, using theoretical arguments from ΛCDM to construct a prior likelihood for a 'clustering dipole' term (see e.g.Tiwari & Nusser 2016;Cheng et al. 2023;Dam et al. 2023).Since the matter dipole amplitude is consistently in excess across many independent studies, it is worth probing the fraction of the amplitude that can be explained by local clustering, if at all.This is especially the case given the significance that a genuine tension would have for the CP and modern cosmology.Thus, we revisit two of the catalogues which have seemingly contributed to the tension: NVSS and RACS-low.Although many studies have calculated the matter dipole across a spectrum of different samples, their analysis is generally limited to characterising the dipole amplitude and direction and evaluating its significance with respect to a null hypothesis of no dipole signature.In this work, we jointly analyse NVSS and RACS-low, utilising Bayesian inference to compare the relative strengths of competing hypotheses by their marginal likelihoods and to understand if the inferred matter dipole is in agreement with the CMB kinematic dipole.

RADIO GALAXY SAMPLES
At the first stage of processing, both the NVSS and RACS-low source catalogues are binned into equal-area sky pixels using the healpix 1 procedure implemented in the python package healpy (Górski et al. 2005;Zonca et al. 2019).We set the  side parameter to 64, generating 49 152 pixels with an angular size of about 55 arcminutes.

Flux limit and mask
The NVSS survey was conducted using the Very Large Array (VLA) at 1.4 GHz, between 1993to 1997(Condon et al. 1998).It covers about 82% of the sky, above declination  ≥ −40 • ).The full source catalogue consists of 1.8 million sources.In this work, we make use of the integrated flux densities recorded for each source.
NVSS used two configurations of the VLA for different sky regions.Namely, the D configuration covered declinations −10 • ≤  ≤ 78 • , with the DnC configuration covering the remaining portion of the sky.Critically, this means that without a sufficiently bright flux density cut, NVSS shows a strong systematic bias in number density with declination -even though NVSS is claimed to reach 100% completeness by about 4 mJy (Condon et al. 1998).While some authors have chosen to mitigate this with a 10 mJy flux density cut (see e.g.Secrest et al. 2022), we instead select a 15 mJy limit in line with for example Tiwari & Nusser (2016) and Wagenveld et al. (2023).We illustrate the effect of declination on source density in Fig. 1.Moving from lower to higher flux density cuts (top to bottom), the average source density by declination becomes increasingly clustered around the mean density.From 2.5 mJy to 15 mJy, the average percentage 1 https://healpix.sourceforge.io/deviation drops from ≈ 1.8% to ≈ 0.9%.We also show the resulting flux density distribution with the 15 mJy cut applied in Fig. 2. Note that both Fig. 1 and Fig. 2 are generated after incorporating masking and source selection, described below.
NVSS suffers from significant contamination near the Galactic plane, appearing as localised regions of high source density.We therefore choose to mask Galactic latitudes || ≤ 10 • .However, we also identified other regions of spuriously high source counts.Accordingly, we follow some of the masking choices identified in Cheng et al. (2023).First, bright, extended sources can appear as multiple entries in the NVSS catalogue, so we impose an upper flux density limit of 1 Jy, such that the flux density distribution of our NVSS sample spans 15 mJy ≤  ≤ 1000 mJy.As an additional measure, one might choose to mask all pixels containing a bright source with flux greater than 1 Jy, or alternatively construct a circular mask centred around that bright source (see e.g.Blake & Wall 2002a).We tested this, but ultimately found it had not effect on our results.Therefore, we persist with the upper flux density cutoff.We then visually identify a number of localised regions with substantially higher number counts and mask these regions with discs of varying radii.We show the Galactic coordinates of these regions in the corresponding rows of Table 1.A number of known bright radio sources are located in or at the edge of these masked discs, which we specify in the 'A-source' column of the same table.Where there is no relevant nearby source, the cell is left blank.Lastly, we mask an additional degree north of the NVSS declination cutoff to remove pixels with lower source counts at the survey limit.Our continuous region includes  > 40 • , leaving a final count of 341,072 NVSS sources.The main features of the mask can be seen in the left panes of Fig. 5.

Local clustering
It is conceivable that the inhomogeneous distribution of low redshift sources will contribute to the recovered dipole signal in NVSS.Since the cosmological principle is a statement about the composition of the Universe on the largest scales, ideally any contamination from nearby clusters should be avoided.As mentioned, some studies (see e.g.Dam et al. 2023;Cheng et al. 2023  into account predictions from ΛCDM.In this case, the net inferred dipole would include a contribution from a clustering term. Here, we investigate the effect of clustering by cross-matching NVSS with known local radio sources, creating two variant catalogues: one with local sources included, and one in which local sources have been removed.The first is constructed as we have already described above.To construct the second, similar to Cheng et al. (2023) we use the source catalogue from the Two Micron All Sky Survey Redshift Survey (2MRS; Huchra et al. 2012) and radio sources with redshift less than 0.01 identified on the NASA/IPAC Extragalactic Database (NED). 2 The positional accuracy of NVSS sources with flux density greater than 15 mJy is ≲ 1 arcsecond.We therefore adopt a cross-matching radius of 5 arcseconds to capture as many genuine local sources as possible while minimising spurious matches.Specifically, for each local source, we find all NVSS sources within 5 arcseconds and remove the one with the smallest positional offset.We apply the cross-matching regime before masking but after flux-limiting the sample.With the 15 mJy cut, 3049 sources are removed from the NVSS catalogue.After masking, this final sample has 338,222 NVSS sources, 2850 less than the sample with local sources included.For future reference, we refer to sample with local sources removed as NVSS(B), and the sample containing local sources as NVSS(A).

Flux limit and mask
The RACS-low survey (Hale et al. 2021) was conducted with ASKAP at a central frequency of 887.5 MHz.Observations took place between 2019 and 2020.The survey covered declinations  ≤ 41 • .The resolution of RACS-low changes with declination, and so in addition to the main catalogue, a catalogue constructed from data convolved to a common resolution of 25 arcseconds was produced.This is the catalogue we used in this work.However, not all observational tiles could be convolved to the common resolution, and so the effective coverage of the catalogue used here is a continuous region between −80 • ≤  ≤ 30 • , excluding the Galactic plane for || < 5 • .This covers about 67.9% of the sky, and the catalogue contains about 2.1 million sources.
We first probe the impact of our choice of limiting flux density.All  our references to flux with respect to RACS-low refer to the integrated flux density, as was the case for NVSS.RACS-low is asserted to be 95% complete at about 3 mJy, though we note improvements in the homogeneity of the sample with deeper flux cuts up to about 15 mJy.To see this, we plot how the source density changes with declination for each flux limit in Fig 3 .From 2.5 mJy to 15 mJy, the average percentage deviation from the mean density decreases from about 4.3% to 1.1%.We ultimately use a 15 mJy cut for our final sample, and, as for NVSS, use an upper flux density limit of 1000 mJy.We also show the sample flux distribution in Fig. 4. Note that both Fig. 3 and Fig. 4 are made after the masking and source selection that we describe below.
Turning to the choice of mask, we remove under-dense pixels along the edge of the excluded Galactic plane.We also mask the excluded survey tiles at the southern pole with a disc of radius 13 • , and mask an additional degree of declination in the northern hemisphere such that our final sample covers declinations −77 • <  < 29 • .Finally, we note localised regions with a dearth of source counts which persist in our sample, corresponding to under-dense tiles in the original RACS- low merged sky catalogue.These can be seen in Figure 4 of Hale et al. (2021) near the upper declination limit of the survey.Similar to our approach for NVSS, we mask these regions with discs of radius ≈ 3 • .The position, exact radius and any known nearby bright radio sources are specified in the corresponding rows of Table 1.We also mask an over-dense region proximate to the LMC, as shown in the same table.This leaves us with a sample of 462,911 RACS-low sources.The key features of the mask can be seen in the right panes of Fig. 5.

Local clustering
We repeat the same cross-matching procedure as for NVSS (see Section 3.1.2).The typical RACS-low source position is accurate to within about 1 to 2 arcseconds.For this reason, we keep the same cross-matching radius of 5 arcseconds as with NVSS.This identifies 3700 local objects in RACS-low between the 2MRS and NED catalogues, which we remove, leaving a final source count after masking of 459,276.We use the same naming convention as above, referring to the sample with local sources removed as RACS-low(B), and the variant containing local sources as RACS-low(A).

Summary
We give the sky projections of our final samples in the top row of Fig. 5.Note that there we have plotted the B variants of NVSS and RACS-low (the cross-matched samples), since with a visual inspection there are only very minor differences between the A and B variants.We also show a smoothed version of the raw samples in the bottom row of the same figure.Here, for each unmasked pixel, we compute the average density across all pixels within 1 steradian.We also give a summary of the key choices we made in processing both catalogues in Table 2 to aid the reader.

Expected amplitude
The conventional approach to finding the expected matter dipole amplitude is to use the spectral indices of the sources  and the slope of the flux density distribution  at the low flux density limit.These terms are then used in the expression of Ellis & Baldwin (1984) -see equation ( 1) -with  ≈ 369 km s −1 , as determined from the CMB temperature dipole.

Spectral index
Since both the RACS-low and NVSS catalogues list flux densities at a single frequency, we do not at first instance have information about each source's spectral index.Typically, the expected spectral index for synchrotron emission -the assumed emission mechanism -is used ( ≈ 0.75; see e.g.Secrest et al. 2022).Often, overlapping sources between radio surveys at different frequencies can be used to infer the spectral index distribution, the median of which is the origin of the typical value of  ≈ 0.75.For example, Figure 8 of Mauch et al. (2003) gives the distribution of spectral indices for NVSS sources overlapping with the Sydney University Molonglo Sky Survey (SUMSS), showing that the distribution is reasonably broad about the median of ≈ 0.83.Although we find that the expected amplitude is not very sensitive to the choice of , to safeguard against the possibility, we can assume that the median spectral index is 0.75 but has a significant spread.We take the largest dispersion from Figure 8, being about 0.5, and assume that the distribution of NVSS spectral indices can be described by a Gaussian G( = 0.75,  = 0.5).For RACS-low, Hale et al. (2021) found  distributions through comparison with four other surveys at different frequencies, including NVSS.
The typical  was between 0.6 and 0.9, with dispersions at most being ≈ 0.5, and we therefore use the same Gaussian distribution to describe the RACS-low spectral indices.

Source number counts
While we could use the slope of the flux density distribution  to ultimately arrive at the expected amplitude, we instead follow the approach developed in Mittal et al. ( 2024) but with modifications to allow flux density measurement uncertainties to be adequately accounted for.Namely, we randomly regenerate either catalogue a number of times using the recorded flux density and its associated uncertainty, as we describe below.(ii) For each source , we draw a spectral index   from the Gaussian G( = 0.75,  = 0.5).
(iii) We then determine the number of new sources   above some limiting flux density  0 .
(v) We then find the number of boosted sources   above the same  0 and multiply that number by  2 , accounting for relativistic aberration.
(vi) The expected amplitude is then determined via For a given  0 , we can repeat (i)-(vi) a number of times to arrive at a mean and standard deviation for the expected dipole amplitude.However, if we choose  0 to be too close to the actual flux density limit we impose for NVSS and RACS-low (15 mJy), then the dipole amplitude is underestimated.This can be seen by the drop-off near  0 = 15 mJy in Fig. 6.The drop-off in dipole amplitude is a boundary effect caused by the manner in which we resample flux densities.
To see this, suppose we pick  0 = 25 mJy.A source flux density slightly below this limit has a non-negligible likelihood of being resampled above the limit, owing to its measurement uncertainty.This is averaged over the many runs of our procedure.But, if we choose  0 ≈ 15 mJy, sources below 15 mJy cannot be resampled above the limit since, by construction, they are excluded from either catalogue.Meanwhile, the dipole amplitude is expected to rise with  0 , as seen in Fig. 6, since the NVSS and RACS-low cumulative flux density distributions become steeper for larger limiting flux densities.
We can circumvent this issue by extrapolating the gradient in   D CMB versus  0 back to 15 mJy, as seen by the dashed line in Fig. 6.We do so by running our procedure (steps i-vi) 100 times for different values of  0 between 15 mJy and 35 mJy, computing the mean dipole amplitude over those runs.We then determine the gradient using a straightforward least-squares linear fit for  0 > 17 mJy (NVSS) and  0 > 20 mJy (RACS-low).Extrapolating back to  0 = 15 mJy, we arrive at D CMB = 4.31 × 10 −3 for NVSS and D CMB = 4.27 × 10 −3 for RACS-low.Note that in the procedure described above, the B/cross-matched sample variant is used; the change in amplitude using the A variant is negligible.

Overview
Our statistical approach is similar to and explained in depth in Mittal et al. ( 2024), but we recapitulate the main concepts here.The process of Bayesian inference can be discretised into two stages: the level of parameter optimisation and the level of model inference (Mackay 2003).In the first stage, we wish to populate the terms of Bayes's theorem and solve for the posterior distribution of a model 's parameters, where for data D and model parameters , as well as likelihood, prior and marginal likelihood functions L,  and Z respectively (as using the notation of Speagle 2020).At the second stage, models are ranked according to the marginal likelihood, an integral over all parameter space: Z = ∫ Ω  L () × () .This works by reframing Bayes' theorem as the posterior probability for a model directly, such that Thus, noting that the term Z(D) cancels, the ratio of posterior probabilities for two models is where L (D|) is identical to the marginal likelihood Z of equation (4) and    is the Bayes factor: the ratio of model marginal likelihoods.() is identified as the model prior likelihood, representing our beliefs about the relative strengths of each model before any knowledge of the data (see Section 4.2.3).Accordingly, in this work we compare models using the natural logarithm of the Bayes factor, using the null hypothesis as a common benchmark.For example, if models  1 and  2 have marginal likelihood Z 1 and Z 2 , and the null hypothesis has marginal likelihood Z 0 , we can effectively rank the explanatory power of  1 and  2 by comparing ln  10 = Z 1 − Z 0 and ln  20 = Z 2 − Z 0 .This is advantageous where there are many models and it would be cumbersome to directly write down the Bayes factors between all possible pairs of models.The degree to which one Bayes factor is larger than the other is a reflection of the extent to which a model better explains the data over others while conserving prior volume.A qualitative interpretation (the 'Jeffrey's scale') of the Bayes factors can be found in Kass & Raftery (1995); we use those normative categories in this work, although we caution that they use units of 2 ln , whereas here all Bayes factors are presented as ln .
What remains is to fill the terms of Bayes's theorem at equation ( 4).For each model, we construct a likelihood function and determine the prior likelihoods for its parameters, which we explain below.Then, parameter optimisation, as well as the calculation of the marginal likelihood, is handled by the python package dynesty (Koposov et al. 2023a). 3dynesty implements the Nested Sampling algorithm to efficiently sample the posterior distribution in shells of increasing likelihood, also allowing the marginal likelihood to be evaluated (Skilling 2004(Skilling , 2006)).

Likelihood functions
In Mittal et al. (2024), two types of likelihood functions were tested alongside each other: the Poissonian approach and the point-by-point approach.Both were found to give consistent results, and so in this work we chose to use only the point-by-point case, mainly since it is computationally less expensive.Under this approach, we consider some arbitrary function  =  ( p ) which gives a scalar for the -th pixel on the sky as described by the unit vector p pointing to that pixel.Such a function needs to be normalised where for example arbitrary sky masks are chosen.Thus, the normalised value of the function at pixel  is for total number of pixels  pix .This is our sky probability map, the form of which depends on the assumptions about the chosen signal: a dipole, monopole, etc.The likelihood function can then be written as for number of sources   in pixel .This works because each source in a pixel is associated with the same function value f , and so the -th pixel contributes the value [ f ( p )]   to the likelihood function.
With this established, we need only to consider what the functional form of  will be for each model.
(i) Starting with the null hypothesis, we suppose that either sample exhibits uniform density over the sky.In this case, we need to weight all pixels equally, favouring no part of the sky over other regions.This can be done by setting  to 1 for all pixels, and so this model has zero parameters i.e.  null = ∅.
(ii) Next, we introduce a dipole vector D of magnitude D and pointing in some direction ,  in Galactic coordinates.The scalar value for an arbitrary pixel on the sky is the sum of the monopole signal and the dipole signal D • p .That is, we fit a monopole and dipole simultaneously.This leads to where   is the angle between the -th pixel unit vector and the dipole vector.Of course, the value of   depends on the direction of the dipole, which leaves the model parameters as  dipole = {D, , }.
(iii) We may also fix any of the parameters in the dipole model to constant values in order to test alternative hypotheses.Specifically, we fix (, ) = (264.• 021, 48.• 253) -the direction of the CMB dipole -and term this the 'kinematic direction' model.
(iv) We also fix the magnitude D to either of the expectations determined in Section 4.1 while leaving the direction free, terming this the 'kinematic velocity' model.
(v) We additionally investigate the 'kinematic dipole' model in which all parameters ,  and  are fixed to their CMB-derived values.
(vi) Finally, in our joint analysis (described below), for one of our models we set the prior likelihood function for the model parameters to the posterior distribution of Wagenveld et al. (2023).We approximate this distribution, as based on the information provided in the study, by drawing the parameters from Gaussians with mean and standard deviation given by the quoted median value and error.Namely, we have D NVSS = D RACS ∼ G( = 0.0129,  = 0.0018),  ∼ G( = 269.• 5,  = 8 • ) and  ∼ G( = 56.• 2,  = 11 • ).This is the 'W23' model.
We summarise the six models mentioned above in Table 3, in which we have labelled each model  0 through  5 .
In this work, we analyse NVSS and RACS-low individually as well as jointly, so we consider now how our approach differs in the joint case.Some studies, for example Darling (2022), combine multiple radio galaxy catalogues by scaling fluxes according to the source spectral index.There is the possibility that doing so leads to a spurious dipole signal by neglecting varying systematic effects across both catalogues, as was noted in Wagenveld et al. (2023).Here, we need not combine RACS-low and NVSS; rather, we suppose each catalogue is influenced by the same dipole signal, although, since they will not have the exact same flux or spectral index distribution, the dipole magnitudes will be different.That is to say, this joint model has parameters  = {D NVSS , D RACS , , }.We fit a dipole to either catalogue separately, but combine information across both catalogues with the joint likelihood function ln L = ln L NVSS + ln L RACS (10) where either individual likelihood is determined through equation (8).

Prior likelihood functions
Our choice of priors for the model parameters D, ,  is conditioned on the information we have before any knowledge of the data.We describe the choices for each in turn.
(i) For all dipole amplitudes, we sample from a uniform distribution U [0, 0.1].This reflects the spectrum of different results across independent studies (see e.g.Abdalla et al. 2022;Kumar Aluri et al. 2023).
(ii) We sample  and  uniformly over the sky, preferring no direction over another.To do this, since internally our calculations are performed in equatorial coordinates, we sample the right ascension in radians according to  ∼ U [0, 2].We then sample the codeclination according to  ∼ cos −1 (1 − 2) where  ∼ U [0, 1].It is important not to sample the co-declination straightforwardly from a uniform distribution between 0 and  since the area element of a sphere is dependent on the co-latitude; doing so would lead to a bias in the prior likelihood function towards the poles.
(iii) For the prior likelihood of each model (), we set each model to be equiprobable; that is, () = 1/5 where each sample is analysed individually and () = 1/6 in the joint analysis.

RESULTS
We present the results for each of our models, including Bayes factors and posterior distributions where noted.As a quick reference for the reader, we also summarise all inferred dipole amplitudes in Table 7.

NVSS
We give our Bayes factors by model for both NVSS(A) and NVSS(B) in Table 4.For NVSS(A), the model with the greatest explanatory power is the kinematic direction model, which again assumes a dipole aligning with the CMB dipole but free in magnitude.In this case, we find an amplitude of D ≈ (12 ± 5) × 10 −3 with 2 uncertainties.Although the free dipole is not the most favoured model, we represent its posterior distribution for the fit to NVSS(A) with a corner plot, shown in the left pane of Fig. 7 (top row).This is largely for illustration, indicating the dipole signature that is picked up at the first level of inference.Note that the inferred amplitude of the free dipole ( 1 ) is not necessarily the same as that of the kinematic direction model ( 3 ), since these are distinct hypotheses.We show the free dipole fits for all subsequent samples as well.
For NVSS(B), we instead find that the kinematic dipole model best explains the data, which assumes a dipole totally aligning with the CMB dipole in direction and magnitude.This is illustrated in the corresponding highlighted row of Table 4.However, the difference in marginal likelihoods between this model and the next best explanation -the kinematic direction model -is 0.3 log units.This is only worth a 'bare mention of support' for the kinematic dipole (Kass & Raftery 1995).Again, we show the corner plot for the free dipole fit on NVSS(B) in the right pane of Fig. 7 (top row).This illustrates that the free dipole model infers an amplitude of D ≈ (11 ± 6) × 10 −3 .Notably, the median dipole amplitude has diminished by ≈ 2 × 10 −3 moving from the NVSS(A) to NVSS(B) free dipole fit, although the uncertainties remain large.Lastly, if we fix the direction of the NVSS(B) dipole to that of the CMB ( 3 ), we find that D ≈ (10 ± 5) × 10 −3 .

RACS-low
We give the Bayes factors for both RACS-low samples in Table 5.For both RACS-low(A) and RACS-low(B), the prevailing model is a free dipole i.e. a dipole where D, , and  are free parameters ( 1 ).We give the posterior distributions for the free dipole -shown via corner plots -for both samples in the bottom row of Fig. 7. Similar to the case of NVSS, the free dipole amplitude for RACS-low decreases moving from the A to B variant by ≈ 2 × 10 −3 .For the A variant, we find a free dipole amplitude of D ≈ 14 +4 −5 × 10 −3 , and for the for the B variant, we find D ≈ (12 ± 5) × 10 −3 .
A dipole fixed to the kinematic direction but with free amplitude ( 3 ) has the second-highest level of support across both variants.In this case, we infer D ≈ 13 +4 −5 × 10 −3 for the A variant and D ≈ 11 +4 −5 × 10 −3 for the B variant.

Joint analysis
We give the Bayes factors for our joint analysis in Table 6.Here, the prevailing explanation is the W23 model ( 5 ), with ln  50 = 24.5 where both sample A variants are used and ln  50 = 17.1 where both sample B variants are used.Compared to the kinematic dipole model, the W23 model has a log Bayes factor of ln  54 = 10.5 for the A variants and ln  54 = 5.7 for the B variants.This is beyond overwhelming support for the model with prior likelihood given by the results of Wagenveld et al. (2023) over the kinematic dipole model, as using the scale of Kass & Raftery (1995).With this prior likelihood, we infer D NVSS = (12 ± 3) × 10 3 and D RACS = (13 ± 3) ×10 3 for the sample A variants, as well as D NVSS = (12±3) ×10 Free dipole: joint Lastly, we show the model  1 corner plot for the joint analysis with the B variants in Fig. 8, and project the marginal distributions for  and  onto the sky in Fig. 9. There, we find D NVSS = (9±6) ×10 −3 and D RACS = (12 ± 5) × 10 −3 .In the case of the joint analysis with the A variants, for the free dipole fit we find D NVSS = (11±6) ×10 −3 and D RACS = (13 ± 5) × 10 −3 .Table 7. Inferred dipole amplitude by sample and model.All values are reported with a 95% credible interval and given in units of 10 −3 .

Effect of local clustering
The results we presented in the previous section indicate that, while the dipole in RACS-low and NVSS is consistent in direction with the CMB dipole, there is significant support for an amplitude larger than that of the CMB dipole.While the case is not decisive for NVSS(B) alone, with the kinematic dipole being marginally preferred over the other hypotheses, RACS-low(B) has sizeable evidential support for a dipole larger than the CMB expectation (D CMB ≈ 4.0 × 10 −3 ).In addition, the joint analysis -either with the A or B variantsfavours a larger dipole amplitude, even if we neglect the results of Wagenveld et al. (2023) and consider model ( 3 ) to be the dominant explanation.
Nonetheless, one important observation is that the dipole amplitude in the individual catalogues diminishes moving from the A to B variants.Thus the clustering of low-redshift radio sources has a nonnegligible impact on the dipole amplitude.In fact, even though about 3000 sources were removed from NVSS and about 3700 sources from RACS-low, representing only 0.9% and 0.8% of the final source counts of each sample respectively, the inferred dipole amplitudes decreases by ≈10%-15% depending on the sample used.This accounts for why, in both B samples, the kinematic dipole model gains a comparitively larger marginal likelihood with respect to the other competing models (see Tables 4 and 5).

Clustering by redshift
We can further probe the effect of local clustering on the dipole amplitude by noting how the amplitude changes with increasingly deeper redshift cuts.The NED catalogue contains objects at a redshift of 0.01 and below, whereas 2MRS contains objects up to a redshift of about 0.1.Although the homogeneity scale -the distance over which the Universe is expected to average out to FLRW -is presently thought to be smaller than that implied by  ≈ 0.1 (see e.g.Sarkar et al. 2009;Scrimgeour et al. 2012 which place the transition to homogeneity around 100-115 Mpc ( ≈ 0.023-0.027)assuming ℎ = 0.7), there remains considerable uncertainty over its exact number (Kumar Aluri et al. 2023).It is therefore worth seeing how deeper cuts impact any clustering bias in the dipole amplitude.To do this, we simply repeat our model fitting for different redshift bins.In the first case, we remove no sources, and in the second case, we only remove sources with redshift ≤ 0.01.Then, we remove sources with  ≤ 0.02, and so on.
Our results for this process are presented in Fig. 10.Immediately, it can be inferred that the dipole amplitude continues to decrease up until  ≈ 0.04 (comoving distance of ≈ 170 Mpc), where it starts to plateau off.This is broadly the case whether or not one considers the samples individually or the joint analysis, in which both catalogues are assumed to be affected by a dipole pointing in the same direction but with different magnitudes.We note that -as mentioned earlierthe uncertainties on the amplitude are quite significant.None the less, by the  < 0.03 cut, the expected D NVSS (shown by the dashed brown line in the top pane of Fig. 10) is just outside the 2 uncertainties on the inferred dipole amplitude.However, the disagreement between the expectation and the inferred amplitude is much larger for RACSlow, reaching ≈ 3.4 even after local sources have been identified and removed from the sample.
To visualise the local structure being caught here, we plot the spatial distribution of NED and 2MRS sources which have a redshift lower than 0.02 in the top pane of Fig. 11.A significant band of low redshift structure concentrated near the CMB dipole can be seen.This band of structure is being picked up in the cross-matched sources identified in NVSS and RACS-low, shown in the middle and bottom panes of Fig. 11 respectively.This likely explains why the inferred dipole amplitude decreases until the redshift bin of  < 0.04; this local structure is being incrementally removed until the bin is sufficiently deep to cover its spatial volume.
With this established, it is worth mentioning the role that the NED sample has in our study.Our chief aim is to remove as many low-redshift sources as possible.However, since the NED database collates objects from independent studies and resources, there is the possibility that it introduces a bias or some other unknown effect in our results -especially since those studies may have different levels of completeness and coverage.As a consistency check, we repeated our above analysis without the NED sample and only the 2MRS sample, which, unlike the former, is constituted solely from the Two Micron All Sky Survey.We find that our results are essentially unchanged with this modification, likely because the cross-matched NED sources represent only ≈ 10% of the net cross-matched sources.We can therefore say that our findings are not influenced by an unforeseen bias in the NED objects.This supports our observation that clustering is a genuine influence on the dipole amplitude.

Clustering dipole
Another way of understanding the effect of local clustering is by visualising the dipolar contribution of local sources to the overall inferred dipole.We refer to this as the 'clustering dipole'.More explicitly, we take the cross-matched sources for either NVSS or RACS-low -as shown in the bottom two panes of Fig 11 -and fit a free dipole (model  1 ) to them.For the cross-matched NVSS sources, we infer a dipole amplitude of D ≈ 0.25 ± 0.06, and for the cross-matched RACS-low sources, we infer an amplitude of D ≈ 0.22 ± 0.05.We project the resulting marginal posterior distribution for  and  onto the sky in Fig. 12 after performing this fit.Strikingly, the cross-matched sources -both for NVSS and RACS-low -have a well-constrained dipole signal pointing near the CMB dipole.A priori we expect that the clustering dipole might point in any direction on the sky, as it ultimately traces the net dipole component of nearby structure.Meanwhile, the overall CMB dipole includes contributions from, for example, the motion of our heliocentric frame of reference with respect to the Galaxy.This is not necessarily correlated with local structure.
However, it is worth pointing out that the CMB dipole -i.e. the motion of the Sun with respect to the CMB -points in roughly the same direction as the motion of the Local Group (LG) with respect to the CMB (,  ≈ 270 • , 30 • ; see Table 3 in Planck Collaboration et al. 2020).Meanwhile, the motion of our heliocentric frame with respect to the LG is approximately oppositely-aligned with the CMB dipole.Accordingly, the fact that the clustering dipole seems to align with the CMB dipole means that it is also roughly aligned with the axis defining the LG's motion with respect to the CMB.It is therefore possible that the kinematic and clustering dipoles are correlated as implied in Fig. 12, since the cause of the LG's motion with respect to the CMB is gravitational interaction with local structure.Even so, it is curious that the clustering dipole seems to be well-defined with a small number of sources (< 4000) and aligns so closely with the CMB dipole.At least for this study, it explains why the removal of local sources is correlated with a decreasing dipole amplitude; the clustering and matter dipoles act in concert, contributing to a larger inferred amplitude.This is a 10%-15% effect, as was noted earlier.
The possibility of a local void or overdensity contributing to the NVSS net dipole amplitude was considered in, for example, Rubart et al. (2014).There, the authors simulated the effect of a perturbed region with source density below the average ( = −1/3) out to redshifts of 0.07, arriving at a similar value of D clust.≈ 2 × 10 −3 (although their approach is totally distinct from our own).In a similar vein, Bengaly et al. ( 2019) constructed a forward-looking sample, using predictions from ΛCDM to anticipate the radio number count map that will be made available by the upcoming Square Kilometre Array (SKA).The authors found that the dipole estimate in this mock sample is contaminated by local structure, with most of the inconsistency between the true kinematic dipole and the inferred dipole caused by sources out to  = 0.1.Thus, they expect crossmatching radio SKA sources with optical/infrared data to be of key importance in excising low redshift sources, which is what we also conclude in this work.In any case, the two studies lend support to our proposition that local structure contributes non-negligibly to the inferred dipole in NVSS and RACS-low.

Joint analysis
In our joint analysis, we give results both where the catalogue A variants and catalogue B variants are used (see Table 6).To reiterate, the B variants have had the local structure in Fig. 11 removed by crossmatching with both 2MRS and NED.The results with the B variants are representative of our main finding.In this case, we find that the model assuming the results of Wagenveld et al. (2023) as the prior likelihood function yields the highest evidence by a significant margin.If we blind ourselves to the findings of that work and consider the next-best explanation of the data, we find that fixing the dipole to the CMB direction but allowing the amplitude to vary as a free parameter offers the greatest predictive power.
We can understand this result by turning our attention to the marginal distributions of the dipole amplitudes, as shown in Fig. 13.There, one can see that by fixing  and  to the CMB dipole direction (top pane), a dipole with larger amplitude than the CMB expectation (brown line) is preferred.Models which fix D to the the CMB expectation are consequently penalised; they give a poorer explanation of the data.
One interesting point is that, if we assume a free dipole ( 1 ), the posterior distribution of the RACS-low amplitude tends to higher values than for NVSS.This is shown in the middle pane of Fig. 13, although the effect is also discernible for the other models (see the other panes).This is also present in our analysis of the individual samples, as shown in Fig. 7.It could therefore be asked why, on net, the RACS-low dipole tends to be larger than that NVSS -both in the joint and individual analyses.This cannot be explained by variation in the spectral index and/or flux distributions of each sample, since we convert to an implied velocity in Fig. 13 using the dipole expectations.
It is possible that this amplitude discrepancy evinces further systematic effects biasing the amplitude which we have not been able to adequately account for in RACS-low.For example, it is reasonable to suppose that we have not removed all possible local radio sources in our cross-matching procedure.This is especially the case considering the sensitivty of the amplitude to local clustering.Tiwari & Nusser (2016) estimate from mock NVSS catalogues that sources with redshift  < 0.1 contribute to the dipole amplitude by as much as 70%.With a different methodology, the fiducial model of Cheng et al. (2023) estimates the clustering dipole (the component of the inferred dipole arising from locally-clustered sources) in NVSS to be ≈ (4 ± 1) × 10 −3 , 80% of which comes from objects below a redshift of 0.05.This is about twice as large as the clustering component we inferred in the foregoing sections.Therefore, while the top pane of Figure 11 might be a good representation of local structure, it only provides a lower limit on the contribution to the clustering dipole; there may be more sources which have not been accounted for.

Comparison with other studies
Keeping the effects of clustering in mind, how does our result compare to other studies?Speaking broadly, the current state of the literature leans towards a matter dipole pointing in the same direction of the CMB dipole -which is consistent with our joint analysis -but with a magnitude in excess of the CMB dipole amplitude (Peebles 2022;Kumar Aluri et al. 2023).Now, it is important to recognise that our amplitude estimates have significant uncertainties.That being said, at face value our results are consistent with prior observations of an excessive dipole amplitude.For example, our results are evidently consistent with Wagenveld et al. (2023), the prior knowledge of which offers the best explanation of the data.Although, we note that the authors' approach differed slightly from ours, since they assumed a common dipole amplitude across NVSS and RACS-low, whereas we allowed it to vary between the samples.In doing so, we have made the additional observation that the RACS-low sample prefers a larger dipole amplitude.One other point of contrast is that the authors assumed that the effects of clustering are negligible; we have found substantial evidence that this is not the case.Though accounting for clustering only slightly reduces the dipole amplitude for the W23 model ( 5 ; D ≈ 0.013 → 0.012), it has a more pronounced effect on our other models.We interpret this as evidence for clustering playing an important effect in the Ellis & Baldwin (1984) test.
Turning to a separate study, Secrest et al. (2022) also examined NVSS, finding an amplitude of D ≈ (12 ± 3) × 10 −3 .This result is encompassed by the uncertainties on our findings, whether or not one turns to NVSS(A), NVSS(B) or the joint analysis.The fact that the authors did not account for the clustering of local sources can explain why their median amplitude is slightly larger than what we inferred after cross-matching (see NVSS(B) in Table 7).We cannot comment on the joint case, since the authors compared the quasar catalogue of CatWISE2020 (Marocco et al. 2021) and NVSS.
Additionally, in Colin et al. ( 2017), the authors accounted for local clustering by cross-matching with 2MRS.Nonetheless, this did not significantly impact their dipole amplitude estimate, which was around four times greater than the expectation.The lack of sensitivity of the dipole amplitude to local sources conflicts with our findings.Now, the authors cut and combined NVSS with SUMSS to create a merged catalogue with ≈ 600 000 sources, rather than considering both samples jointly with a likelihood function like equation ( 10).Thus, it may be worth revisiting their work in the future with a joint Bayesian approach and seeing if the contribution from SUMSS still leaves an anomalously high dipole amplitude -even when clustering is mitigated.

Future outlook
On the basis of the above, our result raises questions about the nature of the tension between the CMB dipole and the matter dipole.It sits in agreement with historic studies which have reported an anomalously large radio galaxy dipole, as well the result of Secrest et al. (2021), which confirmed this with quasars in the near-IR.However, one outstanding issue is why our finding is mildly in tension with that of Mittal et al. (2024).One possible explanation is that, with the additional dataset afforded by the joint analysis here, the statistical power of the result increases -noting that the analysis of Mittal et al. (2024) dealt with Quaia alone.Further, while the joint analysis points towards an excessive dipole amplitude, there are still additional questions regarding the full effect of local structure to pursue.We therefore do not claim that the state of the tension is unequivocal, clear-cut or well-understood.
Importantly, we have shown that the impact of clustering in both NVSS and RACS-low is significant.Contamination from local radio sources, even if they represent a small part of the sample, disproportionately biases the inferred dipole amplitude to larger values.This necessitates a more detailed profiling of local sources in future studies, possibly including calculations of the clustering dipole expectation as some authors have done previously.We are therefore of the view that the common assumption that clustering is negligible is no longer tenable, and may at least in part account for a fraction of the dipole excess reported in many radio galaxy studies.One correlated question is why our inferred 'clustering dipole' (the dipole in the cross-matched sources) appears to align with the CMB dipole.Is this consistent with the predictions of the CP, or is it an anomaly?
With respect to our methodology, there are important factors to consider about our cross-matching, in which we pair near-infrared (2MRS) and radio (NVSS/RACS-low) sources.Generally, crossmatching radio catalogues with other optical/infrared samples is a non-trivial and well-known problem.This is largely because radio sources are often extended with multiple spatially distinct components, and it can be unclear whether the radio emission is even associated with the optical/infrared counterpart source itself (Fan et al. 2020a; see e.g.van Velzen et al. 2012 andFan et al. 2020b for recent attempts at cross-matching manually and automatically respectively).Roughly 10% of radio sources will have extended, spatially-separated lobes, and it can be challenging to distinguish between these structures (part of the same source) and genuinely distinct sources.It is therefore conceivable that by cross-matching NVSS/RACS-low with 2MRS, we occasionally only remove one of the many lobes of a radio source and miss the extended structure.In a follow up analysis, we tried to probe this issue and remove highly-extended sources by a cut in the ratio of the integrated to peak flux density  int./ peak .We tried cuts from  int./ peak > 10 to as low as  int./ peak > 5, but found that this has a negligible impact on the inferred parameters and model marginal likelihoods.It would therefore be worth exploring this in more depth in future analyses.
One additional point relates to the actual value of the crossmatching radius itself, which we set at 5 arcseconds in light of the positional accuracy of NVSS/RACS-low objects.A more thorough way of setting the radius would involve constructing a 'fake' catalogue with random positional offsets from NVSS/RACS-low, then evaluating the fraction of genuine versus fake matches as a function of radius (see e.g.Mahony et al. 2011;Chhetri et al. 2020).This leaves scope for future work to investigate more robust and consistent ways of removing local galaxies and handling the realities of extended radio structure, beyond what we have already achieved with our choice of mask and cross-matching.
Finally, our joint methodology invites extension with other catalogues.The utility of this approach is that we may add additional terms to equation (10) for each other sample, and these need not be limited to radio galaxies.For example, we could include samples of quasars or Type Ia supernovae.Doing so improves the statistical power of the result by analysing independent, unique objects in which different processing techniques have been used.While an unaccounted for systematic effect could bias the result, it would have comparatively less of an effect amongst the other samples in the joint analysis.In future work, this will help further understand the nature of the claimed tension between the CMB dipole and the matter dipole in the literature, robustly interrogating the role of the CP in our cosmological paradigm.

Figure 1 .
Figure 1.Deviation from NVSS(B) mean source density (%) by declination angle for different flux density cuts.The declinations using the DnC configuration are shown in the shaded light brown region, and the D configuration region is left unshaded.The flux density cuts chosen are labelled on the right of each panel.

Figure 2 .
Figure 2. Flux density distribution of the NVSS(B) sample, shown in red, in which a 15 mJy flux limit has been used.Overlaid onto this in yellow is the integrated flux density distribution i.e. the total number of sources above a limiting flux density.

Figure 3 .
Figure 3. Deviation from RACS-low(B) mean source density (%) by declination angle for different flux density cuts.The flux density cuts chosen are labelled on the right of each panel.

Figure 4 .
Figure 4. Flux density distribution of the RACS-low(B) sample, shown in blue, in which a 15 mJy flux limit has been used.Overlaid onto this in yellow is the integrated density flux distribution i.e. the total number of sources above a limiting flux density.
(i) For each original NVSS or RACS-low source  with flux density   , we generate a new flux value  *  by drawing from a Gaussian G( =   ,  = Δ  ) for flux density measurement uncertainty Δ  .

Figure 5 .
Figure 5. Sky projection (Mollweide) of our final cross-matched (B) samples in Galactic coordinates, binned into healpixels ( side = 64; pixel angular size of 55').Masked regions are indicated in grey.Left: NVSS(B).Right: RACS-low(B).Top row: Number counts for each pixel.The mean source count is given by the middle tick of the colour bar.Bottom row: Averaged counts over a 1 steradian scale.The mean over all averaged counts is given by the middle tick of the colour bar.

Figure 6 .
Figure 6.Computed expected amplitude D CMB for NVSS and RACS-low with different choices of  0 .The desired value of D CMB is the extrapolated intercept at  0 = 15 mJy, as indicated by the dashed line.

Figure 7 .
Figure 7. Corner plots for the free dipole model tested on our individual samples.The dashed lines on the 1D marginal posteriors represent a 95% credible interval i.e. 2 either side of the median.Quoted uncertainties above each 1D posterior also give 2 credible limits on each model parameter.The contours on the 2D marginal posteriors give intervals of 0.5, encompassing 11.8%, 39.4%, 67.5% and 86.4% of the distribution.Top row (red): NVSS(A) and NVSS(B) samples.Bottom row (blue): RACS-low(A) and RACS-low(B) samples.

Figure 8 .
Figure 8. Corner plot for the NVSS(B) and RACS-low(B) analysis.Although the free dipole is not the strongest model, we show its posterior distribution here for illustration.The dashed lines represent a 95% credible interval, and the quoted uncertainties above each 1D posterior also give 2 credible intervals on each model parameter.The contours on the 2D marginal posteriors give intervals of 0.5, encompassing 11.8%, 39.4%, 67.5% and 86.4% of the distribution.D N refers to the NVSS dipole amplitude (D NVSS ), and D R refers to the RACS-low dipole amplitude (D RACS ).

Figure 10 .
Figure 10.Inferred dipole amplitude against maximum redshift removed during cross-matching.The error bars represent a 95% (2) credible interval.Red is NVSS, blue is RACS-low, and pink is the joint sample (the relevant samples are also labelled to the right of each pane).The brown line indicates the expected dipole amplitude in each case, and the grey line shows the inferred amplitude with no sources removed (i.e. the amplitude at  = 0).
Spatial distribution (Galactic coordinates; Mollweide projection) of NED and 2MRS sources ( < 0.02) and those that are cross-matched in NVSS/RACS-low.The black star indicates the direction of the CMB dipole.Some sources may reside in masked NVSS/RACS-low regions since crossmatching is applied before masking.Top: NED and 2MRS sources with redshift less than 0.02.The graticule labels have been removed for clarity.Middle: All 3049 cross-matched NVSS sources removed in NVSS(B).Bottom: All 3700 cross-matched RACS-low sources removed in RACS-low(B).

Figure 12 .
Figure 12.Projection of the marginal posterior distribution ( and ) for the cross-matched source dipole onto the sky (Mollweide).Blue indicates the RACS-low clustering dipole and red the NVSS clustering dipole.The contours give intervals of 0.5, encompassing 11.8% and 39.4%.Additional contours have been omitted for clarity.The black star indicates the direction of the CMB dipole.

Figure 13 .
Figure13.Marginal posterior distributions for the implied velocity of our heliocentric frame with respect to the CMB in the joint analysis (B variants).The distributions have been smoothed with a Gaussian kernel for the sake of visualisation.The top pane shows the inferred amplitudes assuming the kinematic direction model ( 3 ), while the middle and bottom panes assume the free dipole model ( 1 ) and the W23 model ( 5 ) respectively.

Table 1 .
Centre position in galactic coordinates ( • ,  • ) and radius  • of the distinct disc masks used in NVSS and RACS-low.Bright radio sources contained within or at the edge of the discs are indicated in the 'A-source' column.

Table 2 .
≤ 1000 −77 • <  < 29 • || ≤ 5 Summary of final RACS-low and NVSS samples, including flux cuts, masking choices and total number of sources.(A) samples are the catalogues containing local sources, and (B) samples have had local sources removed.

Table 3 .
Wagenveld et al. (2023) (2023)Description of the six models used in this work labelled  0 through  5 .

Table 4 .
Bayes factors by model for each NVSS sample.The highlighted row indicates which model has the highest Bayes factor -hence marginal likelihood -for each sample.

Table 5 .
Bayes factors by model for each RACS-low sample.The highlighted row indicates which model has the highest Bayes factor -hence marginal likelihood -for each sample.