Heliocentric distance dependence of zodiacal light observed by Hayabusa2#

Abstract

Zodiacal light (ZL) is sunlight scattered by interplanetary dust particles (IDPs) at optical wavelengths. The spatial distribution of IDPs in the Solar System may hold an important key to understanding the evolution of the Solar System and material transportation within it. The number density of IDPs can be expressed as $n(r) \sim r^{-\alpha }$, and the exponent $\alpha \sim 1.3$ was obtained by previous observations from interplanetary space by Helios 1/2 and Pioneer 10/11 in the 1970s and 1980s. However, no direct measurements of $\alpha$ based on ZL observations from interplanetary space outside Earth’s orbit have been performed since then. Here, we introduce initial results for the radial profile of the ZL at optical wavelengths observed over the range 0.76$-$1.06 au by ONC-T aboard the Hayabusa2# mission in 2021-2022. The ZL brightness we obtained is well reproduced by a model brightness, although there is a small excess of the observed ZL brightness over the model brightness at around 0.9 au. The radial power-law index we obtained is $\alpha = 1.30 \pm 0.08$, which is consistent with previous results based on ZL observations. The dominant source of uncertainty arises from the uncertainty in estimating the diffuse Galactic light (DGL).

Graphical Abstract

Introduction

The zodiacal light (ZL) is sunlight scattered by interplanetary dust particles (IDPs) at optical wavelengths, and it is a major constituent of the diffuse celestial brightness. A continuous supply of IDPs is necessary to sustain the diffuse brightness because IDP is removed from the Solar System due to the Poynting–Robertson (PR) effect and by radiation pressure from the Sun (Wyatt and Whipple 1950; Burns et al. 1979). Possible sources for this supply are asteroid collisions (Dermott et al. 1984; Schramm et al. 1989; Tsumura et al. 2010) or cometary ejections (Liou et al. 1995; Nesvorný et al. 2010; Yang and Ishiguro 2015), but the relative ratios of the contributions from these sources are still unknown. Dust of interstellar origin also contributes $\sim$10% to the total amount of IDP (Rowan-Robinson and May 2013). Thus, observational constraints that can tell the differences among these sources are important for a better understanding the origin and characteristics of IDPs and of the way planetary and exoplanetary systems evolve with time (Leinert et al. 1998; Lasue et al. 2020).

Historically, extensive ZL observations were conducted from ground-based telescopes at high-altitude sites in the 1960 s and 1970 s (Dumont and Sanchez 1975; Levasseur-Regourd and Dumont 1980), but the accuracy of these ZL observations is limited due to atmospheric emission. In contrast, space-based platforms eliminate atmospheric contamination and provide precise ZL measurements (Murdock and Price 1985; Matsuura et al. 1995; Matsumoto et al. 1996; Tsumura et al. 2010, 2013a; Buffington et al. 2016; Korngut et al. 2022; Takimoto et al. 2022, 2023). The ZL is the only sky-brightness component that is not fixed on the celestial sphere. In general, the ZL is smoothly distributed, and its small-scale spatial structures are only at the level of a few percent owing to the smooth spatial distribution of IDP as a smooth cloud (Pyo et al. 2012). The plane of symmetry of the smooth cloud is slightly inclined to the ecliptic plane because of the Jovian orbit. Seasonal variations in ZL occur for an Earth-based observer due to the orbital motion of the Earth, which changes the heliocentric distance and the position of the observer with respect to the symmetry plane. A detailed IDP distribution model has been established based on the seasonal variation of the ZL (Kelsall et al. 1998; Wright 1998).

The number density (n) of IDP is presumed to be of a form that is separable into radial and vertical terms:

$$\begin{aligned} n(r, \beta ) = n_{0}\left( \frac{r}{r_0}\right) ^{-\alpha } f(\beta ), \end{aligned}$$

(1)

where $n_0$ is the reference number density of IDPs in the symmetry plane at the heliocentric distance $r_0$, and $f(\beta )$ denotes the vertical distribution as a function of an elevation angle $\beta$ from the symmetry plane (Giese et al. 1986). The assumption that the vertical distribution of the IDPs depends only on $\beta$ is suggested by the fact that the PR effect does not affect the orbital inclinations of particles as they spiral into the Sun. The radial power-law is induced by the radial distribution expected for particles under the influence of the PR effect, which results in $\alpha =1$ for dust bound in a circular orbit (Burns et al. 1979). When dust-grain sizes are reduced by sublimation near the Sun, such smaller dust particles are expelled from the Solar System as $\beta$-meteoroids by radiation pressure (Zook and Berg 1975; Wehry and Mann 1999; Krüger and Grün 2014). The radial profile of $\beta$-meteoroids is expected to follow a power law with $\alpha =2$ (Szalay et al. 2020). The relative ratio of these two components remains an open issue and may hold an important key for understanding the evolution of the IDP distribution (Leinert and Grün 1990; Mann et al. 2004).

The ZL brightness $I_{ZL}$ can be modeled as the integral of scattered sunlight along the line of sight:

$$\begin{aligned} I_{ZL} = \int F_{\odot }(r) n(r) A \Phi (\theta ) dl, \end{aligned}$$

(2)

where $F_{\odot }(r) \sim r^{-2}$ is the solar flux at the distance r from the Sun, A is the albedo of the IDP, $\Phi (\theta )$ is the phase function at the scattering angle $\theta$, and dl is an increment along the line of sight. If the scattering properties (size and albedo) of IDPs do not change significantly with heliocentric distance, the heliocentric dependence of ZL toward the antisolar direction on the symmetry plane can be written as $I_{ZL} \sim r^{-(\alpha +1)}$. If the line of sight is not oriented in the antisolar direction, the heliocentric dependence of $I_{ZL}$ becomes much more complex, since it depends on the phase function $\Phi (\theta )$ for which $\theta$ will vary. In addition, a heliocentric dependence of the local albedo of the IDPs has also been reported (Levasseur-Regourd et al. 1991), which makes the heliocentric dependence of $I_{ZL}$ even more complex.

Direct observations of the radial power-law index $\alpha$ based on ZL observations were performed from spacecraft outside Earth’s orbit in the 1970 s and 1980 s. Pioneer 10/11 observations of the ZL at 1$-$3.3 au gave $\alpha =$1$-$1.5. More specifically, a single power-law model with $\alpha \sim 1$ and a cutoff near 3.3 au gives the best fit to the observational data, although a two-component model with $\alpha \sim 1.5$ and increased IDP in the asteroid belt fits the data equally well (Hanner et al. 1976). Helios 1/2 observations of the ZL at 0.3-1 au gave $\alpha =1.3 \pm 0.05$, although $\alpha = 1.35$ gives a better fit for small solar elongations ($<50 ^{\circ }$), and $\alpha = 1.25$ is more appropriate for large solar elongations ($>100 ^{\circ }$) (Leinert et al. 1981, 1982). ZL observations from spacecraft outside Earth’s orbit have not been performed following these missions. The Japanese Venus orbiter Akatsuki tried but could not detect the ZL due to insufficient cooling of the sensor (Satoh et al. 2016).

Some IDP distribution models were developed based on observations of the all-sky ZL brightness and its seasonal variation from geocentric orbit. In particular, observations from the Cosmic Background Explorer (COBE) yielded $\alpha = 1.34 \pm 0.022$ (Kelsall et al. 1998) and $\alpha = 1.22$ (Wright 1998), and observations by AKARI gave $\alpha = 1.59 \pm 0.02$ (Kondo et al. 2016). These observations of the ZL were performed at 1 au, so the accuracy in determining $\alpha$ was worse than that obtained by direct observations from interplanetary space.

The value of $\alpha$ has also been determined based on the observations of the inner ZL or F-corona. Observations of the inner ZL by Clementine from lunar orbit while the Sun was in eclipse behind the Moon yielded $\alpha = 1.45 \pm 0.05$ (Hahn et al. 2002). Values of $\alpha$ from 1.31 to 1.35 were obtained from F-corona observations at elongations ranging from 0.07 to 0.45 au from the Sun (Stenborg et al. 2018) by the Heliospheric Imager-1 (Eyles et al. 2009) onboard the Solar TErrestrial RElations Observatory-A (STEREO-A) orbiting the Sun at approximately 1 au. In addition, $\alpha = 1.31$ was obtained by F-corona observations between 0.1 and 0.4 au (Howard et al. 2019) by the Widefield Imager for Solar Probe inner telescope (WISPER-1) (Vourlidas et al. 2016) onboard the Parker Solar Probe (PSP) when it passed perihelion at 0.16$-$0.25 au. These results are limited to dust distributions close to the Sun.

A technique for studying the distribution and properties of IDPs independent of the ZL observation is in situ dust counting using dedicated dust detectors. The size distribution of IDPs was studied by the in situ dust-counting method and it was suggested that large (10-100 $\mu$m) dust is dominant around 1 au (Grün et al. 1985; Divine 1993). The ZL brightness is indicative of the IDP distribution in the inner Solar System, where the IDP density is substantial, and the IDP distribution derived from these ZL observations is confined to the inner Solar System ($<5$ au). Conversely, dust distribution in the outer Solar System has been investigated by the in situ dust-counting method (Poppe et al. 2019; Bernardoni et al. 2022).

This paper introduces the IDP distribution based on the ZL observations from the Hayabusa2# mission at 0.76$-$1.06 au performed in 2021-2022. These are the first successful observations of the ZL from outside Earth’s orbit in the last 40 years.

Data acquisition and reduction

Hayabusa2# overview

Hayabusa2 is the second Japanese asteroid-sample-return mission. The Hayabusa2 spacecraft was launched in December 2014 and successfully arrived at asteroid (162173) Ryugu in June 2018. After extensive scientific observations for $\sim$1.5 years, it departed from Ryugu in November 2019 and successfully brought the capsule containing Ryugu samples back to Earth in December 2020 (Tsuda et al. 2022; Tachibana et al. 2022). With the successful main mission of the sample return completed, an extended mission named Hayabusa2# (SHARP; Small Hazardous Asteroid Reconnaissance Probe) was initiated to explore new asteroids; it will perform a fly-by of (98943) 2001 CC21 in July 2026 and a rendezvous with 1998 KY26 in July 2031 (Mimasu et al. 2022). Some scientific observations including ZL observations will be performed during this long cruising phase (Hirabayashi et al. 2021).

The Optical Navigation Camera (ONC) onboard Hayabusa2 consists of one telescopic camera (ONC-T) and two wide-angle view cameras (ONC-W1/W2) (Kameda et al. 2017; Suzuki et al. 2018; Tatsumi et al. 2019; Kouyama et al. 2021; Yamada et al. 2023), and it was used for both global and local high-resolution optical observations of Ryugu (Sugita et al. 2019). The ONC was carefully calibrated both before and after launch, and it remains in good condition after contact with the surface of Ryugu during the two touchdowns for sampling. In this study, we used ONC-T for the ZL observations. The field of view of ONC-T is $6.27 \times 6.27$ deg$^2$, which is covered with a $1024 \times 1024$ pixel region of a CCD detector (Kameda et al. 2017). The longest exposure time of ONC-T was 178 s, which we used for the ZL observations in this study. ONC-T has a wheel system that rotates seven color-bandpass filters and one wide clear filter (a panchromatic glass window). We used the wide-band filter ($\lambda$ = 612 nm and $\Delta \lambda$ = 448 nm, see Fig. 1) for the ZL observations, with a v-band filter ($\lambda$ = 550 nm and $\Delta \lambda$ = 28 nm) for stray light subtraction (see "Stray-light subtraction" section).

Fig. 1

A comparison of the relative transmissivity. The solid line shows the bandpass of ONC-T/Hayabusa2# wide-band (Tatsumi et al. 2019), the dashed line shows the bandpass of Gaia G-band (Riello et al. 2021), and the dotted line shows the bandpass of LORRI/New Horizons (Cheng et al. 2008)

Observation fields

The Hayabusa2 spacecraft followed an elliptical orbit over the range 0.76$-$1.06 au before the Earth swing-by on December 2027, as shown in Fig. 2 and 3 (Mimasu et al. 2022). The spacecraft maintained an attitude in which the solar-array paddle (+Z direction) was pointed toward the Sun during this period, and we performed the ZL observations during periods when the ion engines were not in operation. Since ONC-T points toward the -Z direction, the ZL is observed toward the antisolar direction. This is an advantage of our ZL observations over past observations because previous ZL observations in interplanetary space were made at various solar elongation angles, making it difficult to distinguish whether the ZL changes were due to changes in the heliocentric distance or in the solar elongation. In our ZL observations, the change in ZL brightness due to the solar elongation was minimized by observing the ZL at a nearly constant solar elongation (see Table 1).

Fig. 2

Hayabusa2# orbit in the J2000EC inertial frame before the Earth swing-by on December 2027. The regions suitable for ZL observations are shown in red, and unsuitable regions owing to pointing toward the Galactic plane (Galactic latitude $< 20^{\circ }$) and the Galactic center (Galactic longitude $< 20^{\circ }$) are shown in blue and green, respectively

Fig. 3

a Heliocentric distance of the Hayabusa2 spacecraft and b its ecliptic longitude, c Galactic longitude, and d Galactic latitude of the antisolar direction during the period from April 2021 to March 2023. The regions suitable for ZL observations are shown in red, and unsuitable periods owing to pointing toward the Galactic plane (Galactic latitude $< 20^{\circ }$) and the Galactic center (Galactic longitude $< 20^{\circ }$) are shown in blue and green, respectively

Stray light is produced when sunlight hits the radiator that cools the ONC-T detector from a certain range of directions (Suzuki et al. 2018; Tatsumi et al. 2019). Thus, the ZL observations need to be conducted in a “stray-light-avoidance attitude”, in which the -X side and +Y side of the spacecraft are illuminated by the Sun. For this reason, the actual directions of our ZL observations are shifted from the antisolar direction by $\sim$10 degrees. Table 1 summarizes the observed fields and the position of the spacecraft when the observations were conducted.

Table 1 Spacecraft positions and observed fields

Periods when the observable direction (antisolar direction) is pointed toward either the Galactic plane or the Galactic center are not suitable for ZL observations because the Galactic brightness is too strong (see "Integrated starlight" and "Diffuse Galactic light" sections). Thus, we define as unsuitable ZL observation periods those in which the antisolar direction is pointing toward (1) Galactic latitude $< 20^{\circ }$ (the blue zones in Figs. 2 and 3) or (2) Galactic longitude $< 20^{\circ }$ (the green zones in Figs. 2 and 3). The ZL observations were made approximately once a month during the time suitable for ZL observations (the red zones in Figs. 2 and 3).

Acquired images

The ONC-T detector has a $1024 \times 1024$ pixel imaging region, with a $16 \times 1024$ pixel masked regions termed “optical black” on each side as a dark reference (Kameda et al. 2017). The two optical-black images are combined and treated as one $32 \times 1024$ pixel optical-black image. Raw images acquired by ONC-T are processed in a sequence of steps to calibrate the image data. In this work, we used L2a-level images, which are raw FITS images with header information containing the spacecraft system housekeeping data and ONC status data (the temperatures of the detector, lens system, and electronics as well as the voltages of the electronics, etc.) (Tatsumi et al. 2019). The signal from each pixel is provided in 16-bit digital numbers (DN).

One ZL observation dataset includes one bias image B(x, y), three wide-band images $W_i(x,y)$ ($i=$1-3), two v-band images $V_j(x,y)$ ($j=$1-2), and their respective optical-black images (Bb(x, y), $Wb_i(x,y)$, and $Vb_j(x,y)$), as shown in Fig. 4. For the ZL observations on 2021-11-29, 2021-12-06, 2021-12-28, 2022-01-24, and 2022-02-14, as part of the calibration operations we acquired two data sets to monitor the stability of the ONC-T sensitivity after it was turned on.

Fig. 4

Dataset acquired for one ZL observation. One dataset includes one bias image B(x, y), three wide-band images $W_i(x,y)$ ($i=$1-3), two v-band images $V_j(x,y)$ ($j=$1-2), and their respective optical-black images (Bb(x, y), $Wb_i(x,y)$, and $Vb_j(x,y)$)

Dark-current subtraction

Dark-current subtraction is essential for the measurement of diffuse radiation such as the ZL. We estimated the dark current for each image in our dataset from the corresponding dark image, which we obtained from the optical-black image by subtracting a bias image ($Wb_i(x,y)-Bb(x,y)$ and $Vb_j(x,y)-Bb(x,y)$). We then created a histogram of the dark image, and we take its peak position to be the value of the dark current for that image. We fitted the histogram of the dark image with a Gaussian function. The peak position of the histogram corresponds to the mode of the dark image. Using this Gaussian-fitting procedure, we eliminated bad pixels such as those due to leakage of light from the imaging region or due to hot pixels caused by cosmic-ray hits. We expressed the resulting dark-current values for the wide-band and v-band images as $I_{\text {dark}}^{W_i}$ and $I_{\text {dark}}^{V_j}$, respectively. As an example, we found the dark current of the first bias-subtracted optical-black image taken on 2022-08-29 to be $I_{\text {dark}}^{W_1} = 5.92$ DN from the peak position of the histogram, as shown in Fig. 5.

Fig. 5

Dark-current estimation. A histogram of the first bias-subtracted optical-black image, $Wb_1(x,y)-Bb(x,y)$, obtained on 2022-08-29 (black) and its best-fit Gaussian function (red). We obtained the dark current as ${I_{\text {dark}}^{W_1}}=5.92$ DN from the peak position of this histogram

Next, the obtained dark current and the bias image are subtracted from the imaging region, yielding three wide-band subtracted images and two v-band subtracted images in one data set. From the wide-band images, we generated a single reduced image W(x, y) from three subtracted images by taking the median of each pixel:

$$W(x,y) = {\text{median }}\left[ {W_{i} (x,y) - (B(x,y) + I_{{{\text{dark}}}}^{{W_{i} }} )} \right].$$

(3)

This median procedure removes many hot pixels caused by cosmic-ray hits. Since there are only two v-band images, we cannot use the median procedure for them. Instead, we generated a single reduced image V(x, y) by taking the minimum of each pixel to reduce hot pixels caused by cosmic-ray hits:

$$V(x,y) = {\text{min }}\left[ {V_{j} (x,y) - (B(x,y) + I_{{{\text{dark}}}}^{{V_{j} }} )} \right].$$

(4)

Stray-light subtraction

Since stray light occurs in ONC-T images when the spacecraft is at certain attitudes (Suzuki et al. 2018; Tatsumi et al. 2019), we performed all ZL observations in stray-light-avoidance attitudes (see "Observation fields" section). However, weak stray light remains even in this case. Figure 6a and b shows the reduced wide-band image W(x, y) and v-band image V(x, y) obtained on 2021-08-23, respectively, and the stray-light patterns can be seen clearly in these images. It is known that the intensity and pattern of the stray light do not depend on the filter selection (Suzuki et al. 2018). Thus, we subtracted the v-band image V(x, y) as a stray-light reference frame from the wide-band image W(x, y) to remove the remaining stray light, as shown in Fig. 6c and d. Since the ZL signal in the v-band image is estimated to be less than 1 DN, there is little impact on the scientific analysis of the ZL due to this stray-light-removal procedure.

Fig. 6

Data-reduction procedure. a The wide-band image generated by taking the median of the three dark-subtracted wide-band images. b The v-band image generated by taking the minimum of two dark-subtracted v-band images as a reference frame for the stray light. c The stray-light-subtracted image $W(x,y)-V(x,y)$. d Same as c. e A normalized flat-field image (Suzuki et al. 2018). f The flat-field-corrected image SKY(x, y), which is the final reduced image of the sky used for scientific analysis

We found an additional stray-light pattern in the data obtained on 2022-01-24 and 2022-02-14, as shown in Fig. 7. Because this stray light appears only in the wide-band image, it cannot be removed by the v-band subtraction procedure. The source of this additional stray light is thought to be light scattered at the inner wall of the entrance hole of the ONC-T hood, as is indicated by the shape of the stray light (circular pattern). For this reason, we excluded data from these 2 days from subsequent analyses.

Fig. 7

The additional stray light pattern found in the data obtained on 2022-01-24 and 2022-02-14

Flat-field correction

The stray-light-subtracted image $W(x,y)-V(x,y)$ clearly shows a limb-darkening pattern (Fig. 6c and d), which is the same pattern as in the flat-field image (Fig. 6e) (Tatsumi et al. 2019; Kameda et al. 2021b). This fact means that the detector is uniformly illuminated from the front of the optics, showing that the sky brightness has certainly been detected by the ONC-T. We corrected this limb-darkening pattern by dividing the image by the normalized flat-field image FLAT(x, y), as shown in Fig. 6f. Since we have not created a wide-band flat-field image, we used the v-band flat-field image instead. The wavelength dependence of the flat-field image is negligible because the detector is identical for both bands. The image obtained after the flat-field correction SKY(x, y) is the final reduced image of the sky used for scientific analysis:

$$\begin{aligned} { SKY(x,y) = \frac{W(x,y) - V(x,y)}{FLAT(x,y)} }. \end{aligned}$$

(5)

Sensitivity calibration using stars

Degradation of the ONC-T sensitivity was reported after the two touchdown operations on the asteroid Ryugu (Kouyama et al. 2021; Yamada et al. 2023). Consequently, we monitored and calibrated the sensitivity of the ONC-T in our data using the field stars in our images. For this sensitivity calibration, we used the W(x, y)/FLAT(x, y) image and not the SKY(x, y) image. This is because the flux for the bright stars used in this calibration process in case of the SKY(x, y) image is unsuitable for the sensitivity calibration as the v-band signal of the bright stars has been subtracted in the SKY(x, y) image. Although the W(x, y)/FLAT(x, y) image includes the stray light described in "Stray-light subtraction" section, it can be removed by the aperture photometry procedure described below.

The sensitivity-calibration procedure for the wide-band is as follows. First, we solved the astrometry of the images from the distribution of the stars using the astrometry-calculation code Astrometry.net (Lang et al. 2010). Next, we matched the bright stars in the image with those in the Gaia Data Release 3 (DR3) catalog (Gaia Collaboration 2016, 2022). This catalog is suitable for our dataset because the wavelength coverage of Gaia’s G-band, which is an unfiltered, white-light photometric band, is similar to that of our wide-band filter on ONC-T (Fig. 1). The Gaia DR3 catalog contains around $1.806 \times 10^9$ sources, with a limiting magnitude of about $G \sim 21$ mag, with uncertainties of $\sim$0.3 mmag for $G<13$ mag, 1 mmag at $G=17$ mag, and 6 mmag at $G=20$ mag. We selected stars that meet the following criteria:

• (a) The selected stars are in regions with stray-light intensities less than 20 DN in the V(x, y) image,

• (b) The selected stars are in the region with normalized flat-field values greater than 0.8, and

• (c) The selected stars have fluxes between the 6th and 9th AB magnitude in the G-band in the Gaia DR3 catalog.

Criteria (a) and (b) reduce the uncertainty caused by the reduction processes of stray-light subtraction and flat-field correction. The fraction of the area satisfying both criteria (a) and (b) is approximately 44% of the total detector area in the central region of the detector. Criterion (c) reduces the uncertainty in the photometry by selecting stars that have sufficient signal but are not saturated. The ONC-T detector is known to be linear up to $\sim$3000 DN, with $<1$% deviation (Tatsumi et al. 2019), and the signal value of even the brightest pixel in an image of a 6th-magnitude star is approximately within this range.

For aperture photometry, a boxed region of $41\times 41$ pixels centered on the selected star is cut out, and this boxed region is divided into a region centered on the star and the surrounding background region, as shown in Fig. 8. The radius of the circle used to cut out the star at the center is adjusted according to the G-band brightness of this star in the catalog. Next, all pixels greater than 20 DN are masked, as they are considered to be other astronomical objects (stars and galaxies) or hot pixels caused by cosmic-ray hits. We applied additional masks using a $\sigma$-clipping procedure to remove the remaining bright pixels. Then we examined all the masked images by eye and masked any additional remaining bad pixels. Subsequently, we calculated the background brightness and its noise by computing the average and standard deviation of the masked background region, and we subtracted the background brightness from the star region. We then calculated the flux from the central star by calculating the sum of the masked and background-subtracted star region and estimated its uncertainty based on the background noise. Figure 9 shows the relation between the G-band fluxes of the selected stars from the catalog and their detected signal in the 2022-12-12 data. A good linear relation between them exists in all the observed data, and we obtained the sensitivity in units of (DN/sec)/(W/m$^2$/$\mu$m/sr) by taking their ratio.

Fig. 8

Aperture photometry of a star, as used for sensitivity calibration. A boxed region of $41\times 41$ pixels centered on the selected star from the catalog is cut out from the W(x, y)/FLAT(x, y) image (left), and it is then divided into a region containing the central star (center), which is used for the aperture photometry of the star, and the surrounding background region (right), which is used to estimate the background brightness and its noise

Fig. 9

The sensitivity calibration. The linear relationship between the G-band flux of the selected stars and their detected signals in the W(x, y)/FLAT(x, y) image of the 2022–12-12 wide-band data as an example of the sensitivity calibration. The red line shows the best fit to the data, and the dashed lines show $\pm 1\sigma$ deviations

Figure 10 (left) shows the sensitivity calculated from our observed data as a function of time. The sensitivity obtained after the return to the Earth (Yamada et al. 2023) is also shown. This figure shows that the degradation of the sensitivity has stopped and that the sensitivity has remained almost constant since the return to the Earth. The average and standard deviation of the sensitivity after the $\sigma$-clipping procedure is $16095 \pm 490$ (DN/sec)/(W/m$^2$/$\mu$m/sr) (red solid and dashed lines, respectively, in Fig. 10 left), which we applied to all the data to obtain the sky brightness of the wide-band images. We treated the standard deviation of the sensitivity as a systematic uncertainty (see "Field-variance correction" section).

Fig. 10

Sensitivity variation of the wide-band (left) and v-band (right) images. The red solid lines show the average of these sensitivities and the red dashed lines show their $1\sigma$ standard deviations. The black dashed lines show the sensitivity after the return to Earth (Yamada et al. 2023)

We applied the same sensitivity-calibration procedure to the v-band data (V(x, y)/FLAT(x, y) images) to check the consistency of the result because the wide-band data and v-band data share the same detector. Gaia does not have a V-band filter, but it does have a blue band (BP) and a red band (RP), and the V-band magnitude can be estimated from the BP and RP magnitudes (Riello et al. 2021). In the selection of stars for the v-band calibration, conditions (a) and (b) are the same as for the wide-band calibration, and condition (c) selects stars brighter than 8th mag in the converted V-band. Figure 10 (right) shows the sensitivity profile of the v-band, and we confirmed that the sensitivity remained constant in our dataset. The sensitivity and its 1$\sigma$ uncertainty we obtained for the v-band are $906 \pm 37$ (DN/sec)/(W/m$^2$/$\mu$m/sr).

Point-spread function

We obtained a template for the point-spread function (PSF) of the ONC-T wide-band image by adding the images of several bright stars. First, as suitable images for making the PSF template, we selected 56 objects with G-band magnitudes between 6th mag and 8th mag that have clean stellar images, with few bad pixels or without other objects around them. Then, we aligned these images with 0.1 pixel resolution and added them together to obtain the PSF template. Figure 11 shows the resulting PSF template and its radial profile with 0.1 pixel resolution. The full-width half-maximum (FWHM) of the PSF is 2.00 pixels, which is consistent with previous measurements (Kouyama et al. 2021).

Sky brightness

We obtained the sky brightness as the signal value of dark pixels with no stars nor hot pixels by cosmic-ray hits in the SKY(x, y) image, using the histogram method employed to determine the dark current described in "Dark-current subtraction" section. We created a histogram of the SKY(x, y) image, including stars and hot pixels, and obtained the best-fit Gaussian curve. Figure 12 shows the histogram of the SKY(x, y) image and its best-fit Gaussian curve for the 2022-10-17 data. Since the SKY(x, y) image is dominated by dark pixels with no stars nor hot pixels, the peak of the pixel histogram of the SKY(x, y) image represents the sky brightness. We separate the higher signal tail of bright stars and the dark sky signals in the pixel histogram by the Gaussian fitting as shown in Fig. 12, and we treated the 1$\sigma$ error in the peak position as the statistical uncertainty in the sky brightness. We converted the resulting sky brightness from DN units to $\lambda I_{\lambda }$ in nW/m$^2$/sr units by applying the calibration factor obtained in "Sensitivity calibration using stars" section. Using this procedure, the systematic uncertainty in the calibration factor is transferred to a systematic uncertainty in the sky brightness. Table 2 summarizes the obtained sky brightness and its statistical and systematic uncertainties.

At this point, the detection limit has not been determined (this is obtained in the next subsection using the sky brightness and its standard deviation). As it is not known how faint stars should be masked based on the star catalog, we determined the sky brightness by creating a histogram of the entire SKY(x, y) image without masking the stars. This method worked well because the number of pixels observing dark sky is much larger than the number observing stars and other astronomical objects. Note that the stars detected in the images are masked when we determine the ZL as described in "Zodiacal light" section.

Table 2 Obtained sky brightness and backgrounds with their uncertainties

Limiting magnitude

It is important to know the limiting magnitude in our observed images because we need to mask the detected stars to derive the diffuse brightness of the sky. As shown in Fig. 12, the histogram of the SKY(x, y) image has a side lobe caused by the detected stars in the images, which exceeds the best-fit Gaussian curve. The width ($\sigma$) of the Gaussian corresponds to the standard deviation of the fluctuation in the sky brightness, and we define the limiting magnitude as the brightness of stars for which three pixels in the center of the PSF (Fig. 11) exceed $+2\sigma$ of the sky deviation. We set the uncertainty in determining the limiting magnitude to be 1 DN, which is equivalent to approximately 0.2 mag uncertainty in the limiting magnitude. The limiting magnitude of each image and its uncertainty are summarized in Table 2.

Stars brighter than the limiting magnitude were extracted from the Gaia DR3 catalog (see "Sensitivity calibration using stars" section), convolved with the PSF (see "Point-spread function" section), and distributed in the image to create a Gaia bright-star image, as shown in Fig. 13b. The distribution of stars in this Gaia bright-star image reproduces well the distribution of the stars detected in the image observed by ONC-T (Fig. 13a), indicating the validity of the detection limits determined by the method described above. We used the Gaia bright-star image as a stellar mask to conceal stars when obtaining the ZL (see "Zodiacal light" section).

Fig. 11

Point-spread function. The PSF and its radial profile with 0.1 pixel resolution are shown on a linear scale (left) and a logarithmic scale (right)

Fig. 12

Sky brightness. A histogram of the pixel signals of the 2022–10-17 data is shown on the vertical axis with a linear scale (left) and a logarithmic scale (right). The red curves show the best-fit Gaussian functions

Fig. 13

SKY(x,y) image obtained by ONC-T and ISL image. a A zoom-in image taken by ONC-T. b Gaia bright-star image of the same area as panel (a). We obtained this image using stars brighter than the limiting magnitude of the Gaia DR3 catalog, and we used this as a mask to conceal the stars detected in the SKY(x,y) image. c Gaia ISL image of the same area as panel (a). We constructed this image using stars fainter than the limiting magnitude of the Gaia DR3 catalog. The intensity scale of this image (c) is several times larger because the stars are too faint to be seen at the same scale. d Masked image of SKY(x,y) shown in a obtained using the mask image (b)

Background subtraction

Integrated starlight

Stars fainter than the limiting magnitude are not detected as point sources in the observed images, but the sum of the light from those undetected stars, called integrated starlight (ISL), contributes to the sky brightness. Therefore, the ISL must be estimated and subtracted from the sky brightness to obtain the ZL.

The ISL image ISL(x, y) for each field is produced by stars fainter than the limiting magnitude extracted from the Gaia DR3 catalog, as shown in Fig. 13c. The average brightness of each ISL image and its $1\sigma$ statistical uncertainty are listed in Table 2. Since the ISL at optical wavelengths saturates when the contributions from stars down to 20th mag are added (Leinert et al. 1998), the depth of the Gaia DR3 catalog is sufficient. The photometric uncertainty in the catalog is negligible compared with the uncertainty of the limiting magnitude of ONC-T.

Diffuse Galactic light

Diffuse Galactic light (DGL) consists of starlight scattered by interstellar dust in our galaxy (Elvey and Roach 1937), and it also must be subtracted from the sky brightness to obtain the ZL. A method commonly used to estimate the DGL is to use its correlation with the thermal emission from interstellar dust in the far infrared. The intensity map at $\lambda =$100 $\mu$m, which is a reprocessed composite of the COBE and IRAS maps (SFD map, Schlegel et al. (1998)), is commonly used as a template for the interstellar dust distribution. Thus,

$$\begin{aligned} {\lambda I_{DGL}} = \nu \beta _{\lambda } \cdot d(Glat) \cdot I_{SFD}, \end{aligned}$$

(6)

where $\lambda I_{DGL}$ is the DGL brightness in nW/m$^2$/sr, $I_{SFD}$ is the far-infrared intensity at 100 $\mu$m from the SFD map in MJy/sr, d(Glat) is a geometric function of the Galactic latitude (Glat), and $\nu \beta _{\lambda }$ is the DGL correlation factor in (nW/m$^2$/sr)/(MJy/sr). As discussed in Sano et al. (2016b), this geometric function is given by:

$$\begin{aligned} d(Glat) = d_0 (1-1.1g\sqrt{\sin |Glat|}), \end{aligned}$$

(7)

where $d_0$ is a normalizing parameter, and g is the asymmetry factor of the scattering phase function (Jura 1979).

This DGL correlation factor at optical and near-infrared wavelengths has been derived in many previous studies (Arendt et al. 1998; Witt et al. 2008; Brandt and Draine 2011; Ienaka et al. 2013; Tsumura et al. 2013b; Arai et al. 2015; Kawara et al. 2017; Onishi et al. 2018; Symons et al. 2023), but there are two different results. Recently, the values $\nu \beta _{\lambda } = 3.54 \pm 0.91$ (nW/m$^2$/sr)/(MJy/sr), $d_0 = 1.76$, and $g = 0.61$ have been reported from results obtained by the Long-Range Reconnaissance Imager (LORRI) on New Horizons (Symons et al. 2023). These results were obtained at 10–50 au from the Sun, where the ZL is negligible. The bandpass of LORRI is also similar to that of the wide-band filter of ONC-T (Fig. 1). However, this DGL estimate is about 5–10 times smaller than many previous results. For example, Kawara et al. (2017) found $\nu \beta _{\lambda } = 21.0 \pm 0.9$ (nW/m$^2$/sr)/(MJy/sr) and $d(Glat)=1$ (the geometric function was not considered) at 0.65 $\mu$m based on observations from the Hubble Space Telescope (HST). In the present work, we treat the DGL estimate based on the New Horizons result as a low-level DGL estimate, that based on the HST result as a high-level DGL estimate, and the average of these two DGL estimates as a middle-level DGL estimate. We treat the difference between the low-level and the high-level DGL estimates as a systematic uncertainty.

The spatial resolution of the SFD map (6.1 arcmin) is insufficient relative to our data from ONC-T (22 arcsec). Therefore, we used a far-infrared, all-sky diffuse map based on the AKARI all-sky survey (Doi et al. 2015; Takita et al. 2015) in this study. Because these AKARI all-sky diffuse maps cover wider wavelength ranges with finer spatial resolution and better signal-to-noise ratio than the SFD map, they can serve as a new template for the DGL estimate, replacing the SFD map. In this study, we used the AKARI Wide-S ($\lambda =$90 $\mu$m) map, for which the spatial resolution is $\sim 1.3$ arcmin. One difference between these maps is that point sources have not been removed from the AKARI Wide-S map, whereas they have been removed from the SFD map. Therefore, we masked all the point sources included in the AKARI Far-Infrared Bright Source Catalogue Version 2 (Yamamura et al. 2018). Figure 14 compares the SFD map (left) and AKARI Wide-S map (center) in the field of 2022-08-23, which shows that the AKARI Wide-S map has better spatial resolution than the SFD map.

The publicly available AKARI Wide-S map has a ZL remainder that must be subtracted because only the smooth cloud component of the ZL has been subtracted from the raw data, and other ZL components, such as asteroidal dust bands, have not been subtracted (Doi et al. 2015). The contribution from the unsubtracted ZL components is recognizable in the Wide-S map in the low-ecliptic-latitude region that we study in this work. There is a good linear correlation between the SFD map and the AKARI Wide-S map (Takita et al. 2015), but the AKARI Wide-S map at low ecliptic latitudes shows deviations from the linear correlation owing to the residual ZL that remains to be subtracted. Therefore, we used the data that other ZL components are additionally subtracted from the public AKARI Wide-S map based on a ZL asteroidal-dust-band model (Ootsubo et al. 2016). We confirmed the good correlation between the SFD map and the additionally ZL-subtracted AKARI Wide-S map, as shown in Fig. 14 (right). This correlation is fitted by the equation:

$$I_{\text{SFD}}= a \times I_{\text{Wide-S}} + c,$$

(8)

where $I_{\text{Wide-S}}$ is the far-infrared intensity from the AKARI Wide-S map in MJy/sr, and a and c are the fitting parameters. We obtained $a = 1.54 \pm 0.05$ and $c = -1.20 \pm 0.14$ MJy/sr in our observation fields, which are consistent with the result based on the all-sky data (Doi et al. 2015; Takita et al. 2015). We converted the AKARI Wide-S maps in our observation fields into DGL images DGL(x, y) using equations (6) and (8).

Another issue in the AKARI Wide-S map is the sky coverage. The AKARI all-sky map has $>99 \%$ coverage of the whole sky, but our observation fields contain regions with missing data. We therefore used SFD data for the missing regions in the AKARI Wide-S map. The middle-level DGL brightness based on the AKARI Wide-S map and its statistical and systematic uncertainties are listed in Table 2.

Fig. 14

Diffuse Galactic Light. Left: DGL image for the field of 2021–08-23 based on the SFD map at $\lambda =$100 $\mu$m (Schlegel et al. 1998). Center: DGL image for the same field based on the AKARI Wide-S map at $\lambda =$90 $\mu$m (Doi et al. 2015; Takita et al. 2015) after masking the point sources. Right: Comparison of the intensities from the SFD map and the AKARI Wide-S map of our observation fields. The red line shows the best fit to this correlation

Extragalactic background light

Extragalactic background light (EBL) arises from emissions integrated from the first era of star production to the present day. Recent observations have shown that the EBL measured at optical and near-infrared wavelengths has an excess over the cumulative light from galaxies (Tsumura et al. 2013c; Matsumoto et al. 2015; Sano et al. 2015, 2016a; Matsuura et al. 2017; Mattila et al. 2017; Zemcov et al. 2017; Lauer et al. 2022; Symons et al. 2023; Windhorst et al. 2022, 2023), which means that there are unknown light sources in the universe. The sources for this excess are still under discussion, but some candidates that have been proposed include intra-halo light (Cooray et al. 2012; Zemcov et al. 2014), primordial black holes formed by the collapse of the first halos (Kashlinsky 2016), the decay of hypothetical particles (Kohri et al. 2017), nearby black holes observed as faint compact objects (Matsumoto and Tsumura 2019; Matsumoto 2020), and a warm–hot intergalactic medium (Zhu and Wang 2023). We adopted ${\lambda I_{EBL}} = 21.98 \pm 1.83$ nW/m$^2$/sr at $\lambda$ = 0.44$-$0.87 $\mu$m, as observed by LORRI/New Horizons (Symons et al. 2023). We created an EBL image EBL(x, y) with all pixels having this value.

Zodiacal light

The ZL is obtained by subtracting the background emissions from the observed sky brightness:

$$\begin{aligned} { ZL(x,y) = SKY(x,y) - ISL(x,y) - DGL(x,y) - EBL(x,y) }. \end{aligned}$$

(9)

We first subtracted the ISL image (see "Integrated starlight" section), DGL image (see "Diffuse Galactic light" section), and EBL image (see "Extragalactic background light" section) from the SKY image (see "Flat-field correction" section) to obtain the ZL image, ZL(x, y), and we masked the stars detected in the ZL image (see "Limiting magnitude" section). As shown in Fig. 13d, this masking procedure works well for almost all stars, although the peripheries of some of the brightest stars are not masked perfectly and are smeared out. We then created a histogram of the area where FLAT(x, y) $> 0.5$ and stray light $< 20$ DN (see "Sensitivity calibration using stars" section) for each masked ZL image, and we regard its peak position and its 1$\sigma$ error estimate as the brightness and statistical uncertainty of the ZL (Fig. 15). The systematic uncertainty in the ZL comes from the systematic uncertainties in the calibration and the DGL. Comparing the histograms of the masked ZL image (Fig. 15) with the SKY image (Fig. 12) shows that the excess of the side lobe over the Gaussian has been reduced thanks to the stellar masking. There is still a small excess owing to the smeared peripheries of the brightest stars and to some hot pixels caused by cosmic-ray hits, but this small excess has little effect on the peak position of the Gaussian because we performed the Gaussian fitting on data within a $+2\sigma$ range from the peak (the red dashed line in Fig. 15). The resulting estimate of the ZL and its uncertainties are listed in Table 2 and shown in Fig. 16 (black).

Fig. 15

A histogram of the pixel signals from the masked ZL image of the 2022–10-17 data on a linear scale (left) and a logarithmic scale (right). The red curves show the best-fit Gaussian functions, and the red dashed line shows the $+2\sigma$ distance from the peak position. We performed the fitting using data within this range

Fig. 16

The ZL brightness of each field (black). The red points show the model brightness (Kelsall et al. 1998) extrapolated to optical wavelength using the solar spectrum (Gueymard et al. 2002)

Discussion

Absolute ZL brightness

We have compared the ZL brightness we observed with that predicted by the Kelsall model, which is based on observations of the all-sky ZL brightness obtained from COBE observations at infrared wavelengths (Kelsall et al. 1998). We calculated the ZL model brightness using the ZodiPy code (San et al. 2022), which implements the Kelsall model. The shortest wavelength for which the ZL brightness can be calculated with this model is 1.25 $\mu$m, which is outside the range of the ONC-T wide-band data used in this study. However, the ZL at both 1.25 $\mu$m and optical wavelengths (i.e., the ONC-T wide-band filter) results from scattered sunlight, and its spectral shape is the same irrespective of ecliptic latitude (Tsumura et al. 2010). We therefore determined the ZL brightness at optical wavelengths by extrapolating from the model brightness at 1.25 $\mu$m using the solar spectrum (Gueymard et al. 2002). This extrapolation was performed using the ratio of the solar spectrum at 0.612 $\mu$m and 1.25 $\mu$m. Since the IDP reflectance varies by about 10% between 1.25 $\mu$m and optical wavelengths (Tsumura et al. 2010; Matsuura et al. 2017; Matsumoto et al. 2018), we have assumed a 10% uncertainty in the ZL model brightness at optical wavelengths associated with this extrapolation. Figure 16 compares the observed ZL brightness with the model brightness, which shows that they are consistent with each other within the ranges of uncertainties.

Gegenschein appears in the antisolar direction, and the Kelsall model does not include it. However, our observed fields are shifted from the antisolar direction by $\sim$10 degrees (see Table 1), and the Gegenschein is negligible there (Ishiguro et al. 2013).

There are small excesses of the observed ZL brightness over the model brightness at around 0.9 au, as shown in Fig. 16. This structure may be real, but we cannot confirm this at this point due to the paucity of data points. Verification will require the accumulation of data from future observations.

Field-variance correction

Since the ZL brightness measurements were obtained at different ecliptic latitudes and solar elongations, a correction for the field variance is necessary to compare them under the same conditions in order to obtain the radial profile of the ZL. We performed this field-variance correction based on the Kelsall model. We calculated the seasonal average of the ZL model brightness toward the antisolar direction in the ecliptic plane at various heliocentric distances, as shown in Fig. 17 (left, red). The radial power-law index of this calculated ZL model brightness is $\alpha = 1.34$ because this value is used in the Kelsall model. We compared this ZL brightness toward the antisolar direction with the ZL model brightness toward the fields observed by ONC-T (Fig. 17, left, black). We multiplied the observed ZL brightness by the ratio of these two ZL model brightnesses at each position as correction factors in order to obtain the ZL brightness toward the antisolar direction in the ecliptic plane (Fig. 17, right). These correction factors are listed in Table 2.

Fig. 17

Field-variance correction. Left: ZL model brightness of the fields observed by Hayabusa2# (black) and toward the antisolar direction in the ecliptic plane (red) based on the Kelsall model (Kelsall et al. 1998) at 1.25 $\mu$m. A radial power-law profile with $\alpha = 1.34$ is also shown. Right: Observed ZL brightness after subtracting the background components (black) and the corrected ZL brightness toward the antisolar direction in the ecliptic plane (red)

Since this correction relies on the Kelsall model with $\alpha = 1.34$, it is not self-consistent if the obtained value of $\alpha$ deviates significantly from 1.34. As we discuss in the next subsection, however, the value of $\alpha$ we obtained from our observations is close to 1.34, so consequently this correction works. In addition, the field variance was corrected for local brightness difference at the same elongations but different ecliptic coordinates, so the correction is not sensitive to $\alpha$ which represents the global distribution of IPD.

Dependence on heliocentric distance

Based on the corrected ZL brightness, we calculated the radial power-law index $\alpha$ using the following method. The ZL brightness has two types of uncertainties: statistical uncertainties and systematic uncertainties. Statistical uncertainties appear randomly at each data point, while systematic uncertainties appear with a certain tendency at each data point. Our data contain two types of systematic uncertainties, one due to calibration (see "Sensitivity calibration using stars" section) and the other due to DGL (see "Diffuse Galactic Light" section). We therefore calculated the radial power-law index $\alpha$ for a total of $3 \times 3$ cases (three cases for calibration uncertainty and three cases for DGL uncertainty). Table 3 shows the values of $\alpha$ we obtained for each of these systematic-uncertainty cases. As this table shows, the calibration uncertainty does not have a significant impact on the value of $\alpha$, since the data points only go up and down overall. On the other hand, the DGL uncertainty does have a significant impact on the value of $\alpha$. From all of these values, we obtain $\alpha = 1.30 \pm 0.08$ as the final result. Figure 18 shows the radial profile of the ZL and the best-fit power-law function. Again, the excess structure at $\sim$0.9 au can be seen in Fig. 18.

Table 3 Obtained radial power-law index $\alpha$ for each systematic uncertainty case

Fig. 18

The radial profile of the ZL observed by Hayabusa2# on 2021–2022 at 0.76$-$1.06 au on a linear scale (left) and a logarithmic scale (right). The error bars include both statistical and systematic uncertainties. The solid curve shows the best-fit power-law function ($\alpha = 1.30$), and the dashed curves show the power-law functions for the low-DGL ($\alpha = 1.19$) and high-DGL ($\alpha = 1.42$) cases

Table 4 and Fig. 19 compare the value of $\alpha$ we obtained with previous results, and they show that our result is consistent with them. Since it is difficult to determine the radial profile of the IDP density from ZL observations obtained in a geocentric orbit or from F-corona observations, direct observations of the radial profile of the ZL from interplanetary space outside Earth’s orbit are more reliable for this purpose. Our observations are the first successful observations of ZL from interplanetary space in the 40 years since Helios 1/2 and Pioneer 10/11. In addition, the ZL intensity varies with both the heliocentric distance and the solar elongation, both of which varied in the previous observations by Helios 1/2 and Pioneer 10/11. On the other hand, we confined our observations to the antisolar direction (solar elongation $\sim$180 deg), so we can impose changes in ZL brightness on changes in heliocentric distance.

Table 4 Comparison of the radial power-law index $\alpha$

Fig. 19

Comparison of the radial power-law index $\alpha$ as a function of heliocentric distance. The data plotted here are shown in Table 4. The red data point shows the result from this study

Values of $\alpha$ greater than 1 are obtained from all the ZL observations, even though $\alpha =1$ is expected if the orbital evolution of the IDP is dominated by the PR effect (Burns et al. 1979). This difference is caused by dust production due to the collision of dust particles (Leinert et al. 1983; Grün et al. 1985), dust supplied by comets around 1 au (Ishimoto 2000), the finiteness of the dust cloud (van Dijk et al. 1988), or the heliocentric dependence of the local albedo of the IDP (Giese and Kinateder 1986; Levasseur-Regourd et al. 1991). In fact, a heliocentric dependence of the local albedo of the form $r^{-0.3 \pm 0.1}$ has been reported (Levasseur-Regourd et al. 1991), which partially explains $\alpha$ being greater than 1.

The IDPs falling into the Sun due to the PR effect decrease in size owing to evaporation, and such small IDPs are blown away by radiation pressure as $\beta$-meteoroids. The radial profile of $\beta$-meteoroids is expected to follow a power law with $\alpha =2$, and recent in situ direct counting of the flux of IDPs experienced by the PSP does show $\alpha = 2$ in the range $0.17-0.7$ au (Szalay et al. 2020). On the other hand, the reddening (Tsumura et al. 2010) and polarization (Takimoto et al. 2022, 2023) of the ZL spectrum indicate that the majority of IDPs seen as ZL are large ($>1$ $\mu$m); smaller IDPs do not contribute much to the ZL even though they do exist (Krüger and Grün 2014). While the value $\alpha = 2$ obtained by the IDP impact-counting method is sensitive to small IDPs, the values of $\alpha$ between 1 and 2 obtained from ZL observations are sensitive to larger IDPs. The small dust particles become hotter than the larger ones (Ishiguro et al. 2010), and a hot component in the thermal emission from IDPs has been found using mid-infrared spectroscopy at $\lambda =$3–6  $\mu$m (Ootsubo et al. 1998, 2000; Hong et al. 2009; Tsumura et al. 2013a). We would therefore expect to obtain the value $\alpha \sim 2$ from ZL observations at $\lambda =$3–6  $\mu$m carried out outside Earth’s orbit because small dust particles with high temperatures are mainly observed in this wavelength range.

Future observations

Since ZL observations at 0.7–1 au by Hayabusa2# will continue until the Earth swing-by at the end of 2027, the accuracy of the results we have reported here will be improved through the accumulation of additional observational data. In particular, the excess structure at $\sim$0.9 au needs to be verified by accumulating data during this phase. After the second Earth swing-by in 2028, the Hayabusa2 spacecraft will fly to an orbit in the 1$-$1.5 au range (Mimasu et al. 2022), so we will be able to obtain the radial profile in the outer regions of the Solar System.

Before this orbital change of Hayabusa2#, we will have a chance to observe the radial profile of the ZL in the 1$-$1.5 au range from the Martian Moons eXploration (MMX) spacecraft, which is a Japanese sample-return mission from the Martian satellite Phobos (Kuramoto et al. 2022). The MMX is scheduled for launch in 2024 and arrival at Mars in 2025, and we are proposing to conduct simultaneous multi-wavelength (350–1000 nm) ZL observations during this cruising phase using the Optical RadiOmeter composed of Chromatic Imagers (OROCHI) onboard MMX (Kameda et al. 2021a).

In addition, we are developing an EXo-Zodiacal Infrared Telescope (EXZIT), with the aim of installing it on a spacecraft to Jupiter or farther (Matsuura et al. 2014; Sano et al. 2020). If this instrument can be realized, we will be able to observe the radial profile of the ZL at 1–5 au as well as the EBL without the ZL foreground above 3 au. We are also considering adding mid-infrared capabilities to EXZIT, which would allow us to examine our prediction of the $\alpha \sim 2$ index for the ZL radial profile owing to small particles.

In situ direct dust counting is also important for comprehending the IDP distribution in the Solar System because it provides independent estimates of the radial variation of the IDP density. More quantitative comparisons between the ZL observations and in situ direct dust counting are envisioned for future projects, which will provide us with the differences in the distributions according to dust size and parent bodies.

Summary

We observed the ZL brightness at optical wavelengths at 0.76$-$1.06 au with ONC-T on the Hayabusa2# mission. We detected a small excess of the observed ZL brightness over the model brightness at around 0.9 au, but we cannot determine whether or not this structure is real at this stage due to the paucity of data points. The radial power-law index we obtained is $\alpha = 1.30 \pm 0.08$, and the uncertainty in this estimate is dominated by the uncertainty due to the DGL estimate. This result is consistent with previous results based on other ZL observations.