The 408 MHz Haslam map is widely used as a lowfrequency anchor for the intensity and morphology of Galactic synchrotron emission. Multifrequency multiexperiment fits show evidence of spatial variation and curvature in the synchrotron frequency spectrum but there are also poorlyunderstood gain factors between experiments. We perform a Bayesian model comparison across a range of scenarios using fits that include recent spectroscopic observations at sim 1GHz by MeerKAT. A large uncorrected gain factor of about 60 in the Haslam data is strongly preferred partly undermining its use as a reference template.
['Astrophysics - Cosmology and Nongalactic Astrophysics', 'Astrophysics - Astrophysics of Galaxies']
Bayesian evidence for uncorrected gain factors in Galactic synchrotron template maps
{'Bayesian evidence for uncorrected gain factors in Galactic synchrotron template maps': 'Michael J. Wilensky , 1 , 2 , † , ∥ ★ Melis O. Irfan , 3 , 4 Philip Bull , 2 , 4  \n1 Department of Physics and Trottier Space Institute, McGill University, 3600 University Street, Montreal, QC H3A 2T8, Canada  \n2  \nJodrell Bank Centre for Astrophysics, University of Manchester, Manchester, M13 9PL, United Kingdom  \n- 3 Institute of Astronomy, Madingley Road, Cambridge CB3 0HA, UK\n- 4 Department of Physics and Astronomy, University of Western Cape, Cape Town 7535, South Africa\n- † CITA National Fellow  \n∥ TSI Postdoctoral Fellow  \n12 September 2024', 'ABSTRACT': 'The 408 MHz Haslam map is widely used as a low-frequency anchor for the intensity and morphology of Galactic synchrotron emission. Multi-frequency, multi-experiment fits show evidence of spatial variation and curvature in the synchrotron frequency spectrum, but there are also poorly-understood gain factors between experiments. We perform a Bayesian model comparison across a range of scenarios, using fits that include recent spectroscopic observations at ∼ 1 GHz by MeerKAT. A large uncorrected gain factor of about 60% in the Haslam data is strongly preferred, partly undermining its use as a reference template.  \nKey words: methods: data analysis - methods: statistical - (cosmology:) diffuse radiation', '1 INTRODUCTION': "Synchrotron emission from our own Galaxy is a pernicious contaminant of cosmological surveys at radio and microwave wavelengths, particularly those making total intensity (and polarisation) maps of the Cosmic Microwave Background and 21cm line emission. Fortunately, Galactic synchrotron emission is expected to derive from free electron populations with power-law energy distributions, leading to a spectral energy distribution (SED) along each line of sight that closely approximates a power-law in frequency, 𝑆 𝜈 ∝ ( 𝜈 / 𝜈 ref ) 𝛼 . Observed values of the power-law spectral index tend to reside in the range -1 . 2 ≲ 𝛼 ≲ -0 . 5, with some variation as a function of Galactic latitude (Planck Collaboration 2016c). This permits a simple model of synchrotron foreground emission to be constructed and subtracted from radio and microwave data, revealing the cosmological maps of interest up to some corrections for signal loss and residual contamination depending on the foreground removal method that is used.  \nAn important and long-standing ingredient of many foreground removal approaches is the Haslam 408 MHz all-sky survey (Haslam et al. 1982). The Haslam map anchors the overall amplitude of the synchrotron power-law model in each direction on the sky. Assuming a power-law SED and a suitable model of the spectral index variation across the sky, the synchrotron intensity can then be predicted across a wide range of frequencies. Multiple sky models, such as the Planck Sky Model (Delabrouille et al. 2013) and Python Sky Model (Thorne et al. 2017), rely on the Haslam data as a proxy for all-sky synchrotron emission amplitude in this way. Additionally, emission models such as the Global Sky Model (Zheng et al. 2017), which use empirical data across a wide range of frequencies, still rely on  \nthe Haslam data to provide a low frequency, high angular resolution (56 arcmin) synchrotron emission estimate. Similarly, Bayesian fits to empirical data, such as the Commander fits to the WMAP and Planck data (Planck Collaboration 2016a), use the Haslam data as a prior for synchrotron emission amplitude, with a normalization factor that changes this amplitude across their frequency range. SEDs constructed from multiple experiments across frequencies from tens of MHz to tens of GHz have detected a small degree of curvature in the power-law along some lines of sight (Platania et al. 1998; de Oliveira-Costa et al. 2008; Kogut 2012; Mozdzen et al. 2019; Irfan et al. 2022), and a frequency-dependent correction to the spectral index of the form 𝛼 → 𝛼 + 𝐶 log ( 𝜈 / 𝜈 ref ) has been suggested in Kogut (2012). This hints at energy losses through inverse Compton scattering for the relativistic charged particles responsible for synchrotron emission within our Galactic magnetic field (Strong et al. 2011), and also has implications for the accuracy of the synchrotron model at frequencies significantly higher and lower than 408 MHz.  \nThere are known flaws in the Haslam map however. For example, the original map contains unsubtracted (or imperfectly subtracted) point sources, and striping caused by imperfect filtering of correlated noise in the time-ordered data. Many of these issues were addressed in a reprocessing of the original data (Remazeilles et al. 2015), but some hints of residual systematic effects remain. For instance, Planck Collaboration (2016a) note that the synchrotron emission amplitude at 408 MHz requires correcting for use at higher frequencies if the spectral index is to be treated as constant over those frequencies, while Planck Collaboration (2016c) note that a single spectral index value across all frequency ranges can't be used to fit empirical data without generating 'large and spurious gain corrections at 408 MHz.' However, Remazeilles et al. (2015) point out that there is a large possible source of calibration uncertainty associated with the Haslam data in the form of an unknown relationship between the angular scale  \nof the map and the brightness temperature scale, due to power outside of the primary beam. This beam efficiency factor can be upwards of 30% for a radio experiment (Planck Collaboration 2014; Du et al. 2016). An important question is whether such issues could have led to incorrect inferences in some of the (many) studies that use the Haslam map as an anchor, e.g. regarding the presence of curvature in the synchrotron SED, or over-/under-subtraction of foreground contamination.  \nIn this paper, we use the tools of Bayesian model comparison and new spectroscopic radio observations at ∼ 1 GHz to assess how reasonable gain biases may restrict the ability of sub-GHz data alone to constrain the synchrotron spectral index curvature. We restrict ourselves to considering an overall gain factor for each survey as the only systematic effect in the experiments.  \nIn Wilensky et al. (2023), a Bayesian jackknife framework called Chiborg was introduced, wherein an analyst proposes a discrete selection of hypotheses to be compared. Formally, the hypotheses represent hierarchical models (Gelman et al. 2021) with common hyperparameters. In each hypothesis, hyperparameters are adjusted to represent the possibility that a subset of the data are somehow biased. The marginal likelihood, often called the evidence , is then calculated for each hypothesis, and these are normalized to create a posterior probability mass over them. For small data sets, or with carefully chosen hyperparameters, one can perform an exhaustive search over all of the possible bias configurations and determine which are most likely. A null hypothesis is always available where none of the data are biased.  \nWe slightly depart from the exact formalism in Wilensky et al. (2023), proposing that the data in question are susceptible to multiplicative gain biases, rather than additive offsets. Since the Chiborg framework compares hypotheses by their marginal likelihoods and the underlying model for the data prohibits analytic marginalization, we use the nested sampler polychord (Handley et al. 2015a,b) to simultaneously calculate the evidence for each hypothesis and constrain free parameters. The evidence calculation allows us to determine which experiments, if any, are likely have uncorrected gain calibration errors, while the parameter constraints allow us to estimate their size.  \nWhile it has long been known that the Haslam gain errors might deviate from the usual estimate of 10% uncertainty (e.g. Remazeilles et al. 2015; Monsalve et al. 2021), in this paper we show that there is strong statistical evidence that the gain errors may be as high as 60% on degree scales in some regions. This has potentially serious implications for analyses that rely on the Haslam map, such as joint CMB foreground component separation (Planck Collaboration 2014, 2016a,c; BeyondPlanck Collaboration et al. 2023), and foreground removal strategies in 21-cm experiments (Eastwood et al. 2019; Pagano et al. 2024).", '2.1 Synchrotron SED data and models': 'For this analysis we use the same frequency range and regions on the sky as probed in Irfan et al. (2022): the OVRO Long Wavelength Array (OVRO-LWA) 73MHz data, the Haslam 408 MHz data, and the MeerKLASS 971-1075 MHz data within 1 . 8 · regions centered at (RA, Dec) = (161 · , 2.4 · ), (164 · , 2.7 · ) and (167 · , 3.5 · ). In galactic coordinates, ( 𝑙 , 𝑏 ), these are (246.5 · , 50.7 · ), (249.4 · , 53.1 · ) and (252.1 · , 55.8 · ). We denote these regions as fields 0, 1 and 2, respectively. The data have all been smoothed to a common resolution of  \n1 . 8 · with Gaussian beams. Both the OVRO-LWA (Eastwood et al. 2018) and MeerKLASS (Santos et al. 2017) experiments are recent (within the last decade) 21cm intensity mapping endeavors; the former is interferometric, while the latter is used in single-dish mode. For the Haslam values, we use the destriped but not desourced map from Remazeilles et al. (2015), since neither the MeerKLASS nor OVRO-LWA measurements have been desourced (Irfan et al. 2022). We use an aperture photometry method to extract SED data points at each frequency from the available maps. This means each SED data point at each frequency is an average over the pixels in a given field in a way that leaves us insensitive to the absolute zero-level of the maps, but sensitive to multiplicative factors.  \nFollowing Irfan et al. (2022), we assume as a null hypothesis that the data follow a curved power law,  \n𝑆 ( 𝜈 ) = 𝑆 ( 𝜈 0 ) GLYPH<18> 𝜈 𝜈 0 GLYPH<19> 𝛼 0 + 𝑐 ln ( 𝜈 / 𝜈 0 ) , (1)  \nwhere 𝜈 is the frequency, 𝜈 0 is a reference frequency, 𝛼 0 is the spectral index at the reference frequency, and 𝑐 is the curvature. Each experiment quotes a fractional calibration uncertainty. OVRO-LWA quotes between a 5% and 10% gain calibration error based on the flux measurement of 11 known calibrators, including Cygnus A, Cassiopeia A, and Virgo A. The Haslam gain calibration error is often quoted as 10%, which relates to the level of scatter seen when correlating this data set with an absolutely calibrated map at 404 MHz (Pauliny-Toth &Shakeshaft 1962). MeerKLASS calibrate their data with respect to a system temperature model which consists of Galactic and terrestrial contributions; if this calibration were perfect the difference between the data and the model would simply be Gaussian distributed noise centred around zero. As they observe data-model residuals of 0.3K for a 16K system temperature model, they quote their calibration error as good to 2%. The gain uncertainties we choose to carry forward for OVRO-LWA, Haslam, and MeerKLASS are 5%, 10% and 2%, respectively. Each experiment also quotes a thermal noise uncertainty that is very small compared to the calibration uncertainty. Nonetheless, we model the data as  \nd = g · s + n , (2)  \nwhere the 𝑖 th component of each vector refers to the 𝑖 th measured frequency, d is the data, g are the gains, s is the underlying frequency spectrum, n is the noise, and · represents an elementwise product.  \nWe assume that the noise is Gaussian and uncorrelated,  \nn ∼ N ( 0 , N ) , (3)  \nwhere N is a diagonal matrix containing the quoted thermal variance of each experiment. Conditional on knowing the gains and underlying power law, the data are then specified as  \nd | s , g ∼ N ( g · s , N ) . (4)  \nIn the absence of prior information about the gain amplitudes other than their quoted uncertainties, we assume, under the null hypothesis H 0 , an uncorrelated Gaussian prior centered on 1, with diagonal covariance G that holds the squares of the quoted fractional errors,  \ng | G , H 0 ∼ N ( 1 , G ) . (5)  \nWe can marginalize over this prior to obtain  \nd | s , N , G , H 0 ∼ N GLYPH<16> s , N + ss 𝑇 · G GLYPH<17> , (6)  \nwhere s 𝑇 is the transpose of the vector s . We then propose 31 alternate hypotheses. Each hypothesis in the full set of 32 can be thought of in terms of a 5-bit string. One of these bits corresponds to the curvature,  \ni.e. in each hypothesis we either specify that the spectral curvature is exactly 0, or we allow it to range from -0 . 3 to 0 with a flat prior. The other four bits correspond to gain offsets. There is one bit each for OVRO-LWA and Haslam having a possible gain bias. For the MeerKLASS data, we split into a lower and upper subband. The division point was chosen by eye at 1007.34 MHz, roughly 2/5 the way through the full MeerKLASS band. This point was chosen since there is a visually apparent discontinuity in the SED at this point that is consistent in all three fields. We propose one overall gain bias for each subband. Denoting B 𝑗 as the set of hypotheses where the 𝑗 th suband has a non-unity gain factor, i.e. the 𝑗 th gain bit is \'active,\' the corresponding bias 𝜀 𝑗 is assumed to have a prior distribution  \n𝜀 𝑗 |H ∈ B 𝑗 ∼ N ( 0 , 0 . 25 ) . (7)  \nOtherwise, we assume that bias is exactly 0.  \nThis choice of prior width is chosen to balance two concerns. First, for a useful test, hypotheses must be distinguishable in the sense that one must be able to reliably (i.e. on average) tell the difference between data generated from one model versus the other (formally explored in Wilensky et al. (2023) in terms of mutual information ). As a pedagogical example, consider a coin flip experiment where two coins are each flipped 10 times and their results are compared. It is essentially impossible to reliably tell the difference between an experiment where the probability of heads for both coins is exactly 50% and an experiment where the probability of heads might vary independently between 45 and 55% between the two coins. There simply is not enough information to reliably resolve the difference between these two very similar models in such a sparse data set even when results are extreme (e.g. 1 heads and 9 tails for both coins). On the other hand comparing the model of two perfectly fair coin tossing processes to one where each heads probability might vary independently all the way from 0 to 1 can produce a reliable decision rule where extreme results strongly favor the \'potentially unfair coins\' hypothesis and balanced results (heads and tails in roughly even proportions) strongly favor the \'fair coins\' hypothesis, with a small region of ambiguity between these two possibilities. In analogy to this example, it is important to choose gain bias priors that are sufficiently broad to produce this distinguishability, with the understanding that favoring the null hypothesis (gain bias exactly 0) may also be desirable in situations involving small gain biases with no statistically meaningful difference of interpretation. However, the second concern is that choosing an excessively broad prior implies a significant belief in unreasonable (or even unphysical) situations, which is clearly undesirable.  \nThe gain hierarchy is constructed so that  \ng | 𝜺 ∼ N ( 1 + 𝜺 , G ) , (8)  \nwhere 𝜺 is a vector of all the gain biases (tiled appropriately to match the corresponding frequencies for each experiment) and 1 is a vector of all ones. This implies that  \nd | s , 𝜺 , N , G , H ∈ B ∼ N GLYPH<16> s · ( 1 + 𝜺 ) , N + ss 𝑇 · G GLYPH<17> , (9)  \nwhere B is the set of all hypotheses where any bias exists.  \nIn all hypotheses, we adopt a broad, flat prior on the 𝑆 ( 𝜈 0 ) and 𝛼 0 parameters, using 𝜈 0 = 73 MHz as the reference. The prior for 𝑆 ( 𝜈 0 ) extends from half the reported OVRO-LWA SED value to twice its value. The 𝛼 0 prior extends from -1 . 8 to 0.  \nWe examine these hypotheses using Bayesian model comparison, wherein we calculate the posterior probability of the hypotheses, 𝑃 (H| d ) , and analyze the data in light of what is most likely. This offers a formal way of comparing models that also naturally accounts for overfitting (Kass & Raftery 1995; Jaynes 2003; Mackay 2003),  \nthough it is not without complications (Gelman et al. 2021). By Bayes\' theorem,  \n𝑃 (H| d ) = 𝑃 ( d |H) 𝑃 (H) " H 𝑃 ( d |H) 𝑃 (H) (10)  \nThe term 𝑃 ( d |H) is the marginal likelihood, or evidence, of the hypothesis, H . Representing the model parameters in any given hypothesis as 𝜽 , the evidence is given by  \n𝑃 ( d |H) = ∫ d 𝜽 𝑃 ( d | 𝜽 , H) 𝑃 ( 𝜽 |H) , (11)  \nwhere 𝑃 ( 𝜽 |H) is the prior for the model parameters in the hypothesis H . The term 𝑃 (H) represents the prior probability of a given hypothesis. This allows one to disfavor different hypotheses in the selection process if they represent unlikely scenarios according to a particular state of knowledge.', '2.2 Numerical implementation': "The nonlinearity of the model in the power law parameters leads to an analytically difficult evidence integral (Equation 11). We instead use the nested sampler polychord (Handley et al. 2015a,b) to simultaneously compute the evidence of each hypothesis and constrain the remaining degrees of freedom. For the sake of reproducibility, we state our choice of polychord parameter settings. polychord 's evidence estimation is affected by the number of live points during sampling, 𝑁 live , which should be scaled according to the dimensionality of the problem, 𝑁 dim . Different hypotheses have different dimensionality. We use 𝑁 live = 50 𝑁 dim . The polychord parameter 𝑁 repeats controls how many slice sampling steps are performed to generate a new sample, which affects the correlation length of the chain. For reliable evidences, it is recommended to use 𝑁 repeats ∼ 5 𝑁 dim . Wefollowed this recommendation and found the inference to be very similar between several independent chains.", '3 HYPOTHESIS COMPARISON': "In order to form a posterior probability over the discrete hypotheses, we must assign a prior to each hypothesis. The simplest option is to assign equal prior probability to each hypothesis, shown in Figure 1, along with 1 𝜎 Monte Carlo error estimates. Since polychord estimates the logarithm of the Bayesian evidence, the errors are larger for stronger hypotheses. With a flat prior, the posterior is dominated by a few hypotheses in each field, concentrated on 8 hypotheses where Haslam suffers from a gain bias in every one, and curvature is disallowed. Field 1 demonstrates mild confusion about whether the OVRO-LWA datum has a gain bias. At least one biased MeerKLASS band is implicated in all dominant hypotheses. Hypotheses with nonzero curvature are strongly disfavored.  \nWhether (or not) models with and without curvature can be distinguished over this frequency range is an important concern. To answer this question, we simulated a small number of data sets with randomly determined 𝛼 0 , 𝑆 ( 𝜈 0 ) , and gain biases, all from distributions much more concentrated than the Chiborg priors. In half the simulated data sets, we set the curvature to 0, and in the other half we assigned a curvature randomly, concentrated on the value -0 . 09. We ran half the simulations with a range of anchor frequencies significantly higher than the data (3, 10, 30, and 100 GHz) while holding other model parameters fixed so that we could investigate the effect of having high-frequency anchors. We then used the same priors as we did on the actual data, and evaluated the marginal posterior probability of there being nonzero curvature, expressed as an odds ratio. To  \nFigure 1. Posterior probability over the discrete hypotheses for each field, including Monte Carlo estimation errors. Points within a gray band assume a particular hypothesis. Each corresponding label on the horizontal axis indicates which gain biases are active (Haslam, H; LWA, L; MeerKAT sub-bands 1 and 2, M1 and M2), and whether curvature is allowed (C) in the hypothesis.  \n<!-- image -->  \nbriefly summarize the results, we find that Chiborg 's ability to determine the presence of curvature in the simulations (unsurprisingly) increases rapidly if we add one more anchor at a significantly higher frequency. Using only MeerKAT as an anchor and only small gain biases, it is somewhat more difficult to distinguish between cases with and without curvature in simulation. However, we remark that the simulations are highly conservative, and that marginal odds against curvature are on the order of 30 in each field. This suggests that the actual data are unlikely to be represented by the highly conservative simulations, and that an anomaly is likely present in the data.  \nSince in the alternate hypotheses we propose rather large gain biases, one could assign a small prior probability that this occurs in any given experiment to express skepticism, and then assume that separate experiments commit this type of error randomly. For probabilities less than 0.5, this prior disfavors hypotheses at a rate that is exponential in the number of nonzero biases. We find that one must express decisive skepticism in order to make significant changes compared to Figure 1. For instance, assigning a (very low) prior probability of a bias existing of 1% keeps the mass concentrated on the same 8 hypotheses, but the emphasis shifts towards a more balanced weighting.", '4 FITS AND PARAMETER CONSTRAINTS': "A fully Bayesian estimate of the parameters is formed by marginalizing over the hypotheses. Denoting the full set of parameters as 𝜽 ,  \n𝑃 ( 𝜽 | d ) = ∑︁ H 𝑃 ( 𝜽 , H| d ) = ∑︁ H 𝑃 ( 𝜽 |H , d ) 𝑃 (H| d ) (12)  \nwhere 𝑃 (H| d ) is plotted in Figure 1. The full posterior over the parameters is therefore a mixture over the various hypotheses, and results in weakly multimodal parameter estimates in some instances. Field 0 is almost completely dominated by one hypothesis, which incidentally is the most extreme violation of the null hypothesis. We  \nshow an updated fit to the SED data in Field 0 in Figure 2 using just this hypothesis for the parameter estimates. In this hypothesis, the gain bias estimates are rather large for all of the experiments. The OVRO-LWA, Haslam, MeerKAT lower subband, and MeerKAT upper subband gain bias estimates for this hypothesis are shown in the first row of Table 1. Note the slight difference in bias estimate for the two MeerKAT subbands, which helps smooth out the apparent discontinuity in the uncorrected data.  \nIn Figure 2, we show the SED data calibrated by each experiment independently, along with 'corrected' data, which is the original data divided by the posterior mean of the gain bias for each point. Wethen show the 1 𝜎 posterior quantile and median for the power law model values under the dominant (dashed with grey shade) and null (dotted with grey shade) hypothesis. At each frequency with data, we add a colored 1 𝜎 error bar to represent the remaining uncertainty in the gain realization and thermal noise assuming these particular parameter values. In other words, we have regressed for the mean of g , given by 1 + 𝜺 , but not g itself. The quadrature sum of these two error quantiles roughly corresponds to the 1 𝜎 quantile of the full posterior predictive distribution of the corrected data.  \nThe posterior width of the power law model is much larger when we introduce uncertainty in the gain bias parameters. While the null hypothesis expresses relatively extreme certainty about the power law model, particularly at MeerKAT frequencies, over half of the data lie outside the 1 𝜎 quantile indicated by the gain uncertainties, which dominate the posterior predictive uncertainty for the null hypothesis. For the dominant hypothesis, slightly less than 1/3 of the data lie outside the 1 𝜎 quantile of the best fit predictive distribution, which is roughly within expectation in the instance that the chosen parameter values are correct. The model is clearly not perfect, as there appear to be correlated errors in the MeerKLASS data, i.e. the corrected data seem to consistently over-/under-shoot the model over frequency scales of roughly 20 MHz. A more complete analysis might propose a non-diagonal G parameter and regress for the correlation length.  \nWe show a corner plot of the posterior samples for Field 0 in Figure 3. The power law parameters exhibit a strong correlation with  \nFigure 2. Data and fits in Field 0 using the null hypothesis and dominant hypothesis reported by Chiborg . The 'corrected data' are the data divided by the best fit gain. The 'best fit predictive' is the 1 𝜎 interval for the predictive distribution of the data conditional on the best fit gain bias and model parameters. The grey shades show the 1 𝜎 quantile for the posterior distribution of the predicted power law. The inset shows the MeerKAT data in more detail.  \n<!-- image -->  \nTable 1. Summary of posterior parameter constraints for MAP hypothesis in each field.  \none another, as well as with the OVRO-LWA gain bias. The reference flux is also slightly correlated with Haslam's gain bias. Additionally, the two MeerKAT subbands have very tightly correlated gain biases. Other parameter pairings exhibit mild or no degeneracy. Despite these apparent degeneracies, which were to be expected in such a model, the bulk of the posterior mass still lies rather far away from possibilities where the experiments are unbiased. This reinforces the notion that the data have gain biases. For brevity, we only show PPDs and corner plots for Field 0. Other fields were similar, though with much less extreme gain biases for MeerKAT, and still relatively strong gain adjustments for OVRO-LWA and Haslam.  \nTable 1 shows parameter estimates using, for the sake of simplicity, only the maximum a posteriori (MAP) hypothesis for each field. We find somewhat steeper spectral indices compared to Irfan et al. (2022), however this is to be expected since spectral curvature is fixed at zero. Both OVRO-LWA and Haslam are estimated to have significant positive gain biases that are outside the original quoted gain errors for the experiments, however the OVRO-LWA gain bias has only ∼ 2 𝜎 significance in Fields 0 and 1.  \nOf particular noteworthiness is that the Commander synchrotron maps (Planck Collaboration 2016b), which have their own gain factor accounting, also adjust the 408 MHz SED downward in these fields,  \nFigure 3. Corner plot for the MAP hypothesis in Field 0. Dashed lines show the MAP parameter value (computed over the full joint posterior).  \n<!-- image -->  \nas visible in Figure 11 of Irfan et al. (2022). For MeerKAT, Field 0 is estimated to have a somewhat strong gain bias, but rather weaker ones for Fields 1 and 2. While roughly within the original quoted gain error for Fields 1 and 2, at face value this can be interpreted as a statistically significant detection of a small but systematic gain bias that affects the entirety of one of the two subbands in either case.  \nApotential issue with this interpretation is that there is significant variation between the different fields. Depending on how the data were calibrated, this variation may be inexplicable as an experimental gain. The MeerKAT data were calibrated as described in Wang et al. (2021), and while it is not implausible to have field-dependent calibration errors, alternatives exist for spatial variation such as unmodeled sources and RFI. Given these alternatives, it is hard to decide the source of these errors without additional data analysis. Similarly, direction-dependent calibration in Eastwood et al. (2018) and other aspects of the calibration routine (such as errors in the beam model) could plausibly give rise to field-to-field variation in the OVRO-LWA gain bias, though these effects are folded into the original 5% quoted error. We remark that detection of the gain biases in the OVRO-LWA data is less significant in the sense that the constraints are much broader and models without an OVRO-LWA gain bias are less disfavored than (e.g.) models where MeerKAT has no gain biases. In contrast to the other experiments, the inferred Haslam gain bias is stable across the three fields, removing the need to examine small-scale direction-dependent sources of error assuming they do not somehow align by coincidence.", '4.1 Robustness test': 'To test the robustness of these inferences, we performed the same exercise twice again, but with either the lower or higher MeerKAT subband entirely removed from the analysis. With the substantial reduction in data volume, we expect that there should be more confusion among the hypotheses along with broader parameter constraints, but that strongly disfavored hypotheses should remain strongly disfavored. In other words, the broad conclusions should be similar, but less sharply supported. The posterior probability over the hypotheses for each test run are shown in Figures 4 and 5 for including only the lower and upper subbands, respectively. We again find that the posterior is strongly concentrated on hypotheses where Haslam has a gain bias and curvature is disallowed. There are apparent differences between inferences regarding the OVRO-LWA and MeerKAT data, which we discuss in more detail.  \nFocusing on the low subband test, Field 2 has no significant preference for the configuration of gain biases among the OVRO-LWA and MeerKAT data. Comparing to Figure 1, we see that the MAP hypothesis is one where the low subband has no bias, though there is roughly equal probability in hypotheses where the low subband is biased if we marginalize over all possibilities. This ambiguity regarding the low subband is therefore present even when we consider the high subband data simultaneously. The ambiguity among OVROLWAhypotheses is not present in Figure 1, i.e. it is a new feature that emerges when the high subband is excluded. Meanwhile, Field 0 is relatively unsupportive of a gain bias in the lower MeerKAT subband when the higher one is excluded, and vice-versa in Field 1. Again referring to Figure 1, we see that a gain bias is only supported in', 'Low subband only': 'Figure 5. Same as Figure 1, but only using the higher of the two MeerKAT subbands.  \n<!-- image -->  \nH  \nHypothesis (Haslam biased, nonzero curvature, LWA biased, MK1 biased, MK2 biased)  \nFigure 4. Same as Figure 1, but only using the lower of the two MeerKAT subbands.', 'High subband only': '<!-- image -->  \nH  \nHypothesis (Haslam biased, nonzero curvature, LWA biased, MK1 biased, MK2 biased)  \nthe lower subband in Field 0 when one is also included in the higher subband, though there is a small chance that only the higher subband has a gain bias. It is therefore perhaps understandable that a gain bias in the lower subband is somewhat unsupported when the higher subband is entirely excluded from the analysis. The balance of low subband gain bias hypotheses in Field 1 is symmetric between Fig-  \nure 1 and Figure 4. Similar to Field 2, both Field 0 and 1 are relatively ambiguous about whether the OVRO-LWA datum has a gain bias.  \nThe inference using only the higher MeerKAT subband is nearly the opposite. There is a dominant hypothesis in Field 2 wherein all available data have gain biases. This is essentially also reflected in Figure 1, where there is only one competitive hypothesis excluding  \nany of the data used in the higher subband test, and it requires the use of the lower subband data (which was unavailable in this test by construction). Field 0 is in strong support of a MeerKAT gain bias, consistent with Figure 1 where all relevant hypotheses include a gain bias in the upper MeerKAT subband. Field 1 tends more towards excluding a gain bias in the upper MeerKAT subband. Once again consulting Figure 1, the MAP hypothesis involves an absence of a gain bias in the upper MeerKAT subband, though the marginal probability of a gain bias in this subband is roughly evenly split, demonstrating a discrepancy between this inference and the complete-data inference. The probability mass is again roughly equally distributed for and against an OVRO-LWA gain bias in Fields 0 and 1.  \nSummarizing, examining the MeerKAT subbands separately roughly preserves marginal inferences about the presence of MeerKAT gain biases, but introduces an ambiguity about the presence of an OVRO-LWA gain bias in most cases. The combination of the subbands therefore seems important for making inferences about the presence of an OVRO-LWA gain bias. We suspect this is because the broader spectrum of measurements gives the combined MeerKATdataastronger lever arm over the spectral index parameter, which is degenerate with the reference flux amplitude and OVROLWA gain bias parameters. Importantly, the separate examination preserves the conclusion that Haslam presents a gain bias and there is no curvature. In further support, the Haslam gain bias posteriors are centered at similar values as when all MeerKAT data are used, but the constraints are broader due to the presence of less data (not shown).', '4.2 Discussion': "Where might such a large calibration error arise in the Haslam map? With limited data, we can only speculate. An analysis conducted with several degree-scale absolutely calibrated maps of large sections of the sky, rather than just a few square degrees, could prove extremely useful in this regard.  \nIn any case, we mention one plausible source of calibration error that could give rise to an effect of this size. The conversion from brightness temperature to antenna temperature is scale-dependent, as discussed in Jonas et al. (1998). In particular, the solid angle of the source of interest must be taken into account when converting to brightness temperature. To account for this when studying diffuse emission, the data are often placed on the 'full-beam' scale, meaning the antenna temperature must be scaled by the ratio of the 'main beam' solid angle to the 'full beam' solid angle. The full beam scale usually does not refer literally to an integral out to the horizon, but rather out to some cutoff such as a few main beam widths (Haslam et al. 1974; Reich 1982; Jonas et al. 1998). Since the full beam solid angle is greater than the main beam solid angle, the full beam temperature is generally smaller than the main beam temperature. If the map were not put on the full beam scale, or if the full beam scale were significantly underestimated, then one might observe the size (and sign) of gain bias we claim here.  \nHaslam et al. (1974) claims only a 5% difference in the two scales (i.e. 95% of the beam power is in the main beam). This claim is based on 10 · sweeps surrounding Cas A (5 · in either direction). Compared to other instruments, where a 30% difference is more typical, this is a very low conversion factor that usually requires careful design and underillumination to acquire (see §4.4 of Remazeilles et al. (2015) and citations therein, as well as Dickinson et al. (2019)). It therefore seems unlikely that this 5% figure is accurate, which makes this a plausible explanation for why we might see a gain-like effect in the Haslam map as large as we do.", '4.3 Comparison to other results': 'These results uniformly favor the hypothesis that the Haslam map is too bright to be consistent with the other data, and that there is unlikely to be any spectral curvature in these fields once this is corrected for. This result is seemingly inconsistent with Monsalve et al. (2021), where globally scaling the Haslam map up by 21% and subtracting off a small 4.1 Kelvin additive offset improved agreement with absolutely calibrated EDGES data. Furthermore, allowing for spectral curvature had a similar effect on Haslam in that analysis (i.e. it produced agreement even without gain correction), whereas in our analysis we observe no such degeneracy. We argue that these analyses cannot be directly compared in a simple way due to position-dependent complications in the Haslam map and a general incongruity between the analysis in Monsalve et al. (2021) and ours, although this discrepancy does need to be resolved to understand the nature of any correction that would need to be applied to the Haslam map.  \nThe Haslam map is constructed from four different surveys made with different telescopes throughout the 60s and 70s (Haslam et al. 1982; Remazeilles et al. 2015). This makes it seem unlikely that a single rescaling over the entire sky should fix calibration errors that may be arising from different physical sources specific to particular surveys. For example, in the simple case where there is a single rescaling procedure required for each independent survey, regions of overlap will require different rescalings to regions of non-overlap. Since we are analyzing a tiny region of the sky compared to the more global analysis of Monsalve et al. (2021), it is unsurprising that these results may end up different (despite the fact that we observe consistent gain corrections for Haslam in each of our small, closely spaced fields). Moreover, given the nature of the Haslam map, it may be desirable to have a way of surgically searching for systematic calibration errors in a position-dependent way, such as with Chiborg . This will be particularly relevant for foreground modeling methods such as in Pagano et al. (2024) that partition the sky.  \nIn some fields on the sky, the synchrotron component can contribute a meaningful amount of brightness at the ∼ 30 GHz spinning dust peak (Dickinson et al. 2019; Rennie et al. 2022; Harper et al. 2024). As an example of how an analysis like this one could affect inferences about spinning dust, we take our best fit models under the MAP hypothesis (i.e. gain corrected models) for each field and extrapolate them to 30 GHz. We then compare this to an extrapolation with no curvature or gain corrections. The gain corrections tend to produce shallower spectral indices, resulting in a flux density ratio between the MAP model and no-correction model of 2.6, 1.6, and 2.0 for Fields 0, 1, and 2, respectively. In other words, supposing these data were used as low frequency anchors for inference about spinning dust emission at 30 GHz, there would be significantly different results without accounting for these statistically identified gain biases.', '5 CONCLUSIONS': 'Using a joint Bayesian analysis of data from three radio experiments (OVRO-LWA, the Haslam survey, and MeerKLASS), we found significant evidence for uncorrected gain factors in each experiment, with Haslam in particular showing a consistently large bias. The amount of MeerKAT data used notably affects whether we infer a gain bias in the OVRO-LWA data, which we suspect is related to the increased leverage over the spectral index that is obtained when combining all MeerKAT data. The presence and strength of these  \ncalibration errors appear somewhat field-dependent, except in the case of the Haslam map. An in-depth analysis of how the data were reduced for each experiment could show whether this is a sensible result in each case, but initial examination of the calibration methods of each experiment does not exclude this possibility.  \nWe emphasize that our model assumed that the only possible systematic effect was a gain factor, however one could imagine different types of systematic effect masquerading as large gain errors in this model, e.g. an additive bias from strong radio frequency interference. In any case, calibration errors or other departures from the astrophysical model such as these can cause significant discrepancies in the inferred spectral properties of synchrotron emission. Errors in the spectral properties have direct relevance to inferences regarding the underlying astrophysical processes (e.g. the presence of curvature or lack thereof).  \nThese discrepancies are also deleterious for cosmology experiments such as 21-cm intensity mapping surveys, which rely on foreground modeling for accurate calibration solutions. Even if an accurate calibration can be obtained through other means, such as in situ calibration, an incorrect foreground model could cause structure in the foregrounds to be absorbed by parameters related to the cosmological signal unless appropriate precautions are taken, such as jointly modeling the foregrounds with sufficiently broad a priori uncertainty (Anstey et al. 2021; Pagano et al. 2024). Worse yet, we find several experiments spanning at least one order of magnitude in frequency need to be combined in order to root out these gain biases. The possibility of such biases also makes spectral curvature and higher order spectral structure, formally speaking, difficult to assess without a very broad range of frequencies. Applying this Chiborg -based analysis to a selection of compatible maps (e.g. with appropriate angular resolution) formed from spectroscopic measurements spanning a large range of frequencies would allow for the formulation of a radio sky model complete with formal calibration uncertainties derived jointly among the contributing experiments. Such a sky model would help safeguard against inferential pitfalls in cosmology and astrophysics experiments that depend on accurate sky models.', 'ACKNOWLEDGEMENTS': "We acknowledge Clive Dickinson for very helpful discussions. This result is part of a project that has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 948764; MW, PB). MW was funded by a CITA National Fellowship. MI acknowledges support from the South African Radio Astronomy Observatory and National Research Foundation (Grant No. 84156). PB acknowledges support from STFC Grant ST/X002624/1.", 'DATA AVAILABILITY': 'The OVRO-LWA and (destriped) Haslam diffuse radio maps are available in the NASA-hosted LAMBDA archive. 1 Access to the raw MeerKAT data used in the analysis is public (for access information please contact archive@ska.ac.za). Data products are available on request. The python code we wrote for our analysis is available at https://github.com/mwilensky768/sed\\_jackknife .', 'REFERENCES': 'This paper has been typeset from a T E X/L A T E X file prepared by the author.'}