New Horizons Observations of the Cosmic Optical Background

Scattered sunlight contributes significantly to the LORRI background at solar elongation angles (SEAs) < 90°. Many DKBO observations have been taken in such "bright" conditions at Sun angles well interior to this limit, but there is no calibration methodology that would cleanly isolate the astronomical sky from their considerably stronger scattered-light backgrounds.

For SEAs > 90°, LORRI is in the shadow of the spacecraft, so sunlight can no longer enter its aperture directly. However, spacecraft components extending out from the LORRI-side panel, such as the star trackers, may still be in sunlight, depending on the particular geometry of the observation, and thus may scatter sunlight indirectly into the LORRI aperture at a shallow angle. This problem becomes less important as the SEA increases past 90°, but because the star trackers are mounted close to the bottom edge of the LORRI-side panel (see Figure 2), surprisingly large SEAs are required, if the sunlight is coming from below the spacecraft, to have them completely shadowed. At the same time, because the New Horizons trajectory out of the solar system is directed toward the Galactic center, SEAs approaching 180° will always have low Galactic latitudes. This means that using LORRI to measure the COB at high Galactic latitudes requires evaluating the effects of scattered sunlight for the likely SEAs, even when the LORRI aperture itself is shadowed.

To explore the importance of indirect sunlight for our selected fields, we used ray tracing on a model of the spacecraft generated by one of the authors (D.D.D.). The model was developed for use in artistic renditions and is faithful to the geometric configuration of the various spacecraft components while making plausible assumptions about their surface properties and reflectivities. We thus consider the renderings useful for an overall investigation into whether scattered sunlight might be important at the order-of-magnitude level, while being wary of using them for highly precise corrections.

With these caveats, we generated illuminated renditions of the full spacecraft appropriate to the spacecraft attitudes used to observe the fields discussed at the end of this section. ²² The LORRI-side panels are shown in Figure 2. As can be seen, the two star trackers do remain slightly illuminated for all fields, creating glints that can be observed from the entrance of the LORRI aperture.

Evaluation of the effect of scattered sunlight was done by filling the aperture with a white screen and measuring its average illumination from the star tracker glints. To calibrate this measure we then illuminated the screen with direct sunlight at the elevation angle and direction of the star tracker glints as seen at the aperture entrance. From DKBO images obtained at 70° < SEAs < 90°, we knew the sky level generated by direct sunlight as a function of its incidence angle with respect to the LORRI aperture. This sky level was then simply scaled by the relative glint / direct sunlight ratio to provide the estimated scattered-sunlight contribution to the total-sky levels measured.

For example, the n3c61f field has the smallest SEA and incurs the strongest scattered-sunlight effects. From the ray-tracing tests appropriate to the spacecraft attitude used to observe this field, we observed a glint that had at most 1.8 × 10⁻³ of the solar flux and a 9° elevation as observed from the LORRI aperture. ²³ At SEA = 81°, direct sunlight produced a 21 DN background in 30 s. Multiplying this by the relative flux of the glint, we estimated that the glint produced a scattered-sunlight background of 0.04 DN in 30 s. This is ∼20% of the final dCOB signal that we recovered by the analysis to be described later. The ZL 138c and ZL 138d fields have over an order of magnitude smaller scattered-sunlight components, and the remaining four fields have no detectable scattered sunlight entering the camera at upper limits two orders of magnitude smaller than that in the n3c61f field. Because only one field has a scattered-sunlight component at even relatively small significance, we have chosen not to correct the average sky for this effect.

2.1.1. Glittering Ammonia Crystals?

The need for a long run of images in a data set is to counter a recently discovered background effect within the camera (Weaver et al. 2020). The effect manifests as elevated bias and sky levels, which are most pronounced in the first exposure in a sequence but require a ∼150 s interval (Figure 4) to decay fully to a constant level. This effect appears to be associated with the power-on activation of the camera, which typically occurs 30 s before the start of an imaging sequence. Notably, it is not seen at the start of sequences that occur well after the initial sequence but still within the same operational interval over which LORRI is continuously powered. With sequences of 30 s exposures, we chose to discard the first five exposures following the activation of the camera to avoid the effects of this phenomenon.

Zoom In Zoom Out Reset image size

We identified seven LORRI data sets that satisfied all the constraints. The sequence parameters are given in Table 1. The field locations with respect to a Galactic extinction map derived from the Infrared Astronomical Satellite (IRAS) 100 μm all-sky survey (Schlegel et al. 1998) are shown in Figure 5. Three of the seven data sets were provided by DKBO observations. The DKBO data sets comprise several imaging sequences of three KBOs, each taken over an interval of 14–48 hr. The LORRI pointing changed over each set as needed to track the KBO but varied by only about half a LORRI field. The background 100 μm flux did not vary significantly over each set as the pointing changed (as will be discussed in detail in the analysis section). Each set can thus be combined into a single measurement.

Zoom In Zoom Out Reset image size

Table 1. Sky Fields and Imaging Sequences

			Solar				I100	Date	r
Program	R.A.	Decl.	Elong.	b	β	E(B − V)	MJy sr−1	(UT)	(au)	${N}_{\mathrm{seq}}/{N}_{\mathrm{im}}$
OE 394	1.7790	−17.7780	110.1	−76.1	−17.0	0.029	1.02	2018-08-20	42.2	6/78
OJ 394	348.0611	−41.6360	125.0	−65.1	−33.2	0.010	0.53	2018-08-22	42.2	9/120
n3c61f	33.4069	−50.7529	96.3	−61.8	−58.1	0.017	0.89	2018-09-01	45.2	15/90
ZL 138a	358.2428	−0.5183	108.3	−59.8	0.2	0.025	0.73	2019-09-04	45.3	1/8
ZL 138b	0.8066	0.2915	105.6	−60.2	−0.1	0.031	0.97	2019-09-04	45.3	1/8
ZL 138c	224.9875	36.2331	98.1	61.4	50.3	0.012	0.68	2019-09-04	45.3	1/8
ZL 138d	226.4865	35.2979	99.6	60.3	50.0	0.011	0.70	2019-09-04	45.3	1/8

Notes. Column (1) partial program name ("ZL" stands for "zodiacal light"), (2) sky field R.A. (2000), (3) decl. (2000), (4) solar elongation, (5) Galactic latitude, (6) ecliptic latitude, (7) reddening from Schlafly & Finkbeiner (2011), (8) average 100 μm flux for the field from the Miville-Deschênes & Lagache (2005) maps, with correction for residual ZL and the 0.78 MJy sr⁻¹ cosmic IR background (Fixsen et al. 1996; Puget et al. 1996) subtracted (see Section 4.4), (9) UT starting date, (10) New Horizons distance from the Sun, and (11) number of sequences / total number of images. All angles are in degrees. For the first three fields, the coordinates given are an average over all sequences. The E(B − V) and 100 μm flux values were obtained from the IRSA archive dust tool at https://irsa.ipac.caltech.edu/applications/DUST/.

Download table as:  ASCIITypeset image

The remaining four data sets are single-pointing image sets generated by a ZL program designed to provide a quick test for obvious dust within the plane of the cold classical Kuiper Belt. The program comprised four fields at the same high Galactic latitude (∣b∣ ∼ 60°), but with two fields at an ecliptic latitude β ∼ 0° and two at β ∼ 50°. The target fields that were chosen also have low surface densities of bright stars. A sequence of eight consecutive 30 s LORRI exposures were obtained of each field. The reduction of these sequences led to the discovery of the initially high instrument background shown in Figure 4.

The sky levels are less than 1 DN in the nominal 30 s exposures. Reduction of the images thus requires attention to a number of subtle effects that are only important at this level. Rather than using calibrated ("Level 2") images produced by the standard LORRI pipeline operated by the New Horizons project, we designed a custom reduction of the raw ("Level 1") images to optimize accurate recovery of the faint sky signal.

3.1.1. Bias Level Determination

3.1.2. Measurement of the LORRI Dark Current

The bias column also integrates any dark current present in the detector over the duration of the exposure. As such, any constant dark level over the field will be treated as part of the bias level and removed with it. As is discussed in Weaver et al. (2020), the dark current in a LORRI 1 × 1 mode pixel is estimated to be 4 × 10⁻³ e⁻ s⁻¹ pixel⁻¹ or ∼0.1 DN/pixel in a 30 s exposure in 4 × 4 mode, based on the manufacturer's specifications and the CCD operating temperature. A modest contribution to the background at this level should be well characterized by the bias column.

Higher dark levels, however, place greater demands on the accuracy of the underlying assumption that the bias column truly witnesses the average dark current appropriate for the active imaging area of the LORRI CCD. Unfortunately, without a dark shutter we cannot obtain true dark calibration exposures over the full field in flight; however, we have been able to obtain calibration data that demonstrates that the dark level measured by the unilluminated bias column is indeed close to the expected level.

In the ZeroDark65 LORRI sequence, which operated on the spacecraft in 2020 June, we obtained 32 pairs of 65 s blank-sky ²⁴ exposures immediately followed by a 0 s image obtained 1 s later. While the electronic bias level can drift over a long imaging sequence, we have found it to be stable by a fraction of a DN between consecutive exposures. If the electronic bias level were truly constant, the difference between the measured bias plus dark levels recorded by the unilluminated columns in the pair of images would provide a direct measure of the dark level. With 32 such pairs the measurement noise and any random bias variations should average out. Indeed, the bias level drifted systematically by no more than 0.5 DN over the full duration of the ZeroDark65 sequence.

Figure 6 shows the dark current measures derived by differencing the thirty-two 65 s and 0 s image pairs. The average level measured is 0.334 ± 0.039 DN in 65 s, or 0.154 ± 0.018 DN in 30 s. Referenced to the LORRI 1 × 1 mode, this is 6.2 × 10⁻³ e⁻ s⁻¹ pixel⁻¹,  which is within the range of the CCD manufacturer's specified performance. This modest dark level should be completely corrected by the measured bias level subtraction described in the previous section.

Zoom In Zoom Out Reset image size

3.1.3. Correction for Bias Structure

A variety of coherent and repeatable patterns are evident at low signal levels in LORRI images. The standard LORRI pipeline subtracts a "super-bias" image constructed from the average of a large set of bias frames obtained shortly after the New Horizons spacecraft was launched but before the LORRI aperture was opened. We instead used a revised super-bias image based on extremely short in-flight images obtained with the camera, processed to correct only large-scale features in the bias frame.

LORRI frames exhibit a low-level "jail-bar" pattern in which the average level of the even-versus-odd columns differs by 0.5 DN (Weaver et al. 2020). This pattern is not treated in the standard pipeline and can strongly affect determination of the sky level. Since the bias column has odd parity (it is read out as column 257, following the 256 illuminated columns), the bias level determined from it is strictly valid only for the odd columns. A key aspect of the pattern is that the even columns vary from sequence to sequence in being either a −0.5 or a +0.5 DN offset from the odd columns. The sign of the offset appears to stay constant over a given LORRI power-on interval but changes randomly from different operational cycles (Figure 7). Without correcting for this effect, the measured sky level will vary between being −0.25 DN too low and being +0.25 DN too high. This effect is evident as systematic differences of 0.5 DN between different sequences obtained of the same field, as is shown in Figure 8. An error of this size is close to the full value of the sky signal itself.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The solution to the jail-bar offset is simple. Given the faint background levels in the present image sets, it is easy to determine the sign of the 0.5 DN offset of the even columns from the odd columns and thus apply the appropriate correction. Figure 8 shows the same data set before and after the jail-bar correction was applied. The agreement of the average sky level between visits with opposite jail-bar parities is now excellent.

Lastly, the bias level appears to show random low-level variations over the duration of the image readout, which are evident as horizontal streaking in the image backgrounds. The bias column forces this pattern to have zero mean at the high–column number margin of the images, but the bias may also exhibit a slow low-amplitude drift along the CCD rows, such that the overall mean of the bias background alone may have a small nonzero value. At present this pattern is not corrected, and it may contribute to some of the random exposure-to-exposure variations in the sky level, such as those visible in Figure 8.

To validate these reduction steps, we tested them on two 0 s LORRI 4 × 4 exposures obtained on 2019 July 13 for routine performance monitoring of the camera. Using the histogram estimation methodology to be described in Section 3.2.2, we measured the average signal level to be 0.02 ± 0.06 DN, demonstrating that we measured zero signal in data expected to have a sky level of zero.

3.1.4. Charge Smearing Correction and Flat-fielding

After the average bias level is subtracted from the image and the bias structure corrections described above are completed, a "charge smearing" correction is applied (Weaver et al. 2020). In brief, because LORRI does not have a camera shutter, light from an astronomical source will still deposit charge in the LORRI CCD as it is being read out and as the CCD is clocked in advance of the exposure to erase or scrub out charge deposited prior to the start of the exposure. The exposure is made by simply stopping the charge clocking for its duration. The effect of charge being deposited both before and after the nominal exposure interval is to generate charge trails in the image rows above and below a bright source. The distribution of smeared charge is corrected prior to flat-fielding.

The amplitude of the smeared charge in any given pixel in the same column as a bright source can be estimated from the flux of the source, scaled to the time required to clock out a single row of the image. For the LORRI 4 × 4 mode this is only 0.047 ms. In the present case, for a bright star that just saturates at 4095 DN in 30 s, the smeared charge is only 0.005 DN/pixel, which is only ∼1% of the typical sky levels that we are concerned with. Even so, we still applied smear correction using the new algorithm discussed in Weaver et al. (2020). We note that we flagged cosmic-ray hits in any given image in advance of the smear correction, as these are instantaneous sources that will not be smeared. Including cosmic-ray hits in the correction will induce a small negative bias to other pixels in the affected column.

The last image reduction step is flat-fielding sensitivity correction, which is done with the standard pipeline procedures and products. This is a minor correction. The response of the CCD is highly uniform over its area, with rms sensitivity variations of only 0.9%. The final reduced images are shown in Figure 9.

Zoom In Zoom Out Reset image size

3.2.1. Exclusion of Non-sky Flux Sources

3.2.2. Determining the Peak Location of the Intensity Histogram

The sky level for any given image is determined from the peak of its pixel intensity histogram. The histograms were generated in a way that best preserved the information content of the images (Figure 10). For pixels sampling the low sky levels most of the image reduction steps only slightly altered the raw integral DN values generated by the camera; thus care had to be taken in how the signal was binned to generate the histogram. The most important consideration was to preserve the 0.5 DN offset correction applied to the jail-bar pattern, and the fractional DN portion of the bias level measurement. The procedure adopted was to use 0.5 DN wide bins for the histogram and to phase the centers of the bins to retain the fractional part of the bias level subtracted from the initially integer DN values. Lastly, the average intensity of all the pixels within a given bin was used to define the center of the bin, rather than the simple midpoint of the bin interval.

Zoom In Zoom Out Reset image size

We explored three algorithms for measuring the precise intensity of the histogram peak. All approaches used the bin with the maximum occupation number and bins on either side of it that had occupation numbers falling just below half of the peak value. The first procedure was to fit a parabola to the occupation numbers as a function of bin average intensity. The second was to fit a parabola to the logarithm of the occupation numbers, which is equivalent to fitting a Gaussian to the bins. Lastly, a simple centroid of the peak bins was calculated. All three measures agreed at the 0.01 DN level. We adopted the Gaussian fit measures as they produced the smallest scatter in measurements over a large sample.

The sky levels for each field are given in Table 2. They are averages of the individual levels of all images available for a given field, excluding the first five images obtained following the powering up of LORRI. The errors are the statistical errors of the mean. The typical error in any single image is 0.15 DN or 7.5 nW m⁻² sr⁻¹  but does appear to vary somewhat between sequences. Conversion to flux/solid-angle values is provided by Equation (8) in Weaver et al. (2020), assuming the RSOLAR zero-point defined in that paper. Direct conversion of a sky level, S, in DN/(LORRI 4 × 4 pixel) in 30 s, assuming a pivot wavelength of 0.608 μm, is

$$\begin{eqnarray}&&\lambda {I}_{\lambda }=49.5\,S\,\mathrm{nW}\,{{\rm{m}}}^{-2}\,{\mathrm{sr}}^{-1}.\end{eqnarray} \tag{ 1 }$$

Table 2. Measured Total-sky Levels

Program	N	DN	e−	nW m−2 sr−1
OE 394	63	0.758 ± 0.025	14.7 ± 0.49	37.52 ± 1.24
OJ 394	105	0.613 ± 0.012	11.9 ± 0.23	30.34 ± 0.59
n3c61f	15	0.803 ± 0.034	15.6 ± 0.66	39.75 ± 1.68
ZL 138a	3	0.737 ± 0.053	14.3 ± 1.03	36.48 ± 2.62
ZL 138b	3	0.804 ± 0.044	15.6 ± 0.85	39.80 ± 2.18
ZL 138c	3	0.641 ± 0.045	12.4 ± 0.87	31.73 ± 2.23
ZL 138d	3	0.721 ± 0.028	14.0 ± 0.54	35.69 ± 1.39

Note. Column (1) partial program name, (2) number of images used, (3) average sky-level DN/pixel in 30 s, (4) sky in electrons/pixel in 30 s, and (5) sky in intensity units.

Download table as:  ASCIITypeset image

Background light from stars within the fields below the detection limit of the LORRI images was estimated using version 1.6 of the TRILEGAL Milky Way model (Girardi et al. 2005, 2012). For each LORRI field, we generated a one square degree simulation centered on the boresight coordinates given in Table 1. The simulations were performed using a limiting magnitude of J = 30 and the default parameters for the four main Galactic components: the thin disk, thick disk, halo, and bulge. Given the high Galactic latitudes of the observations used here, the simulations only include stars from the disk and the halo. As the LORRI field covers 029 × 029, running a one square degree simulation effectively generates ∼12 images worth of stars for each pointing, which gives us ample mitigation against sample variance. The variation seen in the simulated star counts when running TRILEGAL for a given sky field multiple times using identical model parameters is less than 0.7%.

We derived the number counts of stars as a function of magnitude for the UBVRIZJ passbands to encompass the LORRI sensitivity range. We applied the following (Vega–AB) magnitude conversions, as needed: −0.7194 (U), +0.123 (B), −0.017 (V), −0.211 (R), −0.450 (I), −0.573 (Z), and −0.940 (J). The star number counts in the V band are shown in Figure 11. We also show in this figure the actual star counts from the Gaia DR2 catalog (Gaia Collaboration et al. 2016, 2018) in the one square degree regions centered on each LORRI dark-sky pointing. We used the transformation provided in the Gaia DR2 photometric validation paper (Evans et al. 2018) to convert the Gaia photometry to the V band:

$$\begin{eqnarray}\begin{array}{rcl}G-V & = & -0.01746-0.006860({G}_{\mathrm{BP}}-{G}_{\mathrm{RP}})\\ & & -0.1732{\left({G}_{\mathrm{BP}}-{G}_{\mathrm{RP}}\right)}^{2}.\end{array}\end{eqnarray} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

The agreement between the Gaia observations and the TRILEGAL simulations is excellent. This agreement is, in part, due to the fact that the TRILEGAL algorithm has been calibrated against actual star counts, providing confidence that the extrapolation to fainter flux levels is likely to be quite reliable. The average difference between the simulated star counts and the observed Gaia star counts for our target fields is 5.5% over the range 10 ≥ V > 19 mag. The variation in star counts from field to field is primarily due to variation in the Galactic coordinates of the fields, as demonstrated in Figure 12.

Zoom In Zoom Out Reset image size

We generated an estimate of the total V-band surface brightness of undetected stars in the LORRI fields by integrating the simulated star counts from 30 mag to the detection limit of the LORRI images using the expression

$$\begin{eqnarray}&&\begin{array}{l}{\mu }_{* }(\mathrm{mag}\ {\mathrm{arcsec}}^{-2})\\ =\,-2.5\ {\mathrm{log}}_{10}\left[\displaystyle \frac{1\ {\deg }^{2}}{({3600}^{2}\ {\mathrm{arcsec}}^{2})}{\displaystyle \int }_{{m}_{\mathrm{Faint}}}^{{m}_{\mathrm{Lim}}}N(m)\ {10}^{-0.4{\rm{m}}}{dm}\right],\end{array}\end{eqnarray} \tag{ 3 }$$

where N(m) is the differential number of stars per unit magnitude per square degree predicted by the TRILEGAL model, dm is the magnitude interval used in the integration (0.10 mag), m_Faint is 30 mag, and m_Lim is the detection limit (∼19.1 mag in V for the typical LORRI image). We used the N(m) derived from the simulated star catalogs for each of the seven fields. The faint-object surface brightness results for each New Horizons dark-sky field are presented in Table 3. The errors for the surface brightnesses of the undetected faint stars were derived from the 1σ $\sqrt{N}$ uncertainties in the simulated TRILEGAL number counts. The surface brightness levels in Table 3 define the normalization for our spectral energy distributions (SEDs) that were used to derive the estimated contribution of undetected sources to the total LORRI signal.

The SED of the stellar population varies with the limiting magnitude due to the increasing fraction of late-type stars and low-mass stars at fainter magnitudes. Given the very broad passband of LORRI (Figure 3), it is essential that we model the change in average stellar SED as a function of magnitude in order to compute a robust estimate of the signal from undetected stars. We estimated the dependence of the SED of the stellar population on apparent magnitude by computing the mean UBVRIJHK photometry for stars in the TRILEGAL simulations in eleven 1 mag wide bins from V = 19 to V = 30. In each bin, we used the mean photometry to generate an average SED, F_λ (λ), for that bin, normalized to the V-band surface brightnesses given in Table 3. The field-to-field variation in the mean stellar colors as a function of magnitude in the TRILEGAL simulations of the New Horizons target regions is negligible (<0.04 mag).

For each magnitude bin, we then computed the expected LORRI count rate per unit solid angle, C_L, for the undetected faint stars by integrating the appropriate SED of the sources over the LORRI response function, where

$$\begin{eqnarray}&&{C}_{{\rm{L}}}(\mathrm{DN}\ {{\rm{s}}}^{-1}\ {\mathrm{sr}}^{-1})=\displaystyle \frac{1}{{{\rm{\Omega }}}_{\mathrm{pixel}}}\ {\int }_{{\lambda }_{1}}^{{\lambda }_{2}}R(\lambda ){F}_{\lambda }(\lambda )d\lambda ,\end{eqnarray} \tag{ 4 }$$

where R(λ) is the LORRI absolute responsivity function (shown in Figure 3), F_λ (λ) is the SED expressed as a flux density, and Ω_pixel is the solid angle subtended by a single LORRI 4 × 4 binned pixel (3.91264 × 10⁻¹⁰ sr). The flux density, F_λ (λ), for each magnitude bin was generated by fitting a cubic spline to the derived mean UBVRIJHK magnitudes (in the AB system) using a wavelength step size, dλ, of 10 nm and a LORRI 4 × 4 pixel area of 16.6464 arcsec² and then converting to the appropriate cgs units (erg cm⁻² s⁻¹ nm⁻¹). The integration was performed over the full LORRI sensitivity range from 0.30 to 1.00 μm. The LORRI counts in a 30 s exposure were then computed and converted to an intensity expressed in units of nW m⁻² sr⁻¹ via Equation (1). The total signal over the full magnitude range 19 < V ≤ 30 is then just the number count–weighted sum of the C_L values in each magnitude bin. The intensities for the faint field stars derived in this way are listed in Table 4. The turnover in star counts at magnitudes fainter than V = 24 mag (see Figures 11 and 13) means that the derived sky intensity from undetected faint stars did not change significantly so long as the faint limit used in Equation (3) was V > 24 mag, at which point the integral reached >98% of the value obtained with our chosen integration limit of V = 30 mag.

The errors in the derived intensities from the TRILEGAL simulations are a combination of statistical $\sqrt{N}$ errors in the star counts and systematic variations that depend on model-specific parameters. As demonstrated in both Figures 11 and 12, the default TRILEGAL parameters reproduce the observed Gaia star counts at V ≤ 20 extremely well. But to allow for plausible model uncertainty, we estimated the systematic changes on the derived sky intensity due to changes in the model parameters that shifted the derived model star counts by ±2σ from the observed star counts. Shifts of this amplitude are achieved by systematically changing the disk scale heights, halo effective radii, and halo oblateness by ±10%. Allowing the model parameters to vary over a larger range than this would result in simulations that were in strong disagreement with the observations. Over the range 19 < V ≤ 30, the TRILEGAL models predicted that, for the New Horizons fields, 20%, 24%, and 56% of the stars, respectively, are part of the thin disk, the thick disk, and the halo. Hence, the parameters for all three components are key to setting the predicted star counts at high Galactic latitudes, although at magnitudes fainter than V = 25, halo stars account for ∼68% of the population.

In sum, the statistical error in the faint-star sky contribution accounts for less than 10% of the uncertainty in this signal; the error in the sky component due to undetected faint stars is dominated by the systematic uncertainties associated with variation in the TRILEGAL models. The combined fractional error from both the statistical and systematic uncertainties is 14%–17% of the derived faint-star signal, and furthermore, the signal due to those undetected faint stars contributes, on average, just ∼7% of the total-sky signal.

We estimated the integrated galaxy light (IGL) in a manner similar to that used to compute the contribution of faint stars below the LORRI detection limit. We used a compilation of well-measured galaxy number counts in the UBVRIz passbands. We took galaxy number counts from various literature sources with incompleteness corrections applied to the raw number counts. Specifically, we used the observed galaxy number counts from the following works: Hogg et al. (1997), Yasuda et al. (2001), Grazian et al. (2009), Laigle et al. (2016), and Sawicki et al. (2019) for the U band; Gardner et al. (1996), Williams et al. (1996), Benítez et al. (2004), and Laigle et al. (2016) for the B and V bands; Hogg et al. (1997) and Laigle et al. (2016) for the R band; Gardner et al. (1996), Postman et al. (1998), Williams et al. (1996), Benítez et al. (2004), and Laigle et al. (2016) for the I-band; and Laigle et al. (2016) for the z band.

We fit the observed differential counts per unit area and per unit magnitude, N(m), with both a quadratic curve and a sequence of four power laws covering the magnitude range 14–28. Both functional forms give reasonable representations of the galaxy number counts and are shown along with the observations in Figure 13 for the UBVRIz passbands. The IGL, expressed as a V-band surface brightness, is given in Table 3. Our estimate for the surface brightness of faint galaxies was derived using Equation (3) with integration limits running from mag = 30 to mag = 19.1, and we used the best multi-power-law fits to the observed, completeness-corrected number counts from various surveys. The errors in the galaxy surface brightness limit listed in Table 3 include uncertainties in number counts, uncertainties in the best-fit parameters, reasonable variations in the form of those best-fit functions, and an estimate of the cosmic variance (see next paragraph). The uncertainty in the faint-end slope of the log(N)–magnitude relation is the dominant source of error in the derived IGL. For the LORRI images used in this study, incident light from Milky Way stars below the detection limit is subdominant (by a factor of ∼5) compared to that from the IGL.

Zoom In Zoom Out Reset image size

In contrast to the stellar population, the galaxy population tends toward bluer SEDs at fainter magnitudes. We used the COSMOS multiband survey (Laigle et al. 2016) to generate mean galaxy SEDs in bins 1 mag wide from V = 19 to V = 25 using the UBVRIJHK absolute magnitude data provided in the COSMOS catalog. Below V = 25, the catalog is less complete and we just assumed the SEDs for galaxies with V > 25 are identical to those in the 24 < V ≤ 25 bin. We then followed the same prescription discussed in Section 4.1 using Equation (4) to compute the signal in each bin using the appropriate SED and then generate the final IGL signal from the number count–weighted sum of the C_L values for each bin. The IGL intensities derived in this way are listed in Table 4. We assumed that the IGL signal does not vary significantly between LORRI fields. This is a reasonable assumption since the $\sqrt{N}$ variations in the cumulative galaxy number counts between V = 19.1 and V = 30 are <0.2% and the effect of cosmic variance on the galaxy number counts on 03 scales down to such depths is estimated to be ±9.4% integrated over the redshift range 0 < z ≤ 6 based on the approach presented in Trenti & Stiavelli (2008). Cosmic variance is included in the ±0.32 mag arcsec⁻² uncertainty in our surface brightness estimate (Table 3) and in the total uncertainty of  ±2.20 nW m⁻² sr⁻¹ in the derived IGL intensity (Table 4). The IGL intensity is only weakly dependent on the integration limit used provided that limit reaches at least V = 30. If we extend the integration of the galaxy number counts to fainter limits, the predicted IGL intensity will increase, but the rise depends on the faint-end slope of the galaxy counts (see Figure 19). If the faint-end slope values observed at V < 29, which typically lie in the range 0.25–0.35, continue to V > 30, then the derived IGL intensity would increase only by 4%–13% even if we extend the integration down to V = 34. Such an increase is covered by the ∼30% errors on our present IGL estimate. A further discussion of this issue is found in Section 5.2.

LORRI is sensitive to sunlight scattered into the detector even out to solar elongations of 90°, beyond which the LORRI aperture is in the spacecraft's shadow (Cheng et al. 2008; Weaver et al. 2020). This implies that SSL may also be an important contribution to the total-sky level, as is demonstrated in the next section. Figure 14 shows the complete composite point-spread function (PSF) / scattering function for LORRI, which describes the radial distribution of scattered light from the pixel scale to large angles. The inner part of the function was determined by LORRI images of stars in an open cluster, augmented with exposures of Vega and Arcturus to extend the PSF out to the angular limits of the LORRI field of view. Prelaunch calibration extended the distribution to the few-degrees scale (with a small gap that was interpolated over). The large-angle portion of the scattering function was based on the amplitude of scattered sunlight as a function of solar elongation, where the light intensity at any angle was taken as the median flux level over the full LORRI field. We allocated a 10% uncertainty for the PSF amplitude at any location, based on observed variations in the scattered-light amplitude with the azimuth at a constant solar elongation and fine-scale noise in the profile as a function of the radius. In passing we note that the LORRI photometric calibration, which we use throughout this paper, was established by integrating the light of standard stars within a 10'' aperture on the CCD. Integrating the PSF over the full sky yielded a total flux 7.8% larger than the integral flux within this aperture.

Zoom In Zoom Out Reset image size

4.3.1. Estimation of the SSL Contribution

The first step in estimating the SSL contribution was to assess the maximum angular scale we needed to include in our calculations. While Figure 14 indicates scattered light can be seen as much as 90° off-axis, the amplitude of the signal is decreasing rapidly with increasing angular distance from the field center. An initial assessment of each New Horizons field indicated that we would need to account for scattered light from stars as far as 20°–45° from the field center. This is highlighted in Figure 15, which shows the relative contribution of scattered light for stars in the vicinity of the n3c61f field. We found that in our seven fields the stars with angular separations in the range 10°–20° contribute up to 6% of the total SSL signal. Stars in the range 20°–45° contribute no more than 1.3% of the total SSL signal. Beyond 45° even the brightest known stars contribute negligibly (<0.001%) to the detected sky level.

Zoom In Zoom Out Reset image size

The next step was to assess how faint a star and how complete a star catalog we would require. The most robust star catalog currently available is the Gaia DR2 catalog (Gaia Collaboration et al. 2016, 2018), which is complete over the entire sky down to V ∼ 20 mag. Extracting stars to this depth out to 45° for each field would involve a very large number of sources (∼13 million stars for each field in our study). Hence, we performed a test to see if we could use stars to this depth out to a smaller angular scale and supplement that list with a shallower catalog out to the full 45°. We used the ZL 138b field as a test case. We extracted all Gaia stars down to V = 19.5 out to an angular separation of 20°, supplemented with any missing bright stars (2 ≤ V ≤ 8) from the Tycho2 star catalog (Høg et al. 2000) and, for V < 2, from the Yale Bright Star catalog v5.0 (YBSC5; Hoffleit & Warren 1995). We computed the scattered-light signal (details given below) and compared that with the scattered-light signal derived from a sample with Gaia stars down to the same depth but only extending out to 10° and supplemented with the Tycho2 catalog down to V = 10.75 over the angular separation range 10° < θ ≤ 45°. The scattered-light signals from the 20° deep Gaia–Tycho2–YBSC5 sample and the 45° Gaia–Tycho2–YBSC5 catalog agree to within 0.35%. We thus conclude the hybrid catalogs (deep Gaia data out to 10° supplemented with shallower Tycho2 and Yale Bright Star data) would allow for a robust estimate of the scattered-light signal and significantly reduce the amount of star data needed (about 730,000 stars per field instead of >10 million). Our adopted procedure is described in detail below.

To perform the SSL estimation, we identified the stars for each field with V ≤ 19 that lie within an angular distance of θ ≤ 10° from the mean LORRI pointing (see Table 1). We supplemented this list with all stars with V ≤ 10.75 that lie at 10° < θ ≤ 45° from each field center. We used the Gaia DR2 catalog to identify stars with θ ≤ 10° and 8 < V ≤ 19. To ensure we had complete coverage of all known bright stars, we used the Tycho2 catalog to identify stars with θ ≤ 10° and 2 ≤ V ≤ 8 or with 10° < θ ≤ 45° and 2 ≤ V ≤ 10.75. We merged the data from the two catalogs with care to reject common objects. We used the established flux transformations to convert the Gaia and Tycho2 photometry to the Johnson V system. Any Tycho2 and Gaia DR2 entries that lie within 3'' of each other and have an apparent magnitude difference of ΔV ≤ 0.20 mag were considered to be duplicate entries. For any duplicate entries, we adopted the Gaia DR2 values for the position and flux. Lastly, we used YBSC5 to include stars with V < 2 and θ ≤ 45°. This procedure generated a star catalog for each field that was complete to V = 19 out to an angular distance of 10° from the field center and to V = 10.75 for angular distances from 10° to 45°. We also derived an estimate of the (B − V) color for each star using the published color transformations. For Gaia data, we used the DR2 transformation to compute (V − R) and then used the TRILEGAL simulations to derive the transformation from (V − R) to (B − V). The total uncertainty in our (B − V) colors is ∼0.1 mag.

For stars within 1.5 LORRI field widths of the center (261) we accounted for SSL variation over the field as well as for faint extended wings of stars within the fields that extended beyond the 12'' masking radius applied to the cores of the stellar PSFs. For this inner circle around each field center, we constructed an image of the star field from the combined star catalogs and convolved it with the composite PSF (Figure 14) and then measured the integrated SSL contribution over the field. For our adopted photometric detection limit of 19.1 V mag, all the stars masked within the fields are included in the Gaia catalog. The SSL value due to stars within 1.5 LORRI field widths was determined from the mode of the pixel intensity histogram in the above convolved image, with star centers masked out using the same process as was done for the observations. This procedure accurately allowed us to account for any gradients in the SSL across the LORRI field of view. Color corrections, needed due to the broad-wavelength response of the LORRI instrument, were applied for stars with different (B − V) values using the methodology discussed below.

For stars with angular separations larger than 1.5 LORRI field widths, we simply calculated the scattered-light contribution of any star relative to its separation from the center of the LORRI field under the assumption that any SSL gradients across the LORRI field of view from these more angularly distant stars were negligible. ²⁵ To compute the SSL contribution for stars with angular separation larger than 1.5 LORRI field widths, we estimated the total V-band surface brightness reaching the central LORRI 4 × 4 pixel due to SSL as a function of the (B − V) color, μ_V,SSL(B − V), using the following expression:

$$\begin{eqnarray}&&\begin{array}{l}{\mu }_{V,\mathrm{SSL}}(B-V)\\ =\,-2.5{\mathrm{log}}_{10}\left[\frac{1}{{\left(4\mathop{.}\limits^{{\prime\prime} }08\right)}^{2}}\displaystyle \sum _{i=1}^{i=N}\mathrm{PSF}(\theta )\ \mathrm{dex}(-0.4\,V)\right]\\ \quad +\,18.88\end{array}\end{eqnarray} \tag{ 5 }$$

where the sum is over all stars within a given (B − V) range, PSF(θ) is the value of the composite PSF for an angular distance, θ, between the star and the LORRI field center, V is the V-band magnitude of the star, and 18.88 is the appropriate zero-point for a 4 × 4 binned LORRI pixel (Weaver et al. 2020). We computed the value of μ_V,SSL(B − V) for five (B − V) color ranges: [1] (B − V) < 0.00, [2] 0.00 ≤ (B − V) < 0.25, [3] 0.25 ≤ (B − V) < 0.50, [4] 0.50 ≤ (B − V) < 1.00, and [5] (B − V) ≥ 1.00. We used the TRILEGAL simulations to generate SEDs for these same five color ranges for stars with 5 ≤ V ≤ 19. We renormalized these SEDs to match the V-band surface brightnesses computed using Equation (5) and then used Equation (4) with the appropriate renormalized SED to compute the predicted LORRI count rates for each (B − V) range. The final step was to sum up the results for all five color bins to get the total predicted SSL signal for each field. This procedure allowed us to account for the sensitivity of the broad LORRI response function to stars with significantly different SEDs. Our estimated SSL signal for each field is listed in Table 4. The uncertainty in the estimated SSL signal is 10%, due to the uncertainty in the determination of the composite PSF.

4.3.2. Estimation of the Scattered Galaxy Light Contribution

Scattered light from bright (V ≤ 19) galaxies outside of the LORRI field of view can also be computed but is expected to be negligible. To make the scattered galaxy light estimate, we first computed the mean V-band surface brightness for bright galaxies in the range 10 ≤ V ≤ 19, μ_V,BG, using Equation (3). We then computed the surface brightness of scattered light from this galaxy population using the expression

$$\begin{eqnarray}\begin{array}{rcl}{\mu }_{V,\mathrm{SGL}} & = & -2.5{\mathrm{log}}_{10}\ \left[\displaystyle \sum _{\theta =0\mathop{.}\limits^{^\circ }03}^{\theta =45^\circ }\mathrm{PSF}(\theta )\right.\\ & & \times \,\left.[\pi ({\theta }_{i+1}^{2}-{\theta }_{i-1}^{2})]\ \mathrm{dex}(-0.4{\mu }_{V,\mathrm{BG}})\right]+18.88\end{array}\end{eqnarray} \tag{ 6 }$$

where the sum is over angular separations from 003 to 45° in annuli 002 in width and the other expressions are the same as defined in Equation (5). We derived the SED for galaxies in this magnitude range using the UBVRIZYJHK data from the COSMOS survey (Laigle et al. 2016) and renormalized this SED to match the μ_V,SGL surface brightness derived with Equation (6). We then used Equation (4) to derive the predicted signal from scattered light from galaxies outside the LORRI fields. The results are presented in Table 4. As anticipated, the contribution to the total scattered light from galaxies is no more than 1/100 of that from stars.

The final optical background correction that we applied was to account for Milky Way starlight scattered by interstellar dust into our line of sight. Following Zemcov et al. (2017), this "DGL," or DGL component, was estimated for any given LORRI field from I₁₀₀, the strength of the 100 μm thermal emission of the IR cirrus present in the field. The IRIS reprocessing of the IRAS full-sky thermal–IR maps (Miville-Deschênes & Lagache 2005) provided the needed 100 μm measures. The conversion between I₁₀₀ and DGL assumes a linear scaling between the two, which should be valid for modest dust optical depth.

This simple proscription, however, belies a number of uncertainties in how it is actually applied to provide an accurate quantitative estimate of DGL, given I₁₀₀. We examined a variety of estimators, as well as searched for systematic biases in the input I₁₀₀ maps. Our initial approach was to evaluate the approach of Zemcov et al. (2017), who attempted to derive a DGL estimator based on the theoretical scattering properties of plausible models of the dust grains. Allowing for a slowly varying phase function dependent on Galactic latitude, they gave the DGL scattered-light signal in the LORRI passband as

$$\begin{eqnarray}&&{I}_{\mathrm{DGL}}={C}_{100}\,\displaystyle \frac{{I}_{100}}{1\,\mathrm{MJy}\,{\mathrm{sr}}^{-1}}\,\displaystyle \frac{1-0.67\sqrt{\sin | b| }}{0.376}\,\,\mathrm{nW}\,{{\rm{m}}}^{-2}\,{\mathrm{sr}}^{-1},\end{eqnarray} \tag{ 7 }$$

where the term on the right gives the dependence of the conversion on Galactic latitude, b, and is normalized to unity at b = 60°. Zemcov et al. (2017) noted that it is important to subtract 0.78 MJy sr⁻¹ from the IRIS flux values to correct for the contribution of the cosmic IR background (CIB; Fixsen et al. 1996; Puget et al. 1996). The specific value of the leading conversion constant that they derived is C₁₀₀ = 9.8 ± 3.9. ²⁶

The nearly ∼40% relative error in the coefficient, however, translates to a large source of uncertainty in our final COB measurements; thus we considered other approaches to estimating DGL. Intriguingly, Brandt & Draine (2012) provided a direct observational measurement of the total DGL flux as a function of wavelength and Galactic latitude, based on correlating residuals in Sloan Digital Sky Survey background spectra with 100 μm flux. Conveniently, they presented their DGL spectra as a scaling between optical and 100 μm flux over a wide optical bandpass that encompasses the LORRI bandpass. Calculating a LORRI-specific conversion scaling requires only integrating the LORRI response function over their DGL spectra. Using the Brandt & Draine (2012) 50° < ∣b∣ < 90° DGL spectrum yields C₁₀₀ = 5.5 ± 1.3, after rescaling the spectrum by their recommended 2.1 ± 0.4 bias correction factor. DGL estimates produced by this methodology have smaller amplitudes than do those provided by the Zemcov et al. (2017) estimator. That said, unfortunately, the particular high–Galactic latitude spectrum may require an even larger bias correction, given the small number of observations used to generate it (T. D. Brandt 2020, private communication). It is thus likely that this methodology underestimates the appropriate DGL correction. Its nominally smaller formal errors do not reflect this potential systematic error.

In passing, we note that we also attempted to estimate DGL corrections directly from our data by correlating the residual-sky level after all corrections besides the DGL correction had been applied with 100 μm flux. Unfortunately, given the amplitude of the errors in the residual-sky levels and the small range of 100 μm flux over our sample, we could not derive a significant scaling coefficient. As emphasized by Zemcov et al. (2017), this approach may be useful for future work with a larger sample of fields, but we will not pursue this further in the present work.

4.4.1. Measurement of the 100 μm Fluxes for the LORRI Fields

The IRIS 100 μm fluxes for the seven LORRI fields are given in Table 1. The fluxes were calculated as an average over the LORRI field of view for the ZL sequences and over a slightly larger area for the DKBO sequences to accommodate the shifting positions of the images needed to track the DKBOs. The dispersions in the map values for all fields were negligible compared to other sources of error.

In attempting to evaluate the quality of our 100 μm fluxes we compared the IRIS measures to those of Schlafly & Finkbeiner (2011). The two sources agreed well for the fields at high ecliptic latitudes but not for two of the fields at low ecliptic latitudes. The IRIS maps indicated markedly stronger 100 μm fluxes, and the implied IRIS-based DGL corrections yielded residual-sky fluxes with a larger variance over the sample than that measured prior to the correction. Inspection of the IRIS maps showed artifacts associated with imperfect correction for ZL close to the ecliptic plane, however. By plotting all IRIS 100 μm flux values in the maps at Galactic latitude ∣b∣ ≥ 60° as a function of absolute ecliptic latitude (Figure 16), we found a trend indicative of residual zodiacal dust at low ecliptic latitudes. The Schlafly & Finkbeiner (2011) maps show a qualitatively similar trend, but of lower amplitude.

Zoom In Zoom Out Reset image size

We computed a correction for the residual ZL component in the IRIS data by first finding the median fluxes in 1° wide ecliptic latitude bins running in the range 0° ≤ ∣β∣ ≤ 60° for the 5.6 million pixels in the IRIS 100 μm sky map with Galactic latitude ∣b∣ ≥ 60° and Galactic longitude 55° ≤ l ≤ 355° as shown in Figure 16. The IRIS 100 μm flux and median data were provided to us by C. Bot (2020, private communication). A constant CIB level of 0.78 MJy sr⁻¹ (Fixsen et al. 1996; Puget et al. 1996) has been subtracted from all IRIS flux values. We then fit an inverse sigmoid function to the median flux versus ecliptic latitude data of the form

$$\begin{eqnarray}&&{f}_{100}(| \beta | )={f}_{\max }-\displaystyle \frac{({f}_{\max }-{f}_{\min })}{(1.0+{\left({e}^{\kappa (| \beta | -{\beta }_{0})}\right)}^{\gamma })}\end{eqnarray} \tag{ 8 }$$

where f₁₀₀(∣β∣) is the predicted median 100 μm flux in megajanskys per steradian at ecliptic latitude, β, in degrees. The best-fit sigmoid parameters are f_min = 0.71 MJy sr⁻¹, f_max = 2.09 MJy sr⁻¹, κ = 0.20, β₀ = 2040, and γ = −0.73. The best fit is shown in Figure 16 as a blue dashed curve in the upper plot. The correction for the excess zodi signal in the IRIS data was then just ${f}_{100}(| \beta | )-{f}_{\min }$, which was subtracted from our measured fluxes. The motivation for using a sigmoid function to fit this trend was simply that it nicely models a continuous transition between two constant levels. There may well be alternative functions that could be used to model the observed trend. The corrected IRIS fluxes are given in Table 1 and are plotted in Figure 16.

We calculated the DGL corrections based on both the Zemcov et al. (2017) and Brandt & Draine (2012) conversion coefficients (Table 4), given the corrected IRIS 100 μm fluxes. Figure 17 shows the sky levels as a function of the average 100 μm flux after the total star and galaxy light contributions have been subtracted, before and after the two sets of DGL components have been subtracted. The final residual-sky levels after the DGL subtraction represent an estimate of the dCOB in each field. The range in the final dCOB values serves to emphasize the likely effects of uncertainty in the DGL correction.

Zoom In Zoom Out Reset image size

We did not expect to detect any sunlight scattered by IPD within the Kuiper Belt, based on detailed models of the distribution of dust within the outer solar system (Zemcov et al. 2018; Poppe et al. 2019), normalized by in situ measurements by the New Horizons Student Dust Counter (Piquette et al. 2019). In brief, the models integrate scattered light along any line of sight from the New Horizons spacecraft, using an appropriate phase function at any given solar elongation for an assumed dust particle size/density function. For example, the predicted scattered-sunlight flux at 40 au and 90° solar elongation is  ∼0.1 nW m⁻² sr⁻¹,  well over an order of magnitude below our present sensitivity limit (Figure 1). We do not include this small term in our analysis.

Despite these arguments, as noted in Section 2 the motivation for observing the "ZL" fields was a simple test to see if there was any evidence for light scattered by dust within the Kuiper Belt independently of any models. The weighted-average residual skies of the two zero–ecliptic latitude fields (OE 394 was excluded, given its intermediate ecliptic latitude) less those of the three highest–ecliptic latitude fields ²⁷ are 3.5 ± 2.2 nW m⁻² sr⁻¹ for the Zemcov DGL and 4.0 ± 2.1 nW m⁻² sr⁻¹  when assuming the BD2012 DGL.

An IPD measure at this low level of significance is of potential interest; however, it is highly sensitive to our assumed ZL correction to the 100 μm flux maps. Prior to the ZL correction, we obtained an IPD scattered-light level consistent with zero; in short, reducing the DGL attributed to the low–ecliptic latitude fields implies more residual signals to be accounted for. Paradoxically, correcting the low–ecliptic latitude fields for this IPD signal increases the significance of the dCOB level (to be discussed in the next section) by bringing these fields into better concordance with the remaining sample and hence reducing the amount of random variance encompassed by the full sample. We thus do not include this potential IPD correction in our dCOB analysis.

5.1.1. Measuring the Total-sky Signal with LORRI

5.1.2. Starlight

5.1.3. DGL from the Milky Way

The IGL—integrated light from galaxies that fall within the LORRI fields but are too faint to be detected individually—nominally accounts for ∼50% of our total COB signal. The difference between our measurements of the COB and the IGL is significant at the ∼2σ level. At this level of significance we indeed cannot falsify the hypothesis that our COB measure is solely explained by the IGL without recourse to an additional dCOB component. Conversely, it may be possible that our COB measure admits evidence for a visual-light source not explained by the present inventory of galaxies. We explore these competing hypotheses in this and the next subsection.

In attempting to understand the IGL component we again benefit from an immense corpus of work done on galaxy counts in multiple optical bandpasses that allows us to compute constraints, including extrapolation to fainter fluxes (≥28 mag), where direct observational constraints are only now just becoming available. The two key factors in setting the value of the derived IGL are the normalization and the faint-end slope of the relationship between the logarithms of the differential galaxy counts as a function of apparent magnitude. Figure 19 shows how our derived IGL would change as a function of the faint-end slope assumed for the V-band galaxy count–magnitude relation. The trends in the figure were derived by setting the normalization of the galaxy number counts to match that observed at V = 23 (where the incompleteness corrections are negligible) but then allowing the slope to vary from 0.2 to 0.6 at fainter magnitudes. The IGL, integrated down to V = 30 mag, would reach the observed COB level if the faint-end slope were ∼0.48. If the IGL is derived from an integration down to V = 34 mag, the faint-end slope would need to be at least 0.42 to match the observed COB level. The observational constraints on the faint-end slope are denoted by the vertical blue lines in Figure 19. While such steep faint-end slopes are not conclusively ruled out, the observations and theoretical predictions strongly favor shallower slopes, with values of <0.38. At V > 27 in the deepest existing galaxy surveys, the number count–magnitude relations typically have slopes of <0.25.

Zoom In Zoom Out Reset image size

Alternatively, one can ask what normalization shift in the number counts would be required to have the IGL explain the full COB signal. We estimate that, for the observed faint-end slope of ∼0.32, the differential galaxy counts would have to be underestimated by a factor of 2.0 ± 0.5. Since the differential galaxy count data we use have already been corrected for incompleteness, such a normalization shift would imply that existing deep galaxy surveys are systematically missing about half of the actual galaxies with V < 30.

Conselice et al. (2016) performed a comprehensive study of how many galaxies could occupy the visible universe. They estimated that there may be up to 10 times the number of galaxies currently counted in existing deep surveys. This difference is driven by their extrapolation of the stellar mass function down to 10⁶ M_⊙, and thus most of those galaxies would have apparent magnitudes fainter than 30 mag. They would nonetheless contribute to the total COB. Indeed, Conselice et al. (2016) provided an accompanying estimate of the IGL that is 5.8 × 10⁻⁹ erg s⁻¹ cm⁻² Å⁻¹ sr⁻¹ for an integration down to 30 mag, which, at a wavelength of 0.608 μm, corresponds to 35 nW m⁻² sr⁻¹. Their estimate of the IGL is comparable to our total-sky level even before correction for the known foregrounds. Our observational limits on the COB would rule out that specific value for the IGL at very high significance, but as their integrated number counts appear to have about a 35% uncertainty, their IGL estimate may yet be compatible with our observational constraint. The difference between the Conselice et al. (2016) estimate of the IGL and our observations might, however, suggest that extending the Schechter form for the galaxy mass function down to 10⁶ M_⊙ is not fully warranted.

An independent assessment of the IGL is provided by analysis of data from the High Energy Stereoscopic System (HESS Collaboration et al. 2013). HESS detected very high energy (>100 GeV) γ-rays emitted by a sample of blazars over a 0.03 < z < 0.19 redshift range. At these energies, the space density of optical photons can provide a significant cross section for the production of electron–positron pairs, thus attenuating the γ-ray flux. While this is an indirect constraint on the COB flux density, it truly probes the extragalactic optical photon density and is not vulnerable to solar system ZL or even Milky Way optical foreground emission. Furthermore, observing background sources at a range of redshifts verifies that the strength of the attenuation increases with path length. The integration of our galaxy counts over the magnitude range [12, 30] agrees quite well with the HESS constraints, as shown in Figure 20. If the COB is due solely to the collective light from faint galaxies, then the 2σ difference between our derived IGL and our estimate of the COB would imply a factor of ∼2 undercount of galaxies even at apparent magnitudes well within the grasp of current telescopes.

Zoom In Zoom Out Reset image size

Finally, we note that although we needed to estimate how the SEDs of galaxies vary with galaxy magnitude to account for the wide LORRI bandpass, the available observations provide reasonably good constraints on that variation, and in the end, any uncertainty in this step of the analysis is not as significant as uncertainties in the galaxy number counts discussed above.

The present COB measurement appears to be in good agreement with previous measurements but probes an interesting part of parameter space that should allow for new constraints on a true diffuse background of cosmological origin. Figure 20 shows our measurements in the context of previous COB measurements relevant to the LORRI bandpass. In reviewing the literature, we noted that our emphasis on detecting a dCOB signature differs from that of most studies, which prefer to measure the total COB and compare it to the IGL inferred from deep galaxy counts.

For reference, Figure 20 also includes our estimate of IGL as a function of wavelength based on the galaxy counts presented in Section 4.2. The IGL shown here is independent of the LORRI bandpass. This an important caveat when evaluating Figure 20. The actual IGL value we subtracted from our measures is an integral over the LORRI bandpass and thus is slightly different from the IGL presented in the figure, which is done here for comparison to all previous COB observations. The IGL trace shown was derived by integrating the multi-power-law fits to the galaxy number counts (see Figure 13) in the UBVRIz bands over the magnitude interval [12, 30] using Equation (3) and converting the resulting surface brightness to intensity (λI_λ ) at the effective wavelengths of each passband. The red curve in Figure 20 is a spline fit to these intensities. The position of the COB measurements relative to the IGL contribution indicates the contribution of a possible diffuse component.

Our present observations are obviously most closely related to the earlier work done using LORRI by Zemcov et al. (2017). Our respective results are compatible, even though there are substantial differences in the image set, data processing, and analysis methodologies. We based our study on one hundred ninety-five 30 s exposures spread over seven fields, while Zemcov et al. were only able to identify four single 10 s archival images of four different fields that were suitable for their use; the subsequent improvement in material suitable for deep sky measurements is due to the DKBO and ZL programs initiated after their analysis was published. Our program also benefited from improved understanding of the low light–level performance of LORRI and knowledge of its sensitivity to scattered light. This said, Zemcov et al. represented their measures conservatively as a 2σ upper limit on the extragalactic background light (EBL) of 19.3 nW m⁻² sr⁻¹,  which readily accommodates our COB and dCOB measurements. Zemcov et al. also produced an EBL estimate of 4.7 ± 10.3 nW m⁻² sr⁻¹,  which is not significantly different from our result (not shown in the figure).

Bernstein (2007) used Hubble Space Telescope / WFPC2 images to measure a COB level of 55 ± 28 nW m⁻² sr⁻¹ at 0.55 μm and 57 ± 33 nW m⁻² sr⁻¹ at 0.81 μm. While these flux levels are markedly higher than our COB result, with their large error bars this difference is not highly significant. We note that our total-sky level, 33.2 ± 0.5 nW m⁻² sr⁻¹,  prior to any corrections is itself only ∼60% of the Bernstein (2007) COB signal.

In contrast, the COB measurement of 7.5 ± 5.8 nW m⁻² sr⁻¹ at 0.64 μm provided by Pioneers 10 and 11 photometry (Matsuoka et al. 2011) is ∼50% of our value but again has low significance. This work is of particular interest, however, as the photometry was obtained when the spacecraft were 3.2 to 5.2 au from the Sun, allowing the authors to claim that the effects of ZL were negligible. More recently, however, Matsumoto et al. (2018) re-reduced the Pioneer 10 data and claimed that instrumental artifacts prevent accurate measurement of the absolute sky level.

Matsuura et al. (2017) obtained low-resolution 0.8 μm ≤ λ ≤ 1.8 μm NIR spectra of the sky using the Cosmic Infrared Background Experiment (CIBER) rocket instrument. They argued that their COB measurements were significantly above the expected IGL contribution and had a SED that could not be matched by improperly subtracted ZL. While the center of our bandpass is well to the blue of their spectral range, it does overlap with the CIBER observations, and a dCOB signal at <0.9 μm would contribute to our sky. The dCOB level that we see is compatible with their measurements.

The HESS Collaboration et al. (2013) COB constraints are represented as a gold-colored band in Figure 20. The HESS Collaboration emphasized the close concordance of their results with the known IGL component, which is evident in our figure. At the same time, it is also clear that the HESS EBL measurements would allow a dCOB component, consistent with our result, particularly in the NIR.

The previous COB measurements discussed so far are estimates of an overall constant background level. A different technique is to look for angular fluctuations in the background to assess light associated with sources fainter than the photometric detection limit. An advantage to this approach is that it rejects ZL, which is assumed to be smooth on fine scales. Zemcov et al. (2014) analyzed the NIR CIBER images and argued for a dCOB component (in our parlance) of ${7.0}_{-3.5}^{+4.0}\,\mathrm{nW}\,{{\rm{m}}}^{-2}\,{\mathrm{sr}}^{-1}$ at 1.1 μm, which is consistent with the present results. They argued that this might be evidence for intra-halo light, which would be a relatively local component of stars tidally ripped out of merging galaxies (Cooray et al. 2012). Matsumoto & Tsumura (2019), analyzing optical-band fluctuations in the Hubble Extremely Deep Survey (Illingworth et al. 2013), also argued for a significant dCOB, which would be generated by faint compact objects at z ≲ 0.1. A key argument for sources at low redshift is to accommodate the very high energy γ-ray opacity limits on COB components at z ≫ 0.1.

We finish this discussion on the possibility of a dCOB by returning to our discussion of Conselice et al. (2016) in the previous subsection. The central thesis of that work is that existing galaxy counts fall an order of magnitude short of accounting for the total visible-light output of stars integrated over the history of the universe. While Conselice et al. (2016) in fact predicted a COB substantially in excess of what our observations and the HESS results can accommodate, their conclusion that the present view afforded by galaxy counts may not be all there is to the story remains an observational challenge.