New Synoptic Observations of the Cosmic Optical Background with New Horizons

A truly cosmological optical background should present no more structure than would be expected given the known projected correlation function of large-scale structure marked by galaxies over the age of the Universe. More plainly, we expect the field-to-field variation of observed COB intensity to be more or less isotropic, once known instrumental effects and Galactic foregrounds have been corrected for. An ideal COB survey would attempt to uniformly sample the full sky on all angular scales. However, foregrounds associated with the Milky Way, or even the solar system, are markedly anisotropic over large angular scales. Thus any realistic survey will have additional Galactic and ecliptic coordinate constraints to minimize those foregrounds. Within such constraints, however, we are free to identify survey areas that specifically minimize both the DGL and SSL foregrounds, which are the two dominant foreground signatures for our observing platform.

In the case of the present survey, several considerations limit us to observing in somewhat restricted regions close to the Galactic poles. The strongest restriction is that needed to avoid scattered sunlight entering the LORRI aperture. This issue is discussed in detail in NH21; briefly, we require the aperture to be fully in the shadow of the NH spacecraft. This requirement means selecting fields with solar elongation angle (SEA) > 95°. While SEA > 90° would be sufficient to keep direct sunlight out of the LORRI aperture, the spacecraft bulkhead in which the aperture is positioned also supports other instruments that could potentially scatter sunlight into LORRI; SEA > 95° ensures that these are also shaded by the spacecraft (Weaver et al. 2020; NH21).

The trajectory of NH out of the solar system was completely specified by its primary mission of obtaining the first exploration of Pluto (Stern et al. 2015) and the Kuiper belt object Arrokoth (Stern et al. 2019). At the time of the mission, Pluto as seen from Earth was projected against the bulge of the Milky Way. This means that the "antisolar" hemisphere accessible to NH is roughly centered on the heart of our galaxy. Requiring Galactic latitude ∣b∣ > 40° to avoid dense stellar foregrounds and strong dust absorption thus eliminates a significant fraction of this hemisphere. Lastly, we restrict ecliptic latitude to ∣β∣ > 15°. While NH is not directly affected by ZL, the far-infrared (FIR) intensities that we use to select COB fields are provided by maps made in Earth-space, and thus may incur larger errors near the ecliptic (Matsuura et al. 2011; Carleton et al. 2022; Korngut et al. 2022).

The combination of these three constraints (SEA > 95°, ∣b∣ > 40°, and ∣β∣ > 15°) leaves 4239 deg² of sky available. We randomly selected 60,000 positions within this area, and for each one estimated the DGL contribution using the Improved Reprocessing of the IRAS Survey (IRIS) 100 μm all-sky data and the amount of SSL entering the LORRI field of view. The details of how we derived these preliminary foreground light estimates are discussed in NH22 (see Section 2.1 of that work). We then ranked each position by the sum of these two foreground contributions and selected the 500 positions with the smallest sums for further review. From these 500 positions, we reduced the sample to 15 fields, to ensure we could perform the desired observations within the available spacecraft fuel constraints. We also ensured that there was coverage in both Galactic hemispheres and that the pointings spanned a broad range in Galactic longitude. The final coordinates of each field were adjusted by up to 07 to minimize the presence of bright galaxies or galaxy clusters within the LORRI field of view. ¹⁴ We also used the NH22 test field (NHTF01) as part of the survey since, by definition, it meets all the above requirements. This brings the total number of COB survey fields to 16. With the exception of NHTF01, we denote the COB science fields with the prefix NCOB.

We selected an additional set of eight fields to perform an improved self-calibration of the relation of FIR intensity to optical DGL. The calibration fields for the FIR–DGL relation are denoted with the prefix DCAL. The DCAL fields were explicitly selected to cover fields with progressively higher 100 μm surface brightness, up to a limit of ∼3 MJy sr⁻¹. This limit was selected to avoid dust optical depths large enough that nonlinear behavior between the FIR intensity and scattered light amplitude might start to come into play. All DCAL fields were selected to have SSL levels in a narrow range (6 < SSL < 9 nW m⁻² sr⁻¹) that was similar to the NCOB science fields.

We also selected eight fields to help verify the SSL estimates. These fields are denoted with the prefix SCAL. The SCAL fields were chosen to be closer to brighter stars than we would otherwise permit for the science observations to test the reliability of the scattered light estimates. All SCAL fields were selected to have DGL levels in a narrow range (3 < DGL < 6 nW m⁻² sr⁻¹) that is similar to that in the NCOB science fields.

Lastly, we selected four fields at low ecliptic latitude (∣β∣ ≤ 6°) solely to verify the lack of significant brightness from interplanetary dust (IPD) at the large heliocentric distance of the NH spacecraft at the time of our observations. These four fields are denoted with the prefix IPDF.

The coordinates, observation dates, and spacecraft heliocentric distance for the 16 COB fields, 16 calibration fields, and four IPD fields are listed in Table 2. The mission elapsed time (MET) image identifiers of the first image of each field are also given. We show the field distribution on the sky with respect to the IRIS 100 μm map in Figure 1. The $30^{\prime} \times 30^{\prime}$ regions of sky from the DESI/DECam Legacy Survey (Dey et al. 2019, hereafter DLS) centered on our 16 COB science fields are shown in Figure 2. The intentional lack of bright optical sources in any of the fields is clearly demonstrated. Table 3 provides the Galactic extinction E(B − V), H i column density, Hα emission, dust temperature, and the FIR and cosmic IR background (CIB) intensities for each field. The CIB-subtracted FIR intensities are used to estimate the DGL foregrounds for each field. While we list the E(B − V) (Schlafly & Finkbeiner 2011) values for each field in Table 3 these extinction indicators are not used in any of the analyses in this work. Shull & Panopoulou (2024) have demonstrated that these E(B − V) values are, on average, 12% lower than those derived from Planck FIR observations. All errors listed are 1σ values.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Images of the survey fields were obtained with LORRI. With the exception of those for NHTF01, the observations were obtained over August and September of 2023 when the spacecraft was nearly 57 au from the Sun. For each COB survey field (those designated with the NCOB prefix), 16 LORRI exposures of 65 s each were obtained, while eight exposures of 65 s were obtained of the DGL calibration (DCAL) and scattered starlight calibration (SCAL) fields. The imaging sequences for each of the four IPD fields consisted of 16 LORRI exposures of 65 s each. To avoid the LORRI "background fade" anomaly associated with the activation of the camera (NH21), no images were obtained earlier than four minutes after LORRI was powered on.

LORRI is an unfiltered (white light) 1032 × 1024 pixel CCD imager mounted on a 20.9 cm aperture Cassegrain reflector. The active imaging area is 1024 × 1024 pixels, with the final eight columns being covered by a dark shield. The last four shielded columns are used to measure the combined dark current and electronic bias level. For deep observations, such as those used for COB measurements, the camera is operated with 4 × 4 pixel binning, producing raw images in 257 × 256 pixel format, including a single bias/dark column. The pixel scale in this mode is 408, which provides a $17\buildrel{\,\prime}\over{.} 4$ field. LORRI's sensitivity extends from the blue (0.4 μm) to near-infrared (0.9 μm) and is defined by the CCD response and telescope optics. The pivot wavelength is 0.608 μm. The gain is 19.4e⁻ per 1 DN, and the read noise is 24e⁻. In 4 × 4 mode, the photometric zero-point is 18.88 ± 0.01 AB magnitudes, corresponding to an exposure level of 1 DN s^–1 (Weaver et al. 2020).

The suitability of LORRI for COB observations is discussed at length in NH21 and NH22. In short, we have explored the spacecraft environment for potential extraneous contributions to the background sky level, as well as similar effects in the LORRI camera. As shown in NH21, the spacecraft shadow is adequate for preventing both direct and indirect solar illumination of the LORRI aperture. Analysis in NH21 also shows that the thrusters that control the spacecraft attitude do not generate particulates around the spacecraft that scatter sunlight. Lastly, γ-rays emitted by the spacecraft's radioisotope thermoelectric generator (RTG) or cosmic rays do not generate significant intensities of optical photons through the Cherenkov or fluorescent mechanisms (see NH22).

As LORRI is operated, the electronic bias level and average dark-current level are measured from the overscan column as a combined electronic background level. In NH21 we presented the analysis of a novel calibration sequence of exposures that allowed for direct isolation of the dark level. This demonstrated that the dark current was as expected, based on prelaunch calibration tests. In NH21 we further demonstrated that bias images obtained with LORRI produced a null background. In NH22 we discussed our discovery of a low-level error in the LORRI analog/digital conversion electronics, which caused the bias level to be in error by 0.02 DN. Subsequent analysis shows that it affects other signal measurements at the same level, thus it has no significant net effect and is ignored in the present analysis.

The sky levels in the images are only slightly greater than 1 DN. Accurate recovery of the sky intensity thus requires attention to subtle effects that become important at this level. As detailed in NH21 and NH22, we use a custom image reduction pipeline to optimize accurate recovery of the faint sky signal. This includes estimating the bias level by fitting a Gaussian to the peak of the DN histogram of the bias column. We also include a special step to correct for the "jail bar" pattern, in which the bias level of the even-numbered columns in the CCD is offset by either +0.5 or −0.5 DN from that of the odd-numbered columns (the sign of the correction varies randomly between LORRI power-on cycles). Lastly, we exclude bright cosmic-ray hits and negative amplifier "undershoot" artifacts associated with overexposed stars from the LORRI charge-smear corrections, as they are not smeared.

The procedures for measuring the sky levels are discussed extensively in NH21. In brief, we measure the sky for each individual exposure by first masking out foreground stars, galaxies, hot pixels, and cosmic-ray events, and then fitting a Gaussian to the peak of the intensity histogram of the remaining unmasked pixels. The histogram-fitting algorithm is designed to take into account fine-scale structure in the distribution of pixel intensity values that results from the image calibration operations applied to the initially integer raw pixel values. We emphasize that the final sky value for a given field is the the average of the individual sky levels measured for the 16 or eight images obtained of the field, as opposed to the single sky value measured from a stack of all the images. As noted in Weaver et al. (2020) and NH21, LORRI exhibits a slowly varying pattern of row-wise low-amplitude (<1 DN) streaks in its bias level. This pattern is treated as a random noise source, which is captured in the dispersion of sky values measured for any field.

In NH21 we discovered that images taken shortly after LORRI was powered on had elevated background levels, which appeared to decay away during the initial four minutes of operation. The cause of this background effect is unknown. The NH22 test COB images were therefore obtained only after this "cool-down" interval had elapsed following activation of the instrument. As noted in NH22, the sky levels in the 16 images obtained of the test field appeared to be constant over the sequence, validating this solution. We thus used the same procedure to obtain the present data.

Examining the sky levels measured for each image in our present richer data set, however, we have found that the background decay still continues even after the four minute delay, albeit at a low level. Comparison of the amplitude of the decaying background between the NCOB and DCAL exposures, which covers roughly a factor of two in total sky level, shows that the background excess was not tied to the exposure level. It is thus modeled as an additive effect.

All NCOB exposures comprise the same sequence of eight images of 65 s taken in rapid succession, followed by an 80 s pause to adjust the spacecraft pointing, followed by the final eight images, again taken in rapid succession. Subtracting the mean total sky level from the complete set of 16 images for any field showed that the first images had generally positive residuals compared to the average level over the sequence, with the final images having slightly negative residuals. Figure 3 shows the average residual trend for all the NCOB fields as a function of time.

Zoom In Zoom Out Reset image size

An exponential decay model appears to be an excellent description of the behavior of the residuals with time. We fitted the trace as

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{\rm{sky}}}(t)={{ae}}^{-t/\tau }+b,\end{eqnarray} \tag{ 1 }$$

where Δ_sky(t) is the average sky residual at any position on the NCOB exposure sequence, t is the time since the start of the sequence, and b is a constant background, which accounts for the fact that the initial mean sky for any field will include the background excess. For the NCOB sequence, a least-squares fit recovers a = 0.315 DN, τ = 295 s, and b = −0.075 DN, which corresponds to 1.72 nW m⁻² sr⁻¹ in intensity units. The background term in essence is the correction needed for the initial total sky levels to account for the presence of the decaying background. In practice, we use the model to apply a correction to each image in the sequence, which removes a source of variance in determination of the final average level, thus reducing the random error in the total sky measurements. We have also now applied this correction to NHTF01, with the caveat that its images were taken over a longer interval, thus reducing its net correction to −0.061 DN, or 1.39 nW m⁻² sr⁻¹.

The final decay model fitted to the uncorrected residuals is shown in Figure 3. We also include the sky residuals from the DCAL images in the figure. These fields have only eight images, and the shorter duration of the exposure sequence following the start of the sequence means that they are affected more strongly by the decaying background. The exponential model also provides the correction for the DCAL sky residuals, which are in excellent agreement with the residuals from the NCOB images taken at the same time lag. Use of the decay correction does add a systematic uncertainty of 0.16 nW m⁻² sr⁻¹ to the the sky level for any field, which is included in the total error budget.

We compute the total IGL in two steps: the bright IGL (hereafter BIGL) for galaxies with V < 19.9 that were masked during the sky estimation process and the faint IGL for galaxies below this LORRI detection threshold. The IGL for the bright galaxies (V < 19.9) is estimated by extracting nonstellar objects in our LORRI field of view from the DLS DR10 (Dey et al. 2019). We use NOIRLab's Astro Data Lab interface for retrieving the galaxy catalogs. The transformation to V-mag from the DLS g and r bands (or the g and i bands when the r band is unavailable) is derived from eight templates of galaxy spectral energy distributions spanning the morphology types E, S0, Sa, Sb, Sc, and Ir. We weight the templates by the morphological fractions observed in the field population of galaxies and derive an average (V − g) versus (g − r) relationship (or versus (g − i) when r-band information is unavailable) over the redshift range 0 < z ≤ 1, typical for brighter galaxies. Our best-fit color transformations are

$$\begin{eqnarray}\begin{array}{rcl}V-g & = & 0.6216-2.03(g-r)+0.7221{\left(g-r\right)}^{2}\\ V-g & = & 2.265-4.557(g-i)+2.238{\left(g-i\right)}^{2}\\ & & -0.3624{\left(g-i\right)}^{3}.\end{array}\end{eqnarray} \tag{ 2 }$$

These relations are valid over the ranges 0.6 ≤ (g − r) ≤ 1.5 and 0.4 ≤ (g − i) ≤ 2.8. We then derive the IGL intensity contribution based on the V magnitude and sum up the contributions for all DLS galaxies with V < 19.9 in the LORRI field of view. The statistical error for the bright IGL is ∼0.01 nW m⁻² sr⁻¹ and is derived from the photometric errors given in the DLS catalogs. The systematic error for the bright IGL, typically ∼0.07 nW m⁻² sr⁻¹, is derived from the uncertainties associated with the color transformation to the V band. As noted in Table 1, the BIGL addition for galaxies brighter than the LORRI detection limit was not applied in our NH21 analysis but was applied in NH22 and is applied in this work. The typical BIGL component for the fields used in NH21 is ∼1.2 nW m⁻² sr⁻¹.

The precepts for estimating the faint IGL due to galaxies at or below the V = 19.9 detection threshold are discussed at length in NH21. We estimate the uncertainty in the faint IGL intensity by assessing the specific contribution to the error from the systematic terms (errors in the fits to the galaxy number counts) and from the statistical errors (cosmic variance). The two systematic errors associated with the fits to the galaxy number counts are from the errors in the coefficients to the power-law fits used in NH21 and the error associated with the form of the fitting function. The formal errors in the power-law coefficients yield a fractional error of 13.1% in the IGL intensity. The difference between the IGL derived from the power-law fits and that derived using a quadratic fit to the galaxy counts yields a fractional change in the IGL of 6.6%. Summing these two error components in quadrature yields a combined systematic fractional error of 14.7% in the IGL intensity. The total error in the IGL must also include the statistical uncertainty due to the effects of cosmic variance over a single LORRI field of view (FOV). The cosmic variance error for a single LORRI FOV used in this work is the same as the single-field cosmic variance error adopted in NH21 (Trenti & Stiavelli 2008), which translates to an IGL fractional error of 11.8%. Summing, in quadrature, this statistical error with the above systematic error yields a total fractional error of 18.8% in the faint IGL intensity. Combining the bright and faint galaxy contributions to the IGL gives a total IGL intensity for each of our survey fields. This IGL corresponds to the expected light in the LORRI bandpass from all galaxies brighter than V = 30 mag. The bright and faint IGL values and their associated total errors are provided in Table 4 for each field.

The LORRI instrument accepts scattered starlight from sources well outside its field of view. The NCOB fields were selected to minimize SSL, and the SSL intensity remaining for any given field depends on the specific stars that surround it. As detailed in NH21, the LORRI scattered light function is estimated from prelaunch calibration tests and from in-flight measurements of the scattered sunlight background as a function of angular distance from the Sun. The SSL is estimated from the convolution of the scattering function with stars surrounding the field as provided by the Tycho2 star catalog (Høg et al. 2000), and the Yale Bright Star catalog v5.0 (Hoffleit & Warren 1995) for bright stars, with fainter stars (11 mag ≤ V < 20 mag) provided by the Gaia DR3 catalog (Gaia Collaboration et al. 2016, 2023). The error in the SSL intensity reflects the 10% scatter in the LORRI scattering function, and thus is systematic over all fields. We used the Gaia ESA Archive to retrieve all of the above star catalogs for each field. The SSL values for each field are given in Table 4.

The scattered galaxy light (SGL) term is the analogous scattered light contributed by bright galaxies outside the LORRI field. As with the SSL calculation, the contribution to the SGL is calculated out to an off-axis angle of 45°. Because no uniform all-sky galaxy catalog yet exists to perform this calculation using the actual positions and fluxes of known galaxies, we estimate the SGL as described in NH21. Briefly, we use the galaxy number counts from well-calibrated surveys to compute the mean surface brightness of galaxies with V < 20 and then compute the contribution to each annular bin extending out to a radius of 45°. The flux contributions in each bin are convolved with the LORRI scattering function and are then summed up to provide the final SGL estimate. The same SGL value is adopted for all fields. The surface density of bright galaxies is so low that this intensity, 0.10 ± 0.01 nW m⁻² sr⁻¹, is almost negligible. Hence, even using actual galaxy positions and brightnesses would not make a substantial difference in the final results. As with the SSL, the uncertainty in SGL is taken to be 10%.

As noted in Section 2.1, we observed eight "SCAL" calibration fields to test the SSL corrections. Figure 4 shows the predicted SSL corrections for the NCOB and SCAL fields as compared to the inferred SSL corrections estimated by subtracting all other intensity components from the total observed sky intensity. In detail, this means subtracting the field-specific IGL, SGL, ISL, and Hα components discussed in this section, as well as the DGL corrections (which will be discussed in the following section), such that SSL_inf =TotSky − IGL − ISL − Hα − DGL − S_U. Note that we also subtracted the small anomalous S_U intensity component. This is constant over all fields, and will be discussed in detail in Section 5.

Zoom In Zoom Out Reset image size

As is seen in the figure, there is excellent agreement between the predicted and inferred SSL intensities for the 16 NCOB fields. The three SCAL fields with the largest predicted SSL values, however, have relatively much smaller inferred SSL values. These fields were selected to return strong SSL backgrounds by positioning the LORRI FOV close to bright stars. The implication is that the scattered light function ≲05 radii from the LORRI field overestimates the scattered light contribution. In contrast, the NCOB fields were positioned to avoid angularly close bright stars; the fainter stars remaining close to those fields contribute little to the total SSL integral. Likewise, the five SCAL fields with the smaller predicted SSL intensities also agree well with their corresponding inferred SSL backgrounds.

Based on the SCAL data, we do not see a strong case for changing the present SSL estimation procedure, which appears to work well for the SSL intensity range experienced by our NCOB and DCAL fields. Comparing the predicted and inferred SSL values directly for the 16 NCOB fields, we find the rms relative difference between the two values to be 15%. As the error in the inferred SSL values strongly dominates in this comparison, it means that the error in the predicted SSL value must be yet smaller. The five SCAL fields close to the equality line in Figure 4 in aggregate have markedly higher SSL backgrounds than do the NCOB fields, but here the rms relative difference between the inferred and predicted SSL values rises only slightly to 18%.

ISL is the integrated light of faint Galactic stars within any field that are fainter than the LORRI photometric detection limit. Our approach is to integrate TRILEGAL models (Girardi et al. 2005, 2012) of the expected population of faint stars down to V = 30 within a 1 deg² region centered on our fields, following the procedures presented in NH21. For the present fields, the bright limit of the intensity integral (Equation (3) in NH21) is V = 19.9. The ISL errors (Table 4) are a combination of systematic and random errors due to uncertainties in the TRILEGAL model parameters and the estimated fluctuations in the star counts, respectively. As seen in Figure 5, the TRILEGAL predictions agree very well with the actual star counts from Gaia DR3 (Gaia Collaboration et al. 2016, 2023). We thus have confidence that the use of the TRILEGAL for predicting the ISL contribution down to V = 30 is a reliable approach. The ISL values for each field are given in Table 4.

Zoom In Zoom Out Reset image size

The existence of a full-sky diffuse Lyα background from the Milky Way (e.g., Gladstone et al. 2021) suggests that there might be associated hydrogen two-photon continuum (2PC; Spitzer & Greenstein 1951) at some sky locations. In NH22, we noted that the preliminary analysis of deep 130–180 nm UV spectra taken of NHTF01 suggested that a significant fraction of the continuum in the UV was due to 2PC. The corresponding optical intensity in the LORRI passband was estimated to be 0.93 ± 0.47 nW m⁻² sr⁻¹. Subsequent analysis of the NHTF01 UV spectra greatly reduced the likelihood that 2PC was needed to account for the observed continuum. The strongest argument against this hypothesis, however, was provided by the extremely low level of Hα emission seen over the LORRI field. As such, in contrast to the analysis in NH22, we conclude that 2PC does not contribute to the observed sky levels at anything more than a trivial level.

The generation of 2PC emission occurs during the recombination of ionized hydrogen. The associated Hα emission is a direct predictor of 2PC continuum intensity (Kulkarni 2022; Kulkarni & Shull 2023). The Hα levels seen in the COB survey fields (see Table 3) are all very low (<1 R, the rayleigh unit) based on observations made with the Wisconsin Hα Mapper (WHAM; Haffner et al. 2003). This in turn implies that any 2PC emission present should have intensity <0.1 nW m⁻² sr⁻¹. We conclude that 2PC emission is negligible in the present NCOB fields.

That said, since Hα emission is present in its own right in the LORRI fields, it will contribute directly to the observed total sky level, albeit at an exceedingly modest level. For one rayleigh of Hα emission, the associated intensity is 0.24 nW m⁻² sr⁻¹; the median Hα intensity for the NCOB and DCOB fields is about 2/3 of that (see Table 4).

Both the NCOB and DCAL fields were selected to have relatively low IR cirrus optical depth. Thus for dust at a given temperature we expect its emitted thermal FIR intensity to be linearly related to its surface density. The corresponding optical DGL intensity is thus proportional to the dust surface density, but it is moderated by a geometric factor, g(b), that accounts for changes in the scattering phase angle with Galactic latitude. In terms of the total sky, S_T, and the intensity components discussed in the previous section, we define the DGL+:

$$\begin{eqnarray}&&\text{DGL+}\,=\,{S}_{{\rm{T}}}-{\rm{IGL}}-{\rm{SSL}}-{\rm{SGL}}-{\rm{ISL}}-I({\rm{H}}\alpha ),\end{eqnarray} \tag{ 3 }$$

where DGL+ is the optical intensity that remains in a given field after all other known signals are subtracted from the field's total sky. The terms on the right-hand side of Equation (3) are all based on the measurements made in the optical passband. We refer to this as DGL+ because it consists of the DGL signal plus any anomalous intensity from sources not currently known or considered. For a single FIR band, a DGL+ estimator can be made by performing a least-squares linear fit using as the independent variable the FIR intensity, I(λ), corrected for the cosmic background intensity at the wavelength of the FIR intensity, CIB(λ) (Symons et al. 2023). In detail:

$$\begin{eqnarray}&&\text{DGL+}\,=\,a\,g(b)\left(I\left(\lambda )-{\rm{CIB}}(\lambda \right)\right)+c,\end{eqnarray} \tag{ 4 }$$

where g(b) is the geometric scaling factor (Jura 1979; Zemcov et al. 2017), which is defined as

$$\begin{eqnarray}&&g(b)=\displaystyle \frac{1-0.67\sqrt{\sin | b| }}{0.376},\end{eqnarray} \tag{ 5 }$$

where b is the Galactic latitude. In this definition g(b) is normalized to be unity at ∣b∣ = 60°, and thus only describes the relative effects of latitude. The absolute conversion from FIR intensity to optical DGL is provided by the coefficient, a, in Equation (4), which is determined empirically by fitting the independent variable (g(b) × (I(λ) − CIB(λ))) to the dependent variable (S_T − IGL − SSL − SGL − ISL − I(Hα)) over the combined sample of NCOB and DCAL fields. The intercept, c, provides an estimate of any anomalous intensity contribution.

We used the Planck 2018 HFI maps (Planck Collaboration et al. 2020) available from the NASA/IPAC Infrared Science Archive (IRSA) to provide the 350, 550, and 849 μm FIR intensities. As is evident in Figure 6, the LORRI field is large enough for the FIR background to vary over its extent. We thus took care to compute the average intensity within the full LORRI field of view. Those average FIR intensities for each field are tabulated in Table 3. The errors in the FIR intensities are predominantly systematic because the ∼5% calibration uncertainties in HFI are the dominant source of uncertainty (Planck Collaboration et al. 2016a).

For the low FIR emission fields selected, the CIB and dust thermal intensities are in rough parity. As such, accurate knowledge of the CIB is critical to correctly isolating the dust thermal emission needed to estimate the optical DGL foreground. The CIB has two components. The first is effectively the "intensity monopole" or full-sky average level of the background over the sky after removing dust and other Milky Way emission sources. The second is the field-to-field variations in the background due to CIB anisotropies.

For the monopole intensities at 350, 550, and 849 μm, we draw on Odegard et al. (2019), who derive 0.576 ± 0.034, 0.371 ± 0.018, and 0.149 ± 0.017 MJy sr⁻¹, respectively, from Planck HFI maps (Planck Collaboration et al. 2014, 2016a). Briefly, they use H i column density as an indicator of dust surface density, which for low dust optical depth is linearly correlated with FIR intensity. The CIB monopole is derived by extrapolating this relation to zero H i column density. As the CIB monopole intensities are common to all fields, the errors in the intensities are treated as systematic.

The CIB anisotropies for any field are provided by the Planck Collaboration et al. (2016b) analysis of the Planck HFI maps, which uses a form of power-spectrum analysis, referred to as the generalized needlet internal linear combination (GNILC) method, to separate the structure of the anisotropies from that of the thermal dust emission. These are treated as a field-dependent intensity correction to the overall CIB monopole. Specifically, the CIB-subtracted FIR intensity for a given field is

$$\begin{eqnarray}&&{I}_{c}(\lambda )=I(\lambda )-{\mathrm{CIB}}_{\mathrm{GNILC}}(\lambda )-\mathrm{CIB}\_\mathrm{MONOPOLE},\end{eqnarray} \tag{ 6 }$$

where I(λ) and CIB_GNILC(λ) are the averages over the LORRI field from the HFI and GNLIC CIB maps hosted at IRSA, and CIB_MONOPOLE is the relevant Odegard et al. (2019) value. The anisotropy corrections, CIB_GNILC(λ) + CIB_MONOPOLE, are tabulated in Table 3; the field-to-field variations of the CIB values are treated as random errors. In application, we have found that accounting for the CIB intensity anisotropies in this way removes a significant source of variance in the relations between FIR and DGL+.

Lastly, the Planck HFI did not observe at 100 μm. We thus rely on the IRIS reprocessing of the IRAS full-sky thermal-IR maps (Miville-Deschênes & Lagache 2005) for intensities at this wavelength. The CIB intensity at this wavelength is taken to be 0.78 ± 0.21 MJy sr⁻¹ (Lagache et al. 2000).

Figure 7 shows the four single-band DGL estimators based on 100, 350, 550, and 849 μm intensities, derived by fitting Equation (4) to the NCOB and DCAL DGL+ values (Equation (3)). The linear fit parameters are given in Table 5. The slope of the 100 μm estimator is nearly the same as that in the Zemcov et al. (2017) theoretical estimator but clearly has much smaller errors. The 350 μm and 550 μm estimators have even tighter trends, with the scatter in the 550 μm estimator nearly a factor of two smaller than that in the 100 μm estimator. In contrast, the 849 μm estimator offers the poorest performance of the four bands tested.

Zoom In Zoom Out Reset image size

One possible reason for larger scatter in the estimators at 100 μm versus 550 μm is field-to-field variation in the dust temperature. The average dust temperature for the NCOB and DCAL fields is 20.1 K, with a dispersion of 1.9 K. Variation in dust temperature can cause variations in the strength of the FIR intensity emitted for the same surface density of dust and the same optical light scattered. At 20 K, 100 μm falls on the short-wavelength side of the peak of the blackbody spectrum, and thus 100 μm emission is more sensitive to small temperature changes than is the intensity at 550 μm, which falls on the long-wavelength side of the peak. We thus investigated two-band DGL estimators, finding that using the 350 μm and 550 μm intensities in combination gives the best performance, returning smaller scatter than the 550 μm single-band estimator.

In detail, the two-band estimator fits the DGL+ value for any field as

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{DGL}+ & = & {c}_{1}+g(b)\,\left[{c}_{2}{I}_{c}(550\,\mu {\rm{m}})+{c}_{3}\left(\displaystyle \frac{{I}_{c}(350\,\mu {\rm{m}})}{{I}_{c}(550\,\mu {\rm{m}})}\,\right.\right.\\ & & \left.\left.-\,\left\langle \displaystyle \frac{{I}_{c}(350\,\mu {\rm{m}})}{{I}_{c}(550\,\mu {\rm{m}})}\right\rangle \right)\right],\end{array}\end{eqnarray} \tag{ 7 }$$

where the subscripts on the intensities indicate that the CIB intensity has been subtracted from them. The mean 350–550 μm intensity ratio is subtracted from ratio for each field to strongly reduce covariance of the intensity-ratio term with the overall intercept term. The fit of this estimator is shown in Figure 8, with coefficient values c₁ = 2.60, c₂ = 48.01, and c₃ = 0.96, with rms residuals of 1.39 nW m⁻² sr⁻¹. The mean 350–550 μm intensity ratio for our sample is 3.66.

Zoom In Zoom Out Reset image size

To use the two-band estimator to predict just the DGL value, we subtract the c₁ coefficient, as the predicted DGL must go to zero when the FIR intensity goes to zero. Specifically, our DGL predictor is

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{DGL}(\mathrm{nW}\,{{\rm{m}}}^{-2}\,{\mathrm{sr}}^{-1})\\ =\,g(b)\,\left[48.01\,{I}_{c}(550\,\mu {\rm{m}})\,+\,0.96\,\left(\displaystyle \frac{{I}_{c}(350\,\mu {\rm{m}})}{{I}_{c}(550\,\mu {\rm{m}})}\,-\,3.66\right)\right],\end{array}\end{eqnarray} \tag{ 8 }$$

where the FIR intensities are in units of MJy sr⁻¹ and the predicted DGL value is in units of nW m⁻² sr⁻¹. The estimated DGL values for all NCOB and DCAL fields are given in Table 4. The errors in Table 4 are the quadrature sums of the random and systematic errors. Those individual random and systematic errors are given in Table 6. The errors for DGL are the rms values obtained from the Monte Carlo analysis described in Section 5.

To visualize the relative importance of all sky components we represent the results as a stacked bar chart for each field in Figure 9. While the summed intensity from every field is less than the total sky level, the dominance of systematic errors in most of the intensity components means the uncertainty in the combined data set will not be diminished by the $1/\sqrt{N}$ factor. We present our approach for determining the proper propagation of errors for the combined suite of data in the next section.

Zoom In Zoom Out Reset image size

Having identified all the sources of foreground optical emission known to us in the two previous sections, the task is then to recover an optimal estimate of the COB intensity with accurate errors that reflect the appropriate random and systematic uncertainties (see Table 6), as well as any covariances among the parameters used to estimate the COB intensity. We will use a Monte Carlo approach that randomly generates complete realizations of the COB observations as based on our error model. As noted at the start of Section 3, we do this in two steps. With the exception of the DGL component, our knowledge of all other foregrounds, as well as the total sky level itself, comes from mutually independent information. As such, we start with Equation (3) to estimate DGL+, the DGL intensity plus any anomalous intensity component, which is not affected by covariance among the terms that are subtracted off to isolate it. The second step is to then use the FIR background intensity to isolate the DGL component itself. As this uses the observations to develop a self-calibrated DGL estimator, this step does account for covariance between the observational parameters.

For the first step then, the Monte Carlo routine generates 10,000 realizations of the observations of total sky intensity and the non-DGL foreground components (S_T, IGL, SSL, SGL, ISL, BIGL, and I(H α)) for each field. We also generate random realizations of the independent variables (FIR intensities and cosmic IR backgrounds) that will be used to estimate the DGL. In detail, for any observed component or observational input, D_obs(j), in field j (1 ≤ j ≤ 24; 16 NCOB + 8 DCAL fields), the Monte Carlo routine generates a set of simulated values for 1 ≤ i ≤ 10⁴:

$$\begin{eqnarray}&&{D}_{\mathrm{sim}}(i,j)={D}_{\mathrm{obs}}(j)\,+\,{\sigma }_{\mathrm{RAN}}(j){G}_{{\rm{D}}}(i,j)\,+\,{\sigma }_{\mathrm{SYS}}(j){F}_{{\rm{D}}}(i),\end{eqnarray} \tag{ 9 }$$

where G_D(i, j) and F_D(i) are sequences of Gaussian random variables with zero mean and unit variance. Note that F_D(i) scales the systematic error term σ_SYS(j) in any field, and is thus the same for all fields for a given realization i. For many parameters, σ_SYS will also be constant over all fields, but relative systematic errors may vary from field to field in step with the given parameter (as occurs in the case of SSL, for example). In contrast, G_D(i, j) scales the random errors, which vary freely from field to field in any given realization. Lastly, we note that while Equation (9) is fully general in allowing for both systematic and random errors for any component, in practice the systematic errors strongly dominate the total error budget, and for many components the random error term is negligible.

A two-band DGL estimator of the form given by Equation (7) is generated for every Monte Carlo realization from its set of 24 DGL+ values produced, which is used to predict the DGL backgrounds for that particular realization. As is described by Equation (6), the independent variables used by the estimator are the FIR background intensities at 350 μm and 550 μ m, with the estimated CIB backgrounds subtracted. Fitting the DGL+ values to the FIR background intensities generates the c₁, c₂, and c₃ coefficients in Equation (7). The best-fit DGL-alone estimates for each field are then provided by the c₂ and c₃ coefficients on the assumption that DGL is zero for zero FIR background intensity. This directly highlights the strong covariance between the assumed CIB backgrounds and the DGL predictions noted in Section 4.2. This covariance is included in the Monte Carlo routine, however, simply by including errors in all the FIR parameters as part of the analysis. The coefficients of the DGL estimator are also covariant, but again by having the Monte Carlo routine sample the error distributions of both the independent and dependent variables used by the DGL estimator, this behavior is included implicitly in the analysis as well. Further, because the DGL+ intensities for any field reflect the errors of all components used to generate them, the Monte Carlo analysis also incorporates the full suite of covariances associated with deriving the DGL estimator from the COB observations themselves.

Lastly, we note that while we presented a mean two-band DGL estimator in Equation (8), we actually do not use it as such for the final analysis. Again, we generate the DGL estimator anew for each Monte Carlo realization. Our interest is really just in the distribution of DGL estimates for any given field generated by the full ensemble of Monte Carlo runs. These distributions, in fact, rather than any simple error calculation, provide the appropriate uncertainties in the DGL values needed to derive the errors in the final COB and S_U intensities.

With the DGL component derived for all fields in all 10,000 Monte Carlo runs, it is simple to derive the COB intensity and any residual anomalous intensity background, S_U, for each realization. The optimal estimate of both parameters for the whole survey will be derived by first generating the COB and S_U intensity distributions for each field, then combining them into a full ensemble distribution weighted by the inverse variance of the individual-field distributions.

The anomalous component, S_U, is just the residual intensity remaining after all known sources of light are subtracted from the data. Specifically,

$$\begin{eqnarray}&&{S}_{{\rm{U}}}={S}_{{\rm{T}}}-{\rm{IGL}}-{\rm{SSL}}-{\rm{SGL}}-{\rm{ISL}}-I({\rm{H}}\alpha )-{\rm{DGL}}.\end{eqnarray} \tag{ 10 }$$

We compute S_U for each field using the 10,000 realizations of each parameter on the right-hand side of Equation (10). We then generate the cumulative distribution function of S_U for each field and find its 68.3% confidence limits. We adopt the 1σ value on S_U to be half the difference between high and low 68.3% confidence limits for that specific field. We then generate a summed distribution of all the S_U values for all fields, weighting the S_U values for each field by the inverse square of that field's corresponding 1σ value. The value of S_U for the full sample is then taken to be the mean value of that weighted sum distribution, and the error in the full sample S_U is taken as the half width of the 68.3% confidence limit range of that weighted sum distribution. The resulting summed distribution of S_U values for all 16 COB survey fields is shown in Figure 10.

Zoom In Zoom Out Reset image size

To get the optimal full-survey COB intensity, we repeat the above procedure but apply it to the following expression:

$$\begin{eqnarray}\begin{array}{rcl}{\rm{COB}} & = & {S}_{{\rm{T}}}-{\rm{SSL}}-{\rm{SGL}}-{\rm{ISL}}-I({\rm{H}}\alpha )-{\rm{DGL}}+{\rm{BIGL}}\\ & = & {S}_{{\rm{U}}}+{\rm{IGL}}+{\rm{BIGL}}.\end{array}\end{eqnarray} \tag{ 11 }$$

We remind the reader that we need to add back the integrated light of bright galaxies (BIGL) prior to computing the COB since all galaxies brighter than V = 19.9 had been masked out prior to measuring the total sky level (see Section 3.1 for details). The resulting summed distribution of COB values for all 16 survey fields is shown in Figure 11.

Zoom In Zoom Out Reset image size

Our derived COB and S_U values are

$$\begin{eqnarray}\begin{array}{rcl}{\rm{COB}} & = & 11.16\pm 1.65\,(1.47\,{\rm{sys}},\,0.75\,{\rm{ran}})\,{\rm{nW}}\,{{\rm{m}}}^{-2}\,{{\rm{sr}}}^{-1}\\ {S}_{{\rm{U}}} & = & 2.99\pm 2.03\,(1.75\,{\rm{sys}},\,1.03\,{\rm{ran}})\,{\rm{nW}}\,{{\rm{m}}}^{-2}\,{{\rm{sr}}}^{-1}.\end{array}\end{eqnarray} \tag{ 12 }$$

The COB is detected at a significance of 6.77σ but the anomalous (residual) sky intensity is only present at the 1.47σ level. In other words, S_U is likely consistent with zero.

To determine whether any one of our NCOB or DCAL fields has an outsized influence on either the COB or S_U estimate, we perform a series of jackknife tests, where for each test we exclude just one field from the Monte Carlo simulations, doing this in turn for each field in the full set of 16 NCOB and 8 DCAL fields. The outcomes of these tests are shown in Figure 12. The x-axis labels show the name of the field that was excluded in each 10,000 run simulation. The red and blue data points show, respectively, the resulting COB and S_U values derived when that field was excluded from the analysis. The standard deviation in the COB or S_U for these 24 jackknife tests is just 0.12 nW m⁻² sr⁻¹. We conclude that the COB and S_U results are resilient against the removal of any one of the NCOB or DCAL fields from the analysis.

Zoom In Zoom Out Reset image size

As noted in Section 2.1, we obtained observations of four fields at low ecliptic latitude to verify the assumption that there is no significant ZL emission at heliocentric distances beyond 50 au. We perform this test by computing the COB and S_U intensity values for the four IPD fields and look for any significant deviations from the main survey results listed in Equation (12). The results are given in Table 7. Because the IPDF observations were not part of our rigorous survey selection process, the errors on the COB and S_U values for the IPD fields are estimated analytically rather than with the Monte Carlo procedure but should be approximately similar to what the Monte Carlo process would yield. Columns (3) and (5) in Table 7 give the difference between the COB and S_U given in Equation (12) and the corresponding values for each IPD field, respectively. We find no significant differences between the COB and S_U values obtained for the IPD fields and those derived from the COB science fields, supporting the assertion that ZL is negligible at the distances where NH made the COB observations presented here. The consistency between the COB science observations and the IPDF observations also suggests that there are no major zodiacal subtraction residuals in the Planck HFI maps.

With 16 fields spread over the sky, we can, in principle, look for anisotropy in the COB. Of the 16 NCOB fields, 13 are in the south Galactic hemisphere (SGH) and three are in the north Galactic hemisphere (NGH). We derive the COB for these two subsets to see whether there are any significant offsets between them. We use the 10,000 realizations generated and separate out the north and south Galactic subsamples from them. The COB and S_U values are then calculated. The results are shown in Table 8. We also plot the individual S_U values as a function of ecliptic and Galactic latitude in Figure 13. We do not see any significant difference between the COB or S_U values derived for the NGH versus SGH subsamples. Nor do we seen any significant trend with ecliptic latitude with the current sample. We note that with only three NCOB fields in the NGH and only three DCAL fields in the SGH, we use the full sample to compute the two-band DGL estimator in this isotropy analysis. A larger NGH NCOB sample and more DCAL observations in the SGH would be required to do a more in-depth exploration of any anisotropy in the COB signal.

Zoom In Zoom Out Reset image size