© M. San et al. 2022

1 Introduction

Zodiacal emission (ZE, sometimes called zodiacal light emission or interplanetary dust emission) is a source of radiation observable from the optical to the submillimeter regimes (Leinert et al. 1998 and references therein). Observations of ZE date back millennia and still hold prominence for a host of modern experiments. ZE is caused by scattering and re-emission at infrared wavelengths of sunlight absorbed by interplanetary dust (IPD) within the Solar System. The diffuse ZE is a major contributor to the total sky brightness at infrared wavelengths between 1 and 100 µm, and it dominates the mid-infrared (~10–60 µm) with a spectral peak at around 10–20 µm.

Zodiacal emission has been a driving force in the exploration of the interplanetary medium since the seventeenth century (Cassini 1685), and it was found to be well characterizable in the infrared by the Diffuse Infrared Background Experiment (DIRBE) instrument onboard the Cosmic Background Explorer (COBE; Mather et al. 1994; Hauser et al. 1998). With the purpose of removing ZE from their data, the DIRBE team constructed a geometric model depicting the interplanetary medium and its distinguishable components. These components were then set up with modified blackbody spectra and evaluated through line-of-sight integration to make predictions of the observed ZE. This model, described in Kelsall et al. (1998) (hereafter K98), has proven to be effective for describing the ZE in the infrared and submillimeter regimes. The DIRBE model was recently used by the Planck Collaboration (Planck Collaboration XIV 2014; Planck Collaboration VIII 2016; Planck Collaboration III 2020) in their analysis of the High Frequency Instrument (HFI) data, where it was adapted to be valid at subterahertz frequencies. The Planck update to the DIRBE model is the state-of-the-art method for removing ZE from data in cosmic microwave background (CMB) cosmology.

In CMB studies, the ZE is a local foreground whose structure is highly dependent on the position of the observer within the Solar System. A consequence of this is that every experiment observes a unique ZE signal due to differences in the scanning strategy and telescope position. As such, it is impossible to describe the ZE foreground through a single template, as would be possible for most Galactic foregrounds, and instead ZE must be dynamically modeled on a per-experiment basis.

While new data will refine our understanding of the interplanetary medium, it is necessary for future models to be consistent with both new and archival data. It is precisely this need that motivates the COSMOGLOBE project¹, which aims to create a framework that will allow for the refinement of astrophysical models jointly with the raw data from complementary experiments. This form of joint analysis is already being explored within the framework of WMAP (Watts et al. 2022) and LiteBIRD (Aurlien et al., in prep.), in combination with Planck Low Frequency Instrument data (BeyondPlanck 2022). ZE is an especially promising direction for joint analysis, in part because HFI made observations of ZE at complementary wavelengths to DIRBE. While the DIRBE model was modified during the HFI analysis, no attempt was made to improve upon the geometrical representation of the model components using the larger effective dataset. Our understanding of the interplanetary medium has improved since the development of the DIRBE IPD model (see, e.g., Reach et al. 1997; Reach 2010a,b). The AKARI satellite, which observed in the infrared at wavelengths between 6 and 180 µm, detected small-scale structures in the ZE which were not well-characterized by the DIRBE model (Ootsubo et al. 2016). While the DIRBE model has been successful in describing the large-scale diffuse ZE, modern high-resolution high-frequency and infrared experiments will require IPD models that can more effectively resolve the small-scale structures in the ZE. As such, an update of the community state-of-art ZE model is long overdue. In this COSMOGLOBE framework, such a model refinement is a natural byproduct of a joint analysis of DIRBE, HFI, and other data. To improve upon this model, it is essential to have the model agree with the data at all wavelengths. The modeling of ZE is in the process of being implemented in the Commander framework, which will ultimately be used for the joint processing of HFI and DIRBE data. However, a stand-alone code is useful for agile model building and data analysis. As such, ZodiPy will function as the first step toward this goal in the COSMOGLOBE framework, in addition to making ZE corrections on arbitrary data more accessible to the community.

In this paper, we present ZodiPy, a modern Python implementation of the DIRBE model in a user-friendly fashion which allows for forecasting of ZE by arbitrary Solar System observers at wavelengths of 1.25–240 µm and frequencies of 30–857 GHz. ZodiPy supports custom user-defined model implementations, and it allows for extrapolations of models to wavelengths or frequencies outside of the valid ranges defined by the models. In addition to forecasting ZE for current and future experiments, ZodiPy will act as a stand-alone code for testing and comparing improved models of ZE with the COSMOGLOBE framework.

In Sect. 2, we introduce the DIRBE IPD model and show geometric illustrations of its IPD components. In Sect. 3, we introduce the emission mechanisms of IPD and describe the implementation of the model evaluation in ZodiPy. Additionally, we present instantaneous maps of the ZE, both the total emission and the component-wise emission. In Sect. 4, we illustrate a few use-cases of ZodiPy by applying the software to make predictions about the ZE foreground as observed by DIRBE.

2 Interplanetary dust model

ZodiPy implements the DIRBE IPD model as described in Kelsall et al. (1998). The model is twofold. The first part is a three-dimensional parametric description of the distribution of IPD in the Solar System, consisting of six individual components: the diffuse cloud, three dust bands, the circumsolar ring, and the Earth-trailing feature. A detailed description of each component is presented in Sect. 2.2. The second part of the model describes the emission mechanisms of IPD. In total, there are over 80 free parameters in the DIRBE model. In the following sections, we introduce the parameterized number density functions of the IPD components.

2.1 Coordinates and geometry

All IPD components, except the Earth-trailing feature, are distributed symmetrically around the Sun with respect to the ecliptic plane. The natural coordinate system for describing the three-dimensional density distribution of an IPD component, c, is Cartesian heliocentric ecliptic coordinates (x, y, z). However, the spatial origin of the components is not necessarily aligned perfectly with that of the Sun. Therefore a coordinate offset is introduced for each component (x_0,c, y_0,c, z_0,c). We define these component-centric coordinates as follows: (1)

In addition to this offset, each component is allowed to have some inclination i and ascending node Ω, with respect to the ecliptic. The symmetry of a component with respect to its mid-plane allows us to completely describe its density with the two coordinates (2) (3)

where R_c is the radial distance from the center of the component, and Z_c is the height above its mid-plane. The ratio between these two coordinates, ζ_c ≡ |Z_c/R_c|, is a helpful auxiliary quantity for describing the fan-like distribution that the diffuse cloud and the dust bands exhibit. We note that all distances in the DIRBE IPD model are given in units of AU and implicitly divided by 1 AU.

2.2 Interplanetary dust components

Below we introduce the DIRBE IPD components. We refer the reader to K98 (and references therein) for an in-depth review of the observational and physical motivations behind the model components.

2.2.1 The diffuse cloud

The primary component in the model is the diffuse cloud. As the name suggests, this is a cloud-like component that encapsulates the inner Solar System and accounts for the majority of total IPD density in the model. The number density of the diffuse cloud, n_C, is modeled in a separable form with a radial and vertical component, (4)

The radial term, is a power law with α as a free parameter. The vertical term, f(ζ_C), represents a widened modified fan model, (5)

where (6)

and β, γ, and µ are free parameters describing the vertical shape and widening of the diffuse cloud. The density distribution of the diffuse cloud, as detailed by the DIRBE model, is illustrated in the top panel of Fig. 1, where we show two cross sections of the tabulated cloud density in the xz and xy planes, respectively.

2.2.2 Dust bands

The model includes three dust band pairs. These are bands of IPD concentrations appearing at distinct ecliptic latitudes associated with asteroid families (Sykes 1990). The density of a dust band, B_j, where j ∈ {1, 2, 3}, is modeled in the following form: (7)

where is the density of the band at 3 AU, is a free parameter representing the inner radial cutoff, and , , and are free shape parameters. We note that there is a factor of difference between this formulation and Eq. (8) in K98. This factor is mentioned in Sect. 4.1.2 in Planck Collaboration XIV (2014) and comes from a difference between the text and the implementation in the COBE DIRBE Zodiacal Light Prediction Software². We follow the convention in the code implementation. The three band pairs, which appear at ecliptic latitudes around ±1.4°, ±10°, and ±15°, respectively, are illustrated in Fig. 2, where we show two cross sections of the tabulated band densities in the xz and xy planes, respectively.

Fig. 1 Interplanetary dust density distribution of the diffuse cloud, the circumsolar ring, and the Earth-trailing feature. Left column: cross section showing the xz plane (y = 0). The gray-dashed line represents the ecliptic plane. Right column: Cross section showing the xy plane (z = 0). The gray-dashed circle represents the orbit of Earth. Top row: the density of the diffuse cloud. Middle row: the density of the circumsolar ring. Bottom row: the density of the Earth-trailing feature. The blue dot represents Earth.

2.2.3 Circumsolar ring and Earth-trailing feature

As the orbit of IPD grains decays, they may end up getting trapped in gravitational resonant orbits with planets, resulting in circumsolar rings. The DIRBE model includes one such circumsolar ring, n_R, for Earth’s orbit, (8)

where n_R is the density of the ring at 1 AU, R_0,R is the radius of the peak density, and σ_R,r and σ_R,z are the radial and vertical dispersions, respectively.

Dermott et al. (1994) showed through dynamical simulations of the circumsolar ring that we should expect enhancements in the dust concentration of the ring in the regions trailing and leading Earth. The DIRBE team was able to fit such a density feature, n_F, to the region trailing Earth, (9)

where θ_F is the ecliptic longitude of Earth seen from the feature-centric frame. The Earth-trailing feature is the only component in the model not symmetrical in its mid-plane, which will require us knowing the heliocentric longitude of Earth at all observation times to determine the position of the Earth-trailing feature. We note that there is a factor 1/2 difference in the first term in the exponential of Eq. (8), and the first and last term in the exponential of Eq. (9) and Eq. (9) in K98, which is again due to a difference between the text and the code implementation. The density distributions of the circumsolar ring and the Earth-trailing feature are illustrated in the middle and bottom panels of Fig. 1, where we show two cross sections of the tabulated densities in the xz and xy planes, respectively.

3 Zodiacal emission

Now that we have a model for the spatial distribution of IPD in the Solar System, the next step is to connect the densities of each IPD component to the foreground emission seen by a Solar System observer. In the following sections, we illustrate the characteristic time variation of the ZE and give an overview of the emission mechanisms of IPD as described in K98. We also show how we apply the full DIRBE model to evaluate the ZE in ZodiPy.

3.1 Time-varying emission

The motion of an observer through the IPD distribution introduces a time dependence in the observed ZE. To an observer orbiting the Sun with Earth, the IPD distribution moves along the ecliptic by approximately one degree a day (with respect to the Galactic coordinate frame). As the observer travels within this distribution, any line of sight toward the same point on the distant celestial sphere contains varying amounts of IPD. This effect is illustrated in Fig. 3.

The characteristic time variation in ZE makes it easy to recognize in data, even at wavelengths and frequencies where the observed emission is very weak. We can extract the ZE signal by splitting the data temporally (with similar sky coverage) and subtracting one from the other, which removes any time-invariant Galactic emission. On the other hand, a downside of this time variation is that the emission becomes a function of the scanning strategy of experiments, making it impossible to construct a universal foreground template of the ZE, as is commonly done for Galactic foregrounds. Instead, the ZE observed by an experiment must be dynamically modeled on a per-experiment basis.

3.2 Emission mechanisms

When heated by their surroundings, dust grains emit thermal radiation. In the case of IPD grains, the heat source is radiation from the Sun. If we assume that each dust grain behaves similar to a blackbody, then the emission from a single dust grain at some delta wavelength, λ, is given by the Planck function B_λ(T), (10)

The temperature of the dust grain, T, is assumed to be a function of the radial distance from the Sun, (11)

where T₀ is the temperature of the dust grain at 1 AU, R is the distance from the Sun, and δ is the power-law exponent describing how the temperature falls with distance from the Sun.

Dust grains come in varying shapes and compositions and they are not true blackbodies. We, therefore, model the IPD grains as modified blackbodies instead, where we scale the Planck function, B_λ(T), by an emissivity, E_λ, which measures the deviation of the IPD grain from the Planck function at some wavelength. Additionally, we assume that all the dust in a specific IPD component is of the same composition, but we allow the composition to vary from component to component, such that the emissivity becomes E_c,λ. The thermal emission then becomes (12)

In addition to generating thermal radiation, IPD grains also scatter sunlight at infrared wavelengths. The scattering term, as defined in K98, is (13)

Here A_c,λ is the albedo, which represents the fraction of photons that get reflected by the surface of the grain and is a function of the grain composition, and is the solar spectral irradiance, representing the flux received at some wavelength by the surface of the dust grains at some distance from the Sun. The spectral solar irradiance is given by (14)

where is the reference spectra of the solar irradiance at 1 AU, and R is the radial distance from the Sun. The solar irradiance reference spectra used in the DIRBE model is tabulated in the COBE DIRBE Zodiacal Light Prediction Software.

Finally, Φ_λ(Θ) is the phase function and represents the distribution of scattered light intensity. The phase function as defined in K98 is (15)

where C_k,λ are free parameters. The phase function can be interpreted as a probability density, representing the probability of a photon getting scattered as a function of the scattering angle Θ, and it must therefore be normalized. The DIRBE normalization factor N has the following analytic solution: (16)

At scattering wavelengths, that is to say wavelengths where A_c,λ ≠ 0, we assume some fraction of the total radiated heat to be scattered away from the observer, such that the thermal contribution becomes (17)

Combining Eqs. (13) and (17) gives us the total emission in intensity from a single grain of IPD, (18) (19)

Fig. 2 Interplanetary dust density distribution of the three dust bands. Left column: cross section showing the xz plane (y = 0). The gray-dashed line represents the ecliptic plane. Right column: cross section showing the xy plane (z = 0). The gray-dashed circle represents the orbit of Earth. Top row: the density of dust band 1. Middle row: the density of dust band 2. Bottom row: the density of dust band 3.

Fig. 3 Illustration showing that the integrated IPD along a line of sight toward a point on the celestial sphere as seen from Earth (blue circles) changes as Earth orbits the Sun (yellow circle).

3.3 Evaluating the zodiacal emission

Now that we have expressions for both the number densities of the IPD components and the emission mechanisms, we can turn to the full model evaluation in ZodiPy. For an experiment observing some point p on the sky at time t, the total contribution to the sky intensity from ZE is (20)

where we integrate along a line-of-sight ds from the observer toward p. Equation (20) represents the brightness integral as described in Eq. (1) in K98, but without the color-correction factor, and it is further reduced to Eq. (8) in Planck Collaboration XIV (2014) at wavelengths where A_c,λ = 0. We note that in the DIRBE analysis, the endpoints of all considered line of sights were fixed to a radial distance of 5.2 AU from the Sun, corresponding to the orbital distance of Jupiter. We would prefer to relax this assumption, but it is likely that the DIRBE IPD model parameters are affected by this cutoff. As such, ZodiPy follows this convention by default, but the user may specify an alternative cutoff distance. The line-of-sight integration is performed using Gaussian-Legendre quadrature, with a default of 100 points, but the user may change this number depending on the desired accuracy.

The pointing in Eq. (20) may be given as a single sky coordinate or a sequence of sky coordinates, either in the form of angles on the sky (θ, ϕ) or as integers referring to the indices of pixels on a HEALPix³ pixelization grid with the resolution given by N_side (Górski et al. 2005). The pointing information is internally converted to ecliptic unit vectors on the sphere using the healpy software (Zonca et al. 2019) and further used (alongside the observer position and the radial cutoff) to construct the evaluated line of sights. The observer position is assumed to be constant over the duration of the user-provided input pointing sequence, which means that the pointing must be appropriately chunked depending on the position and velocity of the observer. For instance, observers orbiting the Sun alongside Earth should provide pointing sequences of periods corresponding to a day of observation at most, so as to account for the time variations in the ZE. The predicted ZE corresponding to a pointing sequence may be outputted as either a timestream or as a binned HEALPix map, both in units of MJy sr⁻¹. For code examples of how to use ZodiPy, we redirect the reader to the documentation⁴.

3.4 Instantaneous emission maps

In this section, we explore the structure of the ZE as seen on a particular day with the DIRBE model. We make maps of the full-sky ZE in both ecliptic and Galactic coordinates and full-sky maps in ecliptic coordinates of each IPD component, respectively. These maps illustrate how the parametric DIRBE model relates to the observed ZE and functions as a validation of the model implementation. Figure 4 shows the instantaneous emission at 25 µm seen by an observer on Earth at two different times (January 1, 2020, in the top panel and June 1, 2020, in the bottom panel). These are maps that represent the ZE as seen by an observer able to simultaneously view the full sky at one instant in time.

The emission can also be viewed component-wise, as illustrated in Fig. 5, which shows the emission from Fig. 4 split into the respective model components in ecliptic coordinates. The structure observed in these maps can be better understood by inspecting them alongside the density distributions seen in Figs. 1 and 2. The diffuse cloud is the brightest IPD component at this wavelength, and as such, it looks very similar to the total ZE emission. We see that dust bands 1 and 3 appear wider on the side of the map facing away from the Sun. This widening is due to more high-latitude line of sights looking through the bright parts of the bands in this direction. The subtle monopole present in the map of dust band 2 results from the observer being encapsulated by dust, which we can see from the density distribution in the middle panel of Fig. 2. The asymmetry of the emission with respect to the ecliptic is due to the inclination of the bands.

The circumsolar ring and the Earth-trailing feature exhibit perhaps the most distinct structures in the maps. The bright stripe in the ecliptic wrapping behind the Sun in the map of the circumsolar ring (bottom left of Fig. 4) can be understood as the observer looking past the Sun and through the opposite side of the circumsolar ring. The two bright regions symmetrically around the Sun in the ecliptic represent the observer looking tangential to Earth’s orbit along the longest sections of the ring. We can understand the partial sky brightening in the circumsolar ring emission by noting that the observer is located at the inner edge of the ring (see the gray circle in the middle-right panel of Fig. 1). This means that all directions pointing away from the Sun look through some portion of the circumsolar ring. A similar signature is observed in the Earth-trailing feature map (bottom-right of Fig. 4; however, in this case, only one such bright region appears on the side pointing toward the feature. Since the Earth-trailing feature is much wider than the circumsolar ring in the xz plane, the observer is again encapsulated by the dust in this component which results in another monopole. Satellite missions that observe from the second Sun-Earth Lagrange point (L2, R ≈ 1.01 AU) should expect a slightly higher contribution from the circumsolar ring and Earth-trailing feature due to the observer being located closer to the denser regions of these components.

Fig. 4 Instantaneous maps of the predicted ZE at 25 µm as seen by an observer on Earth on January 1, 2020 (top panel) and June 1, 2020 (bottom panel), in ecliptic (left column) and Galactic (right column) coordinates. The location of the Sun is marked as an orange circle on all four maps.

Fig. 5 Instantaneous emission maps (January 1, 2020) for each IPD component in the DIRBE model. The Sun’s location is marked with an orange circle. Top left: the diffuse cloud. We note that the intensity is logarithmic due to the highly exponential nature of the zodiacal cloud density near the Sun. Top right: dust band 1. Middle left: dust band 2. Middle right: dust band 3. Bottom left: the circumsolar ring. Bottom right: the Earth-trailing feature.

4 Comparison with DIRBE

In the following section, we illustrate a few use-cases of ZodiPy and validate the code implementation by applying the software to produce simulated timestreams using the pointing of DIRBE calibrated individual observations (CIO)⁵. We show that ZodiPy can accurately reproduce the ZE simulations created with the COBE-DIRBE Zodiacal Light Prediction Software used by the DIRBE team in their analysis of the experiment. Additionally, we use ZodiPy to produce binned simulated maps of the ZE at 25 µm as seen by DIRBE (using the CIO pointing) and compare these to binned maps of the DIRBE CIO.

4.1 Timestreams

The user input required to produce a timestream with the DIRBE Zodiacal Light Prediction Software is (1) a wavelength corresponding to one of the ten central band wavelengths of DIRBE; (2) the day of observation represented as an integer counting the number of days since January 1, 1990; (3) angular pointing information; and (4) an optional flag for whether or not to color-correct the output. In comparison, ZodiPy expects the following input arguments: (1) any wavelength or frequency; (2) the name of the observer (defined below) or explicit observer position in units of AU; (3) the time of observation; and (4) pointing information in the form of angles or HEALPix pixel indices. Since ZodiPy can produce timestreams for the arbitrary observer, the observer position is required explicitly. This argument is provided either by specifying the name of the observer, for which ZodiPy computes all required Solar System positions (given the time of observation) using the Astropy Solar System Ephemerides (Astropy Collaboration 2013, 2018) or by explicitly passing in the heliocentric Cartesian coordinates of the observer in units of AU. In order to compare the two codes, we need to consider the location of DIRBE explicitly. The COBE satellite carrying the DIRBE instrument made observations from a near-Earth polar orbit at a distance of about 900 km above the surface. The variation in DIRBE’s position from its orbit is insignificant compared to the scale at which the IPD is distributed. As such, we approximate DIRBE’s location to be the center of Earth.

The comparison of both codes is shown in Fig. 6 where we have generated timestreams over a small subsequence of the pointing from DIRBE’s first day of observations. The timestreams are simulated at 25 µm, which corresponds to the central wavelength of DIRBE’s band 6, which is the band that observed the most ZE. The top panel shows the ZE timestream as predicted by ZodiPy, and the bottom panel shows the difference in percentage between the ZodiPy and DIRBE software predictions. We find the difference between the two codes to be less than 0.2%. A similar comparison was also conducted for the other DIRBE bands, with results consistent with those presented here. We believe that the minor difference in the predictions made by the two codes may be attributed to one of or a combination of the following: (1) a difference in how the input pointing is mapped to unit vectors on the sphere; (2) numerical integration performance for evaluation of the line of sights; (3) a difference in the computed Solar System position; and (4) implementation of the Planck function.

Fig. 6 Zodiacal emission timestreams. Top: timestream simulated with ZodiPy of the ZE as predicted to have been observed by DIRBEs band 6 detector A on December 11, 1989 between ~13:30 and ~14:30 UTC. Bottom: error in percentage between the ZodiPy timestream in the top panel and a corresponding timestream produced using the DIRBE Zodiacal Light Prediction Software.

Fig. 7 Top: binned week maps of the DIRBE CIO. Middle: binned week maps from the ZodiPy simulations. Bottom: DIRBE minus ZodiPy binned week map residuals. The columns correspond to week 1 (left) and the full-survey (right), respectively.

4.2 Binned emission maps

Figure 7 shows two instances of ZE corrections applied to the binned DIRBE CIO, which are comparable with Figs. 2c and 7 in K98. In the left column, we see the corrections applied to the first week of observations, and in the right column, we see the corrections applied to the full mission.

The top, middle, and bottom panels show the binned CIO, the corresponding ZodiPy predictions, and the CIO minus the ZodiPy predictions, respectively. We note the stripes and discontinuities in the full-survey maps. These are artifacts of the binning process and not features of the infrared sky. More specifically, these stripes correspond to differences in the path through the IPD distribution, as discussed in Sect. 3.1, where adjacent positions along the scanning pattern on the sky have been observed over extended time intervals.

We note that the simulations in this comparison were produced using a modified version of the code that included the color-correction factor in Eq. (1) in K98. This was done to make the simulated data more comparable to the CIO, which are color-corrected. The presented residual maps exhibit some weak digitization-like structures. We believe that this pattern is due to errors related to the conversion of the DIRBE pixel indices (natively given in the COBE quadrilateralized spherical cube format) to longitude and latitude rather than an inherent property of the CIO or the simulated emission.

5 Conclusions

We have presented software for modeling the thermal emission and scattered sunlight by IPD, as defined by the DIRBE and Planck models, for arbitrary Solar System observers at both infrared wavelengths and subterahertz frequencies. The software is open-source and available on the Python Packaging Index. The code is designed to allow for new IPD model additions and updates with new developments to state-of-the-art models. Finally, as a validation of the code implementation, we have demonstrated that the software can accurately reproduce the Zodiacal corrections applied in the official analysis of the DIRBE data.

Currently, ZodiPy only models ZE in intensity. Polarization from thermal emission from ZE has never been detected, and it is expected to be less than 1% (Weinberg & Hahn 1980; Onaka 2000; Ganga et al. 2021). The scattered sunlight is known to be polarized as a function of the scanning strategy with the scattering angle and elongation, as shown in Berriman et al. (1994) and Takimoto et al. (2022), for example. Allowing the user to predict the polarized emission at wavelengths where the ZE is dominated by scattered sunlight (1–5 µm) would help with forecasting for future polarized experiments. For a more complete discussion of the intrinsic polarized emission of ZE, readers can refer to Ganga et al. (2021) and references therein.

ZodiPy is an integral part of the COSMOGLOBE framework. Part of COSMOGLOBE’S goal is to perform end-to-end analysis using all microwave and infrared datasets to constrain the properties of the sky and the large-scale properties of the Universe. This code will be an invaluable tool for characterizing the ZE in archival datasets, such as COBE/DIRBE and Planck/HFI, as well as future data from JWST, SPHEREx, and LiteBIRD. A joint analysis will not only help to better decontaminate data with ZE emission, but the joint analysis approach taken by COSMOGLOBE will necessarily improve the IPD model from K98. In combination with the full Gibbs sampling framework with Commander3, ZodiPy will be used to develop new models and visualize the individual components in an intuitive manner, all while providing the astronomical community at large with a simple interface to predict ZE in their own data.

Acknowledgements

We thank Profs. Hans Kristian Eriksen, Ken Ganga, and Ingunn Wehus for useful comments. We also thank the referee for their useful comments that improved this work. The current work has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement numbers 819478 (ERC; COSMOGLOBE) and 772253 (ERC; BITS2COSMOLOGY). Some of the results in this paper have been derived using the HEALPix (Górski et al. 2005) package. We acknowledge the use of the Legacy Archive for Microwave Background Data Analysis (LAMBDA), part of the High Energy Astrophysics Science Archive Center (HEASARC). HEASARC/LAMBDA is a service of the Astrophysics Science Division at the NASA Goddard Space Flight Center.

References

All Figures

Fig. 1 Interplanetary dust density distribution of the diffuse cloud, the circumsolar ring, and the Earth-trailing feature. Left column: cross section showing the xz plane (y = 0). The gray-dashed line represents the ecliptic plane. Right column: Cross section showing the xy plane (z = 0). The gray-dashed circle represents the orbit of Earth. Top row: the density of the diffuse cloud. Middle row: the density of the circumsolar ring. Bottom row: the density of the Earth-trailing feature. The blue dot represents Earth.

In the text

Fig. 2 Interplanetary dust density distribution of the three dust bands. Left column: cross section showing the xz plane (y = 0). The gray-dashed line represents the ecliptic plane. Right column: cross section showing the xy plane (z = 0). The gray-dashed circle represents the orbit of Earth. Top row: the density of dust band 1. Middle row: the density of dust band 2. Bottom row: the density of dust band 3.

In the text

	Fig. 3 Illustration showing that the integrated IPD along a line of sight toward a point on the celestial sphere as seen from Earth (blue circles) changes as Earth orbits the Sun (yellow circle).
In the text

	Fig. 4 Instantaneous maps of the predicted ZE at 25 µm as seen by an observer on Earth on January 1, 2020 (top panel) and June 1, 2020 (bottom panel), in ecliptic (left column) and Galactic (right column) coordinates. The location of the Sun is marked as an orange circle on all four maps.
In the text

Fig. 5 Instantaneous emission maps (January 1, 2020) for each IPD component in the DIRBE model. The Sun’s location is marked with an orange circle. Top left: the diffuse cloud. We note that the intensity is logarithmic due to the highly exponential nature of the zodiacal cloud density near the Sun. Top right: dust band 1. Middle left: dust band 2. Middle right: dust band 3. Bottom left: the circumsolar ring. Bottom right: the Earth-trailing feature.

In the text

	Fig. 6 Zodiacal emission timestreams. Top: timestream simulated with ZodiPy of the ZE as predicted to have been observed by DIRBEs band 6 detector A on December 11, 1989 between ~13:30 and ~14:30 UTC. Bottom: error in percentage between the ZodiPy timestream in the top panel and a corresponding timestream produced using the DIRBE Zodiacal Light Prediction Software.
In the text

	Fig. 7 Top: binned week maps of the DIRBE CIO. Middle: binned week maps from the ZodiPy simulations. Bottom: DIRBE minus ZodiPy binned week map residuals. The columns correspond to week 1 (left) and the full-survey (right), respectively.
In the text