Improved large scales interstellar dust foreground model and CMB solar dipole measurement

The Cosmic Microwave Background anisotropies are difficult to measure at large angular scales. In this paper, we present a new analysis of the \Planck\ High Frequency Instrument data that brings the cosmological part and its major foreground signal close to the detector noise. The solar dipole signal, induced by the motion of the solar system with respect to the CMB, is a very efficient tool to calibrate a detector or a set of detectors with high accuracy. In this work, the solar dipole signal is used to extract corrections of the frequency maps offsets reducing significantly uncertainties. The solar dipole parameters are refined together with the improvement of the high frequency foregrounds, and of the CMB large scales cosmological anisotropies. The stability of the solar dipole parameters is a powerful way to control the galactic foregrounds removal in the component separation process. It is used to build a model for Spectral Energy Distribution spatial variations of the interstellar dust emission. The knowledge of these variations will help future CMB analyses in intensity, and also in polarization to measure faint signal related to the optical reionization depth and the tensor-to-scalar ratio of the primordial anisotropies. The results of this work are: improved solar dipole parameters, a new interstellar dust model, and a large scale cosmological anisotropies map.


Introduction
Measuring accurately the diffuse extragalactic emissions at low frequencies on the nearly full sky is difficult, especially at large angular scales. It indeed implies being able to separate the atmospheric and the telescope emissions from the astronomical ones and, among those, to separate the extragalactic emissions from the galactic diffuse background ones. The Cosmic Microwave Background (CMB) is the truly diffuse background bringing key cosmological information. At the the peak emission of the CMB frequencies, it imposes to go to space, build cryogenic telescopes and instruments, and find compromises between constraints of the size of telescopes giving a limited angular resolution, and the necessity to cover the whole sky.
In the last 25 years, there has been three generations of CMB space missions COBE, WMAP and Planck, the next one not providing results before a decade. The Planck High Frequency Instrument (HFI) observed the sky (at 100, 143, 217, 353, 545, and 857 GHz) with 50 bolometers cooled at 100 mK. The present and future observations from the ground use large detector arrays, but strongly affected by the atmospheric emission above 200 GHz. Future space missions will allow to map the CMB and the galactic foregrounds at higher frequencies. Progress of experiments on the ground will certainly outperform Planck capabilities on the CMB itself up to 150 GHz but not on the high frequency galactic foregrounds, dominating signal at higher frequency.
The large scale polarized B modes is a major goal for cosmology which requires very high sensitivity at CMB frequencies (through more detectors and better control on polarized systematic effects). It also requires better measurements of the low frequency foregrounds with large ground-based surveys, and of high frequency foregrounds hard to achieve from the ground. It is thus a key objective to squeeze out of the HFI data the information on the large scale high frequencies foregrounds. This is the goal of this work.
Although HFI achieved the remarkable result of being cosmic background photon noise limited at the CMB peak frequency at the intermediate and small angular scales (a few degrees to 5 ), at very large angular scales ( < 15), the control of instrument systematics and foreground residuals took a long time to reach the needed level to measure of the cosmological parameter, the reionization optical depth, with a signal to noise ratio of 10 (Pagano et al. 2020). The noise and systematic effects residuals at the largest angular scales are nevertheless still far from the detector white noise level needed to accurately remove the high frequency foregrounds. Removing those from the data of a differential detector requires taking advantage of the solar system kinetic dipole. The solar kinetic dipole, common to all detectors and frequencies, is the strongest measurable differential signal at large scales. Improving its determination is thus a key issue throughout this paper. Importantly, the solar kinetic dipole is extracted from various splits of the HFI data (frequencies, galactic masks, component separation models), which constrains an high frequency galactic foreground model. This provides a tool for intercalibration of different instruments and missions as well as constraining a better model of the high frequency foreground systematics at large scales. These improvements, not achieved in the Planck Legacy release, are introduced in this work where we develop a new interstellar dust model (the dominant foregrounds for the CMB above 100 GHz), self-consistent with a better CMB anisotropies and the associated better solar dipole measurement.
The paper is organized as follow: Sect. 2 describes the context and goals of the paper; Sect. 3 describes the algorithm to extract coherently the parameters and templates maps of the

Context and definitions
The relative motion between the Planck spacecraft and the CMB rest frame generates, a dipolar modulation of the CMB intensity field. The spacecraft motion with respect to the cosmological background is decomposed in two components.
The first component, referred to as the "orbital dipole", is caused by the relative velocity of the spacecraft at the L2 Lagrange point following the earth motion with respect to the solar system rest frame. It is known with a high accuracy at any time. This orbital dipole of order 30 km s −1 does not project on the sky maps. It is used to calibrate the absolute photometric response of each detector with a high signal to detector noise of order 10 −5 at the 100 and 143 GHz CMB dominated channels, and 2×10 −4 at 353 GHz (see Planck Collaboration III 2020, hereafter HFI-Legacy).
The second component, referred to as the "solar dipole", is the consequence of the relative velocity of the solar system with respect to the cosmological rest frame. Although it is an order of magnitude larger than the orbital dipole, it cannot be determined from extragalactic astronomy, and thus cannot help to calibrate the absolute response of the detectors. The dipole signal associated with this second component projects onto the sky maps with the same Spectral Energy Distribution (SED) as the CMB cosmological anisotropies. As such, it is not a cosmological effect, but is a systematic effect in the maps extracted by component separation methods on the basis of their SED. It is, nevertheless, not separable from the intrinsic dipolar term of the CMB primordial anisotropies. Using a ΛCDM model, the CMB anisotropies measurement leads to an expected dipole amplitude of the order of 30 µK. This induces an uncertainty of the order of one percent on the kinetic solar system velocity. This dipolar term is thus included in the definition of the solar dipole, imposing a null dipolar term for the total CMB anisotropies maps.
The solar system velocity measurement is independent of the CMB monopole temperature, known with an accuracy of 3×10 −4 (Fixsen 2009). This introduces additional uncertainty when expressing the intensity amplitude of the solar dipole in δT CMB anisotropies. Nevertheless, for convenience, all intensity contributions to the dipole amplitudes are expressed in CMB temperature. As such, the conversion between velocity and µK CMB depends directly on the CMB temperature, obtained by fitting a Planck function in the COBE-FIRAS data, which has been refined using a combination of the solar dipole direction from WMAP and the FIRAS full spectrum measurements (Fixsen 2009). In this regard, using the solar dipole to estimate the CMB monopole temperature makes to conversion from velocity to µK CMB slightly dependent of the dipole. However, the difference induced by a new dipole determination leads to a negligible correction.
The amplitude of the solar dipole signal is proportional to the response of each detector to a CMB intensity signal, and thus, very useful to test CMB relative photometric calibration of different CMB detectors. Its amplitude of order 3 000 µK CMB is very large with respect to the sensitivity of the 143 GHz map (amplitude of the noise dipolar < 10 −2 µK). This 3 × 10 5 ratio has been used to measure the relative CMB photometric calibration on single detector maps. The consistency of the relative CMB photometric calibration on single detector maps within a frequency has a noise limited accuracy of 10 −5 . The intercalibration between HFI frequency maps is at the 4 × 10 −4 level from 100 to 353 GHz (HFI-Legacy). This work uses these sensitivities, and exploit the fact that all detectors measure the same CMB solar dipole (called hereafter "the consistency argument") to control the residuals from foregrounds on the largest angular scales.
The CMB map is obtained from the frequency maps using its properties: Planck function SED and Gaussian statistics. The extraction of solar dipole amplitudes measured at various frequencies and on different sky fraction is the most powerful tool to detect the foreground residuals. This was used for example for the HFI data showing residuals behaving with frequency like a galactic dust component as can be seen in Tables 7 and 8 of the HFI-Legacy are obviously due to the dust foreground removal, and showed that introducing large scale dipolar and quadrupolar terms to account for the large scale spectral energy distribution (SED) variation reduces very much these drifts for the highest frequencies and sky fraction used to fit the solar dipole.
Here, the solar dipole is extracted with an improved version SRoll2 (Delouis et al. 2019) of the SRoll algorithm used to produce the Planck Legacy HFI maps. This upgraded algorithm, named SRoll2.2, introduces improvements: -Use of 143, 217, 353, 545, and 857 GHz non-polarized Spider Web bolometers (SWB) data, in addition to Polarized Spider bolometers (PSB) data. This allows for reduced systematic effects for intensity only maps, while avoiding being limited by signal to noise ratios; -Removal of polarized parts of the signals due to the small parasitic sensitivity to polarization of SWBs detectors; -Improved correction of the Analog-to-Digital Converters non-linearities, bringing the associated residuals below detector noise levels, even at the largest angular scales; -Computation of single detector maps removing polarization part in the signal; -Use of the full resolution for the bandpass leakage correction where the signal is very strong, keeping the low resolution of the foreground templates on the rest of the maps to avoid introducing extra noise; -Removal the Cosmic Infrared Background (CIB) monopole (Lagache et al. 1999) and use CMB maps with zero monopole; -Setting the initial zero level of the intensity maps to the zero level of the 21-cm emission of the H i galactic gas foreground.
SRoll2.2 builds two types of single detectors intensity maps 1 where: i) the correction of all foreground signal bandpass mismatch between detectors is applied, ii) no correction of the foreground signal bandpass mismatch is applied. Those latest maps, smoothed at a common 1 • , are used in this work.

Foreground sky model and zero level offsets
The solar dipole amplitudes of Planck Collaboration Int. XLVI (2016) Table 2, are drifting away from the amplitude measured at 100 and 143 GHz both for higher frequencies (217 to 545 GHz), and when increasing the sky fraction f sky used to fit the dipole. The increase of the amplitude errors with frequency follows a typical galactic dust SED. Such variations were expected from previous work on the COBE-DIRBE data (Lagache et al. 1998). HFI-Legacy shows that this is mostly cured by introducing, in the model, large scale variations of the dust SED. This result is obtained using the stability of the solar dipole direction, taking the 100 GHz one as reference. The dipole direction (especially the longitude) drifts by 0.1 • but the parameters for the 545 GHz still exhibits large shift in amplitude and direction.
These results indicate that the residuals from the dust emission removal are a major concern when extracting the solar dipole above 100 GHz in two ways: first, indirectly by the dust residuals along the galactic ridge (central parts versus anticentre of the galactic disc). This forces the introduction of a galactic ridge mask M g where the the CMB anisotropies cannot be removed accurately enough. This induces an unavoidable error when the dipolar term of the cosmological anisotropies is set to zero on the full sky (see Sect. 4.2); second, the large scale galactic emission in the part of the sky where the solar dipole is fitted outside the galactic ridge mask affects the residuals of the component separation process. We thus introduce, in our foreground dust emission sky model, terms improving the modeling of the large scales SED variations.
The CMB anisotropies component of the maps, as well as the noise and systematic effects, are defined with a null monopole term. The HFI frequency maps do not have an absolute zero setting calibration. Any error in the zero level changes the contrast between the high and low galactic latitude foregrounds, thus changes the foreground maps, especially where the intensity is very low at high latitude. The large scale residuals of dust emission depends strongly on the accuracy of these zero levels, and we introduce consistency constraints from the solar dipole to refine the accuracy of the foreground model.
To start with, the initial setting of the zero level of the Planck Collaboration HFI frequency maps use the method in which the extrapolated zero intensity emission of the galactic foregrounds gas correspond to the extrapolated zero column density of gas (Planck Collaboration VIII 2014). The first step is performed by regressing the 857 GHz, fully dominated by the dust emission intensity, to the neutral hydrogen interstellar gas column density. Then the other frequencies are regressed to the 857 GHz emission brightness. This sets the initial zero levels of the frequency maps good enough accuracy as long as one does not require highly accurate intensity on very large scales.
To improve on these results, and especially on the dust template map zero level, a more accurate method is developed here 2 .

Dust model definition
We decompose the intensity sky model as: In Eq. 1, there is only one extragalactic component, the CMB anisotropies. The CIB large scales ( < 100) anisotropies are negligible with respect to the galactic dust component, and are included in it. CMB tot is the total CMB anisotropies map containing the solar dipole and the primary cosmological anisotropies. The CMB tot intensity is described by a single map of δT CMB at all frequencies.
The three next terms of Eq. 1 are the low frequency foregrounds: free-free, synchrotron emission, and CO emission lines. In intensity, the free-free and synchrotron emissions are respectively one, and two orders of magnitude smaller than the CMB anisotropies at 100 GHz. The free-free emission, which is only one order magnitude lower than the dust at 100 GHz, benefits from a stable and accurately known SED. Furthermore, the synchrotron emission is lower than other emissions by orders of magnitude at higher frequencies. Thus free-free and synchroron emissions are removed in an open loop, using Planck Legacy results (Planck Collaboration X 2016). In the HFI data, the CO lines emission are derived using only the large differences of response between detectors within a frequency band containing a CO line and the Planck Legacy CO template maps are used. Nevertheless, their amplitude are updated, consistently with the new dust model.
The last term of Eq. 1, the interstellar dust emission, is the dominant foreground at frequencies above 100 GHz. It is described by a new model: where f (ν) = Planck function(ν, T ) × ν ν 0 β , ν 0 beeing the reference frequency of the Dust map.
The temperature T and the index β are very correlated when fitted in the range of frequencies used in this work, thus the temperature is taken as constant, and the SED variations are only described by the variations of the β parameter. Given the available number of degrees of freedom (5 frequencies), this has been shown acceptable in the HFI frequency range from previous work using a broader frequency range from the absolute intensity maps from COBE. The single temperature adopted 18 K has been shown to give good SED fits for the H i dominant phase of the interstellar medium but also for the subdominant H +, that show very similar SED (Lagache et al. 1998).
The novelty of this dust model is the use of an analytical description of the averaged SED, and the introduction of two additive correction modeling the effects of the SED spatial variations.
The first term of Eq. 2, Dust, is the initial dust map (see further), propagated from the frequency ν 0 to each single detector map with a frequency ν. Furthermore, the frequency is taken at the effective frequency that takes into account the difference of the response of the detector to a dust SED with the response to the CMB orbital dipole signal used to calibrate the detector with high accuracy in K CMB .
The two other terms of Eq. 2 describe the spatially variable part of the SED, absorbing all differences to the spatially constant SED first term. The analytical description of the SED f (ν) allows a description of these variations as an expansion in term of its first and second derivatives with respect to β. The model is not strictly described by the analytical expansion of the spectral energy distribution function f (ν) because the second spatial template ∆Dust 2 is independent from the first one, ∆Dust 1 , which describes the first order expansion term. These two independent full sky maps with a null average have a resolution of order 1 • (HEALPix (Górski et al. 2005 The initial dust map Dust is built from a map at high enough frequency (ν 0 =545 GHz) that it only needs to be cleaned from the CMB, the CO, and the free-free emissions in an open loop. The CMB large scales cosmological anisotropies map, taken from the Planck Legacy release, once propagated to ν 0 , are then three orders of magnitude smaller than the dust emission. The effects of the difference between the initial CMB used at this stage with the CMB map obtained ultimately in this work, is thus negligible.
From Eq. 2, we obtain the dust map model D(ν b ) for a single detector (bolometer b) map: This equation is valid for maps expressed in spectral energy density. The 545 and 857 GHz maps are calibrated on planets and are dealt with in spectral energy density expressed in MJy sr −1 . Nevertheless, 100-353 GHz detectors, calibrated on the orbital dipole of the CMB in δT CMB , are expressed in µK CMB , and must be converted to spectral energy density. The lower frequencies CMB maps, offsets and solar dipole amplitudes are finally be converted back to µK CMB for the convenience of comparison with similar maps in the literature. We know that SED variations exist at intermediate to large scales due to the gradients of i) stellar radiation field intensity changing the dust temperature, ii) the relative fraction on each line of sight of the molecular, atomic and ionized gas fraction which are known to have slightly different dust properties. The two SED variation components are thus defined by maps smoothed at 1 • .

Extraction of the CO emission lines, dust model parameters and templates
The SRoll2.2 single detector intensity maps are built after removing the polarized signal and the bandpass mismatch coefficients have been computed using template dust maps different from the ones build in the present work. We thus use the set of maps in which the color correction has not been corrected to be able to extract band pass mismatch coefficients coherent with the dust model using an iterative process. We name I b the single detector intensity map built from the bolometer b data. The effective frequency ν b , computed from the bandpass ground measurements, is used in the first iteration of the processing, and is corrected at each iteration when regressing the single detector map against the dust initial template Dust ν 0 .
The dust emission is the only broad band foreground component not fully described by its SED, but modeled empirically for the frequencies ν ≥ 100 GHz as a map with an average SED, completed with two additive correction maps describing the spatially variable SED. These three maps are extracted in this work, using other constraints than previous component separation ones as discussed above (e.g. consistency of the solar dipole parameters). The CO and dust foregrounds being partially correlated, the CO foreground map is built using a CO template map from the Planck Legacy results (Planck Collaboration X 2016) but whose coefficients are extracted coherently with our new dust model.
Finally, the data contains noise and monopole offsets, not directly measurable by the HFI differential detectors, and which are to be solved for, using consistency arguments and redundancies.
Equation 4 defines the residual R b which only contains the CMB tot and the foreground components to be modeled. It is computed for each single detector map I b from which are removed, in an open loop, the low frequency foreground maps taken from the Planck legacy maps as described above.
The initial offset O HI applied to the single detector map sets the zero level of the single detector intensity maps I b , following Planck Collaboration VIII (2014) that sets this zero level to the extrapolation to zero column density of HI 21-cm emission of interstellar gas at high latitudes. This is a first approximation of the absolute zero level of the dust foregrounds. The addition of the ionized component traced by H α emission changes the zero level by 3% within the uncertainties, and the molecular component is negligible at high galactic latitudes. The offsets are further refined in this work with smaller uncertainties, and the final results should be consistent with the initial offsets. The offset values for detectors, either within the same frequency band or with different frequencies, change the contrast of the large scales at high latitude but are constrained by the solar dipole signal which is a common component to all bolometers. For the 100 to 353 GHz detectors, the single detector map offset O b is adjusted in three steps: 1. The first step uses all pair differences that gives offsets with respect to a zero average offset of all these bolometers. 2. The second step adjusts the average offset by using the fact that the solar dipole vector seen by all the detectors is the same. 3. The third step is the determination of the 545 GHz frequency map offset, performed by minimizing its projected effect by the dust emission SED from 545 GHz to the lower frequencies offsets. The initial dust template is extracted from the 545 GHz frequency map, which is corrected at the end of this step. This provides a new dust template Dust ν 0 , used at the next iteration.
The residual R b also contains the not-recoverable noise and systematic effects residuals N b (including the very small residuals following the low frequency foregrounds removal).
In order to improve the solar dipole measurement, an improved model of the two main partially correlated high frequency foregrounds components is built extracting simultaneously the CO lines emission intensity coefficients γ b and the effective frequencies ν b and templates maps ∆Dust 1 and ∆Dust 2 describing the new thermal emission dust model. The constraint used is, for each sky pixel seen by all pairs of bolometers b1 and b2, to minimize the difference ∆R b1,b2 = R b1 − R b2 , expressed in µK CMB . The CMB tot term cancels and: (5) where D(ν b ) is the dust emission (Eq. 3). The initial dust map D(ν 0 ) is computed from the SRoll2.2 545 GHz map, after subtracting in an open loop CO, free-free and CMB. For bolometer b, the specific projection coefficient f (ν b )/ f (ν 0 ) from the template dust map computed from 545 GHz bandpass map, account for the detector response to the dust SED, through the effective frequency ν b . These effective frequencies, initially set to the values computed from the ground measurements of the bandpass (Planck Collaboration IX 2014), are adjusted simultaneously with the γ b coefficients.
The CO Planck Legacy maps are used as a template to extract a single effective set of coefficients the CO lines response γ b at 100, 217, 353 GHz. The fact that there is no CO line in the 143 GHz frequency band is used to regularize the solution, requiring that the detectors have a negligible response correlated with the CO template. Weaker interstellar molecular lines nevertheless exist in the 143 GHz frequency band, but this affects the CO coefficients by an error of 1.1%, inducing an error on the solar dipole amplitude of 0.06 µK and 10 −4 degree in direction, thus negligible.
The offsets O b are strongly degenerate with the dust large scale SED variation, thus are also solved for simultaneously using the same minimization of the differences ∆R b1,b2 used to compute the dust large scale SED variation. The quantity to minimize, based on the internal stability of the frequency maps, is, for each pixel p: Conditions on the mean parameter values between all single detector are needed to minimize Eq. 6 that is based on map differences. This is possible thanks to several closure conditions.
-For O b , ∆Dust 1 , and ∆Dust 2 , we impose the mean values to be null. -For γ b , we force the mean values to be null at 143 GHz.
-The level of residual dust emission in the mean map over all detectors after removing the dust and the CO model, and an initial CMB anisotropies map, is related to a miss estimation of the mean f (ν b ). Setting this level to 0, provides a closure condition for all ν b . At the first iteration, we use the CMB SMICA map.
From this minimization, we obtain: the correction to the response to a dust SED for bolometers calibrated on the CMB orbital dipole as an effective frequency ν b , the additive SED spatial variation described by the ∆Dust 1 and ∆Dust 2 maps, the CO effective coefficients γ b , the relative offsets O b of the single detector maps.
Furthermore, a global 100-353 GHz offsetŌ must be introduced to adjust the offsets of the 100-353 GHz maps with respect to the unknown 545 GHz map offset. Moreover, considering the large uncertainty on the 545 GHz map offset, a correction has also to be introduced. This will correct the large scales initial dust map template taken from the 545 GHz frequency map. This is performed by adding a constraint on the solar dipoles built for each detector (Sect. 3.3).

Determination of the frequency map offsets
The dust SED spatial variations are strongly degenerate with the average zero levels of the frequency maps. Equation 6 minimization only constrains the relative adjustments of these offsets, and do not affect the initial band average zero levels of the frequency maps which have rather large uncertainties with respect to the relative ones, especially for the 353 and 545 GHz frequency band.
A global (100 to 353 GHz) offsetŌ is introduced in the algorithm. It is determined thanks to the fact that the spatial distortions of the dust foreground initial dust map uncertainties, biases the solar dipole vector directions. Thus, the minimization of the dispersion of the single detector CMB solar dipole vector over the 43 (100 to 353 GHz) HFI bolometers allows to computeŌ.
To remove the solar dipole, the cosmological CMB anisotropies map must be removed from the CMB tot map. Nevertheless this CMB anisotropies map, extracted by component separation methods, is not reliable in the narrow galactic ridge that needs to be masked and filled with a constrained realization within a galactic ridge mask M g before setting the dipole term to zero on the whole sky. The mask M g is defined as the sky contained between two symmetrical latitudes centered on the galactic plane and characterized by the sky fraction f sky it covers.
The CMB anisotropies full sky map noted CMB M g is the CMB tot map from a component separation masked with M g in which the dipole term has been set to zero. It depends on the size of the galactic mask M g , and on the procedure used to fill this mask with CMB anisotropies constrained realizations. Before having a CMB tot anisotropies map consistent with the new foreground model, we use, at the first iteration, the Planck Legacy CMB SMICA map, estimated reliable for cosmology outside a 5% galactic mask. At the second and following iterations, a CMB large scales anisotropies map, coherent with Sect. 3, result can now be built. We remove from all single intensity bolometer maps R b from 100 to 353 GHz, the current dust model and CO template, and the adjustments of the offsets obtained in Sect. 3. The weighted average over all bolometers provides the new current CMB tot .
A new three-components (CMB, dust and CO) separation algorithm, named Bcsep, is developed using difference between single detector maps in nested iterative loops, also using the consistency of the solar dipole parameters between frequencies. This requires to start with good enough approximations and the convergence of the algorithm is the test of this condition. This new component separation is performed by subtracting the dust or the CMB maps from the frequency residual map R b and improving iteratively CMB tot and the three maps of the dust model per frequency: where σ 2 b is the variance of the cleaned maps for f sky =0.9 estimated by simulations.
Equation 8 defines the single detector intensity map R b as the residual map ℵ b,M g after removal of the high frequency foregrounds model, the current CMB anisotropies CMB M g , and the relative offsets extracted previously. ℵ b,M g is thus the map containing only the CMB solar dipole, together with the noise and residuals from the foregrounds removal: To constrainŌ, we impose that all detectors show a minimal dispersion of the solar dipole direction vector. Starting withŌ = 0, we loop on Eqs. 8, 9, and 10 and solve for theŌ value. The solar dipole parameters for bolometer b are fitted within a set of 10 masks M taken in the sky region left outside the galactic ridge M g to smaller and smaller areas identified by their f sky .
The minimization Eq. 9 extracts, for a given detector, a run-ningŌ and mask M g , the dipole amplitude A b,M g ,Ō and direction vector A b,M p ,Ō by minimizing the difference taken over for a mask M between the residual map ℵ b,p,M g map and the model solar dipole map plus the currentŌ: ThenŌ in carried back into Eq. 8 if we change M g or into Eq. 9 if we change only M after the optimization of M g , until convergence.
This should not be affected at the first iteration by the choice of the CMB cosmological anisotropies SMICA map as it converges to a stableŌ value.
At this stage, the 545 GHz map initial offset uncertainty still affects all lower frequencies zero levels indirectly through the distortion of the initial dust template and the minimization Eq. 6 to Eq. 10, and directly by the propagation of the 545 GHz offset residual to lower frequencies through the SED f (ν).
We define the 100 to 353 GHz offsets O ν as the average of (O b +Ō) over the frequency band ν. The correction of the initial offset at 545 GHz, O 545 , which also affects all the O ν (following the f (ν)/ f (ν 0 ) emission law), is the one which minimizes Eq. 11: where σ HI is the rms of the initial zero level determination uncertainties at frequency ν. Thus this correction will affect mainly the 353 GHz. The next section (Sect. 3.4) discusses the coherence of the offset determination.

Absolute zero level correction consistency
The absolute zero levels of the maps are critical for the knowledge of the foregrounds that need to be modeled with the right zero level. In Sects. 3.2 and 3.3, we develop a method to constrain these offsets using the solar dipole strong signal that should be seen in all single detector maps with the same direction vector.
Four map offsets are determined in sequence: zero level offset corrections, averaged per HFI frequency band, are given as cumulative offset values obtained after each correction of the algorithm refining the initial offset O HI . The convergence of the iterative process shows the stability of the algorithm which leads to a minimal dispersion of the solar dipole vector direction for the single detector CMB anisotropies maps in the frequency bands 100-353 GHz. The final absolute zero level uncertainties are computed using a set of ten simulations for the four lower frequencies. Those simulations have representative noise and systematic effects but no attempt is done to quantify the uncertainties due to the dust model by simulations which have no reliable statistical model. These uncertainties at 100, 143, and 217 GHz are reflected in the dispersion of the single bolometers at each frequency. Thus, these errors are probably dominant over the instrumental systematic effect dispersion. The uncertainty at 353 GHz is about one third of the dispersion between bolometers which indicate that the zero level correction reflects the dust foreground effects. The uncertainty at 545 GHz is computed using a Monte Carlo approach based on uncertainties at the lower frequencies.
The O 545 correction induces a marginally significant correction at 353 GHz only. These corrections steps are weakly correlated, and the convergence confirms that the algorithm steps are efficient to improve the offsets. Figure 1 displays the map offset results, the initial ones from Planck in red, and the final ones as a function of the frequency bands. The offset corrections from this work, in blue, are consis- tent with the initial frequency map offsets but with much better accuracies by factors 2.5 to 5.

Dust model results
Figure 2 presents these additive correction maps (as Healpix maps at N side =128) to take into account the spatial variation of the dust emission SED. These spatial variation maps, in the right column, are built from ∆Dust 1 and ∆Dust 2 maps, respectively weighted by the frequency dependent ratios ∂ f (ν)/∂β and ∂ 2 f (ν)/∂β 2 . Those maps are compared with the HFI-Legacy ones (left column), where each frequency SED spatial variation map has been determined independently. The large scales are similar although obtained in different ways, and despite the fact that the initial dust maps were taken from different frequencies: 857 GHz for HFI-Legacy and 545 GHz for this work.
The SED spatial variations are described for by an expansion in ln(ν b /ν 0 ) at the second order, but associated with two independent large scales maps providing more degrees of freedom. The two bottom maps in Fig. 3 shows these additive corrections maps for the 143 GHz frequency band. The top map shows an estimate of an effective β corresponding to the total correction of the two terms.
To be more quantitative, Fig. 4 compares the ratio of the power spectra of the first and second order additive SED correction maps at 143 GHz, with different f sky . There is a strong galactic center anti-center dipole term ( = 1 and 3), and weaker latitude pattern ( = 2 and 4), and then a white power distribution, probably reflecting the random interstellar cloud distribution. The low rise at > 15 shows the noise contribution. For f sky =0.82, the first order SED variation corrections show a strong centre to anticentre dipole term at the 30 to 50% level of the dust initial spectra with opposite sign to the dust map. The second order correction is two orders of magnitude smaller (one order of magnitude in the map). This agrees with the astrophysical assumption that the dust SED gradients are induced by large scale galactic physics, -gradient of gas molecular fraction in latitude and star light energy density gradient with galactic radius -, and not small scales variations of the dust properties which are seen in the flat power spectrum (5 < < 15).
To test that the dust model is valid up above 545 GHz, we describe, in Appendix A, its extension to 857 GHz.

CMB large scales results and solar dipole determination
The solar dipole is extracted from each single detector CMB tot map from which we can subtract CMB anisotropies. The first iteration uses the Planck Legacy CMB anisotropies SMICA map, after masking the galactic ridge keeping 95% of the sky as recommended in Planck Collaboration IV (2020). In later iterations, the Bcsep CMB tot map from Equation 8 is used.

New CMB large scale anisotropies map
To illustrate the accuracy of the CMB cosmological anisotropies map built by the Bcsep method, we use an end-to-end simulation which contains CMB, noise, systematic effects residuals, and also a CO model and a dust model as per Eq. 2. Thus, this simulation does not intend to evaluate the Bcsep ability to solve for the dust model, but only for the noise and systematic effects residuals. Figure 5 shows the difference between the simulated input CMB primary anisotropies map (without solar dipole) and the retrieved one. The narrow galactic ridge towards the inner galaxy is the only visible contributor to a significant dipole term, with a main vector direction toward b = l = 0 • . To avoid this bias, the solar dipole is extracted outside the galactic ridge mask M g . The error induced on the solar dipole amplitude is 0.103 µK with f sky = 0.90, and 0.108 µK with f sky = 0.95. This shows that the noise and systematic effect on the CMB maps removal outside M g , within this range of f sky , is not a direct problem for the solar dipole determination. The dust foreground residuals remains the main issue. Figure 6 displays the difference maps between the CMB tot maps obtained by Bcsep and the four Planck and Bcsep at b > 20 • does not show significant residuals either. All differences containing either the Commander or the SEVEM maps show zodiacal dust emission residuals along the ecliptic equator. We note that the CMB tot Bcsep map is only aimed at a significant improvement at < 30 and outside the bright galactic ridge M g where it has smaller large scale dust residuals than the Planck legacy ones. Bcsep does not intend to build a CMB map aimed at the analysis of cosmological parameters mostly based on CMB power spectra at > 30, and is not characterized for that purpose.

Optimization of the CMB anisotropies removal
Lacking a model to build statistically representative and controlled dust foreground maps, the foreground residuals cannot be fully evaluated by simulations. A semi-quantitative approach is nevertheless possible by comparing data and simulations. The behavior of the solar dipole parameters with frequency and sky fraction are sensitive probes to test both the CMB anisotropies removal method and the galactic foreground residuals.
Large scale foreground residuals in the part of the sky where the dipole is fitted, induce a source of systematic effects on the dipole parameter retrieval. These effects are introduced in the simulation by adding an empirical estimate of the foreground residuals. This residual is defined as the difference map from two component separation method : SMICA-Bcsep shown in Fig. 6.
We now use the CMB tot Bcsep map built in Sect. 4.1. When extracting the solar dipole from the full sky CMB tot , the solar dipole parameters are biased by the galactic ridge residuals seen in Fig. 6. This leads to the use of the galactic ridge mask M g where the CMB tot is affected by large galactic residuals. This mask M g is filled with a simulated CMB constrained at the boundaries (Wandelt et al. 2004;Thommesen et al. 2020). By definition, the all sky CMB cosmological anisotropies map has no dipolar term, and shall be built by removing the dipolar term from the CMB tot map. The solar dipole is thus fitted in a mask M taken as the complement of M g . On the one hand, the masked sky replaced by a constrained CMB, induces an uncertainty on the dipole parameters. This uncertainty is evaluated by the dispersion of the retreived solar dipole parameters from 1000 Monte Carlo realizations with various M g . On the other hand, the bias, induced by the dust residuals in the central galactic ridge, decreases when increasing M g . The combination of these two effects presents a minimum, which defines the optimal. The large scale dust residuals in the CMB tot map also affect slightly the dipole measurement (Sect. 4.1). Any detectable dust residual effect will then appear as a frequency and/or an f sky dependence of the dipole parameters when changing M g .
We want to compare these two effects: the large scale SED spatial variations of the dust emission, and the CMB tot filling of M g . These behaviors, observed with the data, are simulated. Figure 7 shows the computed solar dipole parameters for the 100-353 GHz frequencies, when varying the galactic ridge mask M g in which the CMB tot Bcsep map is filled with constrained CMB realizations. The dipole parameters are plotted as a function of the sky fraction left outside M g used to fit the dipole.
At first, the systematic effects affecting the retrieved dipole parameters are tested both on simulations and data to ensure that the simulations are a good representation of the data. The two first rows if Fig. 7 use the processing algorithm without SED variation correction. The input solar dipole parameters for the simulated data are taken from the final determination from this work, and given by the black lines with grey bands. In these two first rows, we see the remaining effects of the dust SED variations effects (included in the simulations and not corrected for) well simulated when compared to the results obtained on the data. The magnitude and sign of the residuals are qualitatively demonstrating the quality of the simulation even though the dust behavior is not quantitatively exact as expected, lacking a reliable physical model of the dust SED variations. In the third and fourth rows of Fig. 7, the algorithm now includes the dust SED spatial variation correction. Both show a spectacular reduction of the previous trends. This demonstrates that the simulation takes into account the SED variations with the right order of magnitude, and that also the algorithm captures them well.
The dipole term has to be removed from the CMB tot map, well known only outside a galactic ridge mask M g where there are significant galactic dust residuals. This masking induces an error due to leakages between the very large scale anisotropies and the dipole term. Figure 8 shows the results on the recovered parameters of the solar dipole as a function of the sky fraction left outside of the varying M g where the dipole is fitted. The dispersion of the solar dipole parameters is evaluated by simulation with 1000 realizations of the CMB filling of M g . The plain gray wedges show the dipole term dispersion when nothing is done to compensate this leakage. The dashed (resp. doted) gray wedges show the dispersion of the dipole term after filling M g with a smoothed 1 • (respectively 6 • ) CMB realization of the large scales only (up to = 5) constrained by the cosmological parameters HFI best-fit (Planck Collaboration V 2020), and continuity constraints at the boundaries. Once the dipole term is removed from the CMB tot map, as per the adopted definition, it is used to remove the CMB anisotropies from the full sky map. Then, the solar dipole is fitted in M. This removal and filling procedure introduces a non-recoverable uncertainty, increasing with M g estimated using results shown in Fig. 8.
The simulations with the qualitatively representative dust component residuals defined above, show a bias as expected in the amplitudes and in direction (longitude bias as expected larger than the latitude one) of the solar dipole parameters. The biases are maximal when the galactic ridge mask M g is null, and decreases with decreasing f sky following the predicted behavior. For no masking of the galactic ridge, the direction is shifted by about 2 of longitude. When f sky decreases to 0.95, the longitude bias is reduced by a factor of 2, and crosses the longitude dispersion for 1 • smoothing. For this sky fraction, the amplitude parameter shows a bias of 0.25 µK. The dispersion for 6 • smoothing crosses the bias for f sky =0.92. We thus conclude that, for the CMB anisotropies removal map, a M g mask leaving f sky =0.95 available for the fit of the dipole is optimal for the CMB anisotropies removal with a safer 1 • smoothing. Figure 9 is a zoom-in of Fig.7 third row. The dipole parameters for all frequencies show the behavior described in Sect. 4.2 with a bias increasing for increasing from f sky =0.8 to 1.0 (purple line), and a dispersion, as the gray wedge (1 • smoothing), decreasing rapidly with increasing f sky . The four frequencies show very similar patterns. The longitude, the most sensitive parameter, shows a minimal error for the optimal M g size to f sky =0.95 as expected. The behavior of the three parameters shows a correlated, but non monotonic, behavior as expected when dominated by the dispersion induced by the large scales foregrounds for f sky < 0.95, within the range predicted by the gray wedge. It is at present impossible to draw high quality statistically representative galactic foreground maps, and thus fully evaluate the foreground residuals. The uncertainties are given by the dispersion between detectors in each frequency band, and increase with fre-quency. Calibration errors resulting from noise and systematic effects, both evaluated by simulations, are visible as a variation in amplitude between frequencies.

Solar dipole results
We want now to test the effect of the dust residuals on the stability of the solar dipole parameters, when increasing frequencies, and varying M keeping the optimal M g constant ( f sky =0.95). Figure 10 shows the stability of the solar dipole parameters when the fitting mask M spanning from f sky =0.19 to 0.9 (see Sect. 4.2). The solar dipole parameters stability is impressive, and do not show any obvious trend with the dipole fitting mask M nor with frequency. The most noticeable discrepancies are in amplitude and nearly all within 1σ. The final results of this paper (average over the three lower frequencies) is shown as black lines and grey band uncertainties. We thus conclude that there is no dust foreground residuals detectable above the noise. Table 2 shows the final results of the solar dipole parameters, computed as the weighted averaged per frequency over M with f sky =0.3 to 0.8. For amplitudes, the dispersion between detectors within a frequency band is a more relevant measure of the whole calibration process, including systematic effects.HFI-Legacy uses 100 full end-to-end simulations of the overall absolute calibration process based on the recovery of the solar dipole input. These simulations retrieve the rms uncertainties of the cal- ibration measured by the solar dipole amplitude rms: 1.5 × 10 −4 at frequency from 100 to 217 GHz increasing to 4 × 10 −4 at 353 GHz. Our results in amplitude are coherent with the HFI-Legacy calibration mismatch estimates. Thus, the weighted average (AVG) of the 100, 143, and 217 GHz results provides the best dipole parameters estimate. The AVG amplitude uncertainty (0.36 µK), computed as the dispersion between the three frequencies, is significantly larger than the intra-frequency uncertainties, but reflects the calibration error due the residual systematic effects.
For direction, the uncertainties are the statistical errors. The HFI-Legacy results are obtained by alignment of the 143-353 GHz frequencies solar dipole directions on the 100 GHz solar dipole direction, assuming it was the best reference, because of its best combination of high CMB sensitivity and low fore- Galactic ridge residual CMB large multipole leakages   grounds. In the present work, this assumption is replaced by the constraint of minimum dispersion of the dipole directions for all 100 to 353 GHz detectors, a more refined dust model of the SED spatial variation, and coherent large scales CMB anisotropies. We also introduce a better masking-filling procedure for the removal of the galactic ridge residuals in the CMB cosmological anisotropies. Figure 10 compares the present results with the HFI-Legacy results, shown as the pink lines and bands. The results for the amplitude and latitude are well compatible within their error bars. Nevertheless, the longitude is not, although very stable when varying f sky , and consistent with SMICA, and NILC.
An empirical test of the systematic effects related to the component separation is to use different CMB maps. This is done also in Fig. 10, and the four component separation method used for the CMB anisotropies are shown with thin lines. Although using CMB from Commander, NILC, SEVEM, and SMICA is not fully consistent with the new and better foreground dust model developed in this work, the comparison remains interesting. When other component separations are used to get the CMB anisotropies, the stability with varying f sky (M) parameter is not as good as the one achieved by the Bcsep one. It strengthens the statement that the CMB removal, when done with an optimal galactic ridge mask, has little effect on the solar dipole parameters fitted in a mask M covering a very broad range of sky fraction, thus showing no sign of large scales dust residuals. Thommesen et al. (2020) andPlanck Collaboration et al. (2020) increasing the galactic mask up to 80% of the sky, is far from optimal. The dispersion of the amplitude in these papers grows from ±2.5 µK for M g in the range f sky 0.95 to 0.75 and reach ±10 µK for f sky =0.30. These large uncertainties are in-duced when setting the dipole term to zero in a larger and larger galactic mask. Conversely, the method described in the present work, optimizes this mask and reduce by large factors the uncertainties making full use of the CMB anisotropies accurately measured up to 95% of the sky. We also note that the 4σ calibration discrepancy between 100 and 143 GHz mentioned in Table 10 if Planck Collaboration et al. (2020) is of the same order as the absolute calibration reflected by the shift of the amplitude of the solar dipole by 4 µK in their paper showing a photometric calibration problem. Figure 11 shows three combinations of the solar dipole parameters, results of this work, and of the HFI-Legacy results (in red  Fig. 11. Three combinations of the solar dipole parameters as a cross together with uncertainty ellipses. Red color shows the HFI-Legacy results. The blue color is the present work but without filling the galactic ridge mask M g , dotted when ignoring the systematic effect on the CMB anisotropies dipolar term, full line when taking into account the systematic effect. The black color is for the results of the present work with galactic ridge M g filled with a constrained CMB. shifted by −3.6 (2.5 on the sky) well outside the error ellipses. There is a key difference between the extraction method used in this work and the one used in HFI-Legacy. The latter fit the dipole parameters outside of a galactic ridge mask for which the zero dipolar term of the CMB anisotropies has been removed without filling the mask with constrained CMB realizations non contaminated by dust residuals. This induces a bias on the dipolar component in the part of the sky outside the galactic ridge mask M g . This bias has been simulated on the present final maps, by removing the step of the gap filling, and fitting the dipole parameters outside the optimal mask M g . The dipole direction is displayed as the blue cross and dotted error ellipse in the top panel of Fig. 11. This point falls near the HFI-Legacy red point, as expected as both of them have not the optimal gap filling. The bias on the sky is 1.8 when the shift observed is 2.5 . This simulation shows the bias induced when M g is not filled with a constrained CMB, and we added this bias to the uncertainties of the legacy shown by the full blue line ellipse.
Averaging the three lowest Planck-HFI CMB frequencies (Table 2) with uncertainties including the systematic effects leads thus to the final solar dipole parameters: This result is consistent with HFI-Legacy for amplitude and latitude. The longitude of the dipole axis is shifted by −2.5 on the sky. This is fully explained by the difference of procedure, having introduced the filling of the galactic ridge mask with constrained CMB anisotropies realizations before setting its dipolar term to zero.

Conclusion
This work provides the best measurement of the Solar dipole parameters as reminded in Table 3, showing an unprecedented stability. This measurement, together with the large scale CMB cosmological anisotropies map can be used for intercalibration of future CMB experiments, and testing the removal of foregrounds as done in this work. The dust emission is the dominant large scale anisotropies component above 100 GHz affecting strongly the CMB at the largest angular scales. Such non-Gaussian foreground critically depends on the zero level of the frequency maps which can only be constrained for differential experiments by the coherence of the dipoles. We have developed here a method which does improve these zero levels which have much better accuracy than the extrapolation to zero column density of interstellar gas. Moreover, the component separation between the CMB and the dominant dust foreground is performed iteratively and allows to build intermediate and large scales maps of the dust SED variations.
The data products associated to this paper are thus: improved zero levels of the HFI single detector maps, refined CMB solar dipole vector parameters, a new CMB anisotropies map at large scales, -100 to 857 GHz dust intensity maps of our new dust model including the SED spatial variation properties Figure 12 illustrates this dust model, which can be used to better understand the polarized dust emission, that should be affected by similar spatial variation.  The Planck HFI data will stay, for at least a decade, the best all sky data at the high frequencies. At frequencies higher than 200 GHz, the sky signal is strongly affected by the atmosphere. The present study products are thus unique and useful to combine with balloon borne and ground based CMB data waiting for the next CMB space mission. present dust model (based on the 545 GHz). This suggests that the dust model can be extrapolated to 857 GHz.
The 857 GHz map is the sum of the dust emission, an unknown offset correction, noise and systematics. The CMB anisotropies and low frequency foregrounds are fully negligible here. Thus, the 857 GHz map is used as a dust template which is propagated at 545 GHz using Eq. 3: Dust ν 857 − ln ν 857 ν 545 ∆Dust 1 − ln 2 ν 857 ν 545 ∆Dust 2 A 545 GHz dust model D (ν 545 ) is built from the 857 GHz map by propagating the SED variation dust model. This model is removed from the 545 GHz single detector map, only leaving a CMB solar dipole map plus noise and systematics residuals. To constrain the 857 GHz map offset O 857 , we compute, for each single detector 545 GHz map, the solar dipole parameters following the method described in Sect. 4. The 857 GHz offset is adjusted by minimizing the difference between the 545 GHz solar dipole directions and the 100-217 GHz average (Table 2), and found to be O 857 =0.124 MJy sr −1 . Figure A.1 shows the CMB dipole parameters extracted from the three single detector 545 GHz maps for both a null 857 GHz offset and for the solved one. The convergence and stability (amplitude dispersion about 50 µK) at such high frequencies demonstrate the quality of the dust model and the fact that the deviations do not come from the dust removal (like in the HFI-Legacy) but more from the photometric calibration errors associated with instrument systematic effects. Figure A.2 shows the consistency of the two dust maps at 545 GHz (top row). The differences are shown in the bottom row for the full resolution, and at 5 • smoothing. Even if there are some differences near the galactic ridge, both maps are very consistent at high latitudes, demonstrating that the dust model stands at large scale and high latitude up to 857 GHz.