L Mousset
email: mousset@apc.in2p3.fr

M M Gamboa Lerena

E S Battistelli

P De Bernardis

P Chanial

G D'alessandro

G Dashyan

M De Petris

L Grandsire

J.-Ch Hamilton

F Incardona

S Landau

S Marnieros

S Masi

A Mennella

C O'

Sullivan M Piat

G Ricciardi

C G Scóccola

M Stolpovskiy

A Tartari

J.-P Thermeau

S A Torchinsky

F Voisin

M Zannoni

P Ade

J G Alberro

A Almela

G Amico

L H Arnaldi

D Auguste

J Aumont

S Azzoni

S Banfi

B Bélier

A Baù

D Bennett

L Bergé

J.-Ph Bernard

M Bersanelli

M.-A Bigot-Sazy

J Bonaparte

J Bonis

E Bunn

D Burke

D Buzi

F Cavaliere

C Chapron

R Charlassier

A C Cobos Cerutti

F Columbro

A Coppolecchia

G De Gasperis

M De Leo

S Dheilly

C Duca

L Dumoulin

A Etchegoyen

A Fasciszewski

L P Ferreyro

D Fracchia

C Franceschet

K M Ganga

B García

M E García Redondo

M Gaspard

D Gayer

M Gervasi

M Giard

V Gilles

Y Giraud-Heraud

M Gómez Berisso

M González

M Gradziel

M R Hampel

D Harari

S Henrot-Versillé

E Jules

J Kaplan

C Kristukat

L Lamagna

S Loucatos

T Louis

B Maffei

W Marty

A Mattei

A May

M Mcculloch

L Mele

D Melo

L Montier

L M Mundo

J A Murphy

J D Murphy

F Nati

E Olivieri

C Oriol

A Paiella

F Pajot

A Passerini

H Pastoriza

A Pelosi

C Perbost

M Perciballi

F Pezzotta

F Piacentini

L Piccirillo

G Pisano

M Platino

G Polenta

D Prêle

R Puddu

D Rambaud

E Rasztocky

P Ringegni

G E Romero

J M Salum

A Schillaci

S Scully

S Spinelli

G Stankowiak

A D Supanitsky

P Timbie

M Tomasi

G Tucker

C Tucker

D Viganò

N Vittorio

F Wicek

M Wright

A Zullo

QUBIC II: Spectro-Polarimetry with Bolometric Interferometry

Keywords:

Bolometric Interferometry is a novel technique that has the ability to perform spectroimaging. A Bolometric Interferometer observes the sky in a wide frequency band and can reconstruct sky maps in several sub-bands within the physical band. This provides a powerful spectral method to discriminate between the Cosmic Microwave Background (CMB) and astrophysical foregrounds. In this paper, the methodology is illustrated with examples based on the Q & U Bolometric Interferometer for Cosmology (QUBIC) which is a ground-based instrument designed to measure the B-mode polarization of the sky at millimeter wavelengths. We consider the specific cases of point source reconstruction and Galactic dust mapping and we characterize the Point Spread Function as a function of frequency. We study the noise properties of spectro-imaging, especially the correlations between sub-bands, using end-toend simulations together with a fast noise simulator. We conclude showing that spectroimaging performance are nearly optimal up to five sub-bands in the case of QUBIC.

1 Introduction

This article is the second in a series of eight on the Q & U Bolometric Interferometer for Cosmology (QUBIC), and it introduces the spectroscopic imaging capability made possible by bolometric interferometry. QUBIC will observe the sky at millimeter wavelengths, looking at the Cosmic Microwave Background (CMB). The CMB, detected by Penzias and Wilson [START_REF] Penzias | A measurement of excess antenna temperature at 4080 Mc/s[END_REF], has a thermal distribution with a temperature of 2.7255 ± 0.0006 K [START_REF] Fixsen | The temperature of the cosmic microwave background[END_REF]. It is a partially polarized photon field released 380,000 years after the Big Bang when neutral hydrogen was formed. The polarization can be fully described by a scalar and a pseudo-scalar fields: E and B-modes respectively. B-mode polarization anisotropies are generated by primordial gravitational waves occurring at the inflation era. The indirect detection of these waves would represent a major step towards understanding the inflationary epoch that is believed to have occurred in the early Universe. Tensor modes in the metric perturbations are a specific prediction of Inflation. The measurement of the corresponding B-mode polarization anisotropies would reveal the inflationary energy scale, which is directly related to the amplitude of this signal. This amplitude, relative to the scalar mode, is parametrized by the so called tensor-to-scalar ratio r.

Currently, there are several instruments aiming at measuring the primordial B modes. These include SPTPol [START_REF] Sayre | Measurements of B-mode Polarization of the Cosmic Microwave Background from 500 Square Degrees of SPTpol Data[END_REF], POLARBEAR [START_REF] Ade | A measurement of the cosmic microwave background B-mode polarization power spectrum at sub-degree scales with POLARBEAR[END_REF], ACTPol [START_REF] Louis | The Atacama Cosmology Telescope: two-season ACTPol spectra and parameters[END_REF], and BICEP2 [START_REF] Ade | Constraints on Primordial Gravitational Waves Using Planck, WMAP, and New BICEP2/Keck Observations through the 2015 Season[END_REF]. Planned experiments include CLASS [START_REF] Dahal | The class 150/220 ghz polarimeter array: Design, assembly, and characterization[END_REF], POLARBEAR 2 + Simons Array [START_REF] Suzuki | The Polarbear-2 and the Simons Array Experiments[END_REF], Simons Observatory [START_REF] Ade | The Simons Observatory: science goals and forecasts[END_REF], advanced ACT [START_REF] Henderson | Advanced ACTPol Cryogenic Detector Arrays and Readout[END_REF], PIPER [START_REF] Gandilo | The primordial inflation polarization explorer (piper)[END_REF], upgrade of the BICEP3/Keck array [START_REF] Grayson | Bicep3 performance overview and planned keck array upgrade[END_REF], LSPE [START_REF] Addamo | The large scale polarization explorer (LSPE) for CMB measurements: performance forecast[END_REF], CMB-S4 [START_REF] Abazajian | CMB-S4 Science Case, Reference Design, and Project Plan[END_REF] and LiteBird [START_REF] Suzuki | The litebird satellite mission: Sub-kelvin instrument[END_REF].

The claim in 2014 of the discovery of primordial B-modes [START_REF] Ade | Detection of B-mode polarization at degree angular scales by BICEP2[END_REF] was contradicted by the observation of a significant level of dust contamination in the BICEP2 field by the Planck satellite [START_REF] Ade | Joint analysis of BICEP2/Keck Array and Planck data[END_REF] using the 353 GHz polarized channel from the High Frequency Instrument. This unfortunate event showed in a very clear manner how important is the control for contamination by polarized foregrounds in the search for B-mode polarization. The main foreground contaminant at high frequencies is the polarized thermal emission from elongated dust grains in the Galaxy [START_REF] Akrami | Planck 2018 results-XI. Polarized dust foregrounds[END_REF]. At lower frequencies, emission from synchrotron [START_REF] Fuskeland | Spatial variations in the spectral index of polarized synchrotron emission in the 9 yr WMAP sky maps[END_REF] is expected to be significant, even in a so-called clean CMB field, if r is below 10 -2 . Current estimates depend strongly on the assumption that synchrotron is well described by a simple power law with a steep spectral index (but spectral curvature might well be there). While free-free is not expected to be a major contaminant due to its vanishingly small degree of polarization [START_REF] Macellari | Galactic foreground contributions to the 5-year Wilkinson Microwave Anisotropy Probe maps[END_REF], spinning dust, whose impact is addressed in [START_REF] Remazeilles | Sensitivity and foreground modelling for large-scale cosmic microwave background B-mode polarization satellite missions[END_REF], should not be neglected.

The control of contamination from foregrounds can only be achieved with a number of frequencies around the maximum emission of the CMB relying on the fact that the spectral distribution of the CMB polarization is significantly different from that of the foregrounds so that low and high frequencies can be used as templates to remove the contamination in the CMB frequencies.

A bolometric interferometer such as QUBIC has the ability to perform spectro-imaging, thanks to its very particular synthesized beam. The instrument beam pattern, given by the geometric distribution of an array of apertures operating as pupils of the interferometer, contains multiple peaks whose angular separation is linearly dependent on the wavelength. As a result, and after a non-trivial map-making process, a bolometric interferometer such as QUBIC can simultaneously produce sky maps at multiple frequency sub-bands with data acquired in a focal plane containing bolometers operating over a single wide frequency band.

This article is organized as follows. In section 2 we describe the working principle of a bolometric interferometer taking the example and characteristics of QUBIC. In section 3 we describe the spectral dependence of the synthesized beam and show under which conditions the spectral information can be recovered at the map making level. In section 4 we test spectro-imaging on simple cases using the QUBIC data analysis and simulation pipeline. Finally, in section 5 we present the performance of the spectral reconstruction. We compare a simulated sky to the reconstructed one using the QUBIC data analysis pipeline. Tests with a real point source were carried out in the lab and results are presented in [START_REF] Torchinsky | QUBIC III: Laboratory Characterization[END_REF].

Detailed information about QUBIC can be found in the companion papers: Scientific overview and expected performance of QUBIC [START_REF] Hamilton | QUBIC I: Overview and Science Program[END_REF], Characterization of the Technological Demonstrator (TD) [START_REF] Torchinsky | QUBIC III: Laboratory Characterization[END_REF], Transition-Edge Sensors and readout characterization [START_REF] Piat | QUBIC IV: Performance of TES Bolometers and Readout Electronics[END_REF], Cryogenic system performance [START_REF] Masi | QUBIC V: Cryogenic system design and performance[END_REF], Half Wave Plate rotator design and performance [START_REF]QUBIC VI: cryogenic half wave plate rotator, design and performances[END_REF], Feedhorn-switches system of the TD [START_REF] Cavaliere | QUBIC VII: The feedhorn-switch system of the technological demonstrator[END_REF], and Optical design and performance [START_REF] O'sullivan | QUBIC VIII: Optical design and performance[END_REF].

Bolometric Interferometry as Synthesized imaging

A bolometric interferometer is an instrument observing in the millimeter and submillimeter frequency range based on the Fizeau interferometer [START_REF] Fizeau | Sur les hypothèses relatives à l'éther lumineux et sur une expérience qui parait démontrer que le mouvement de corps change la vitesse avec laquelle la lumière se propage dans leur intérieur[END_REF]. A set of pupils are used at the front of the instrument to select baselines. The resulting interference pattern is then imaged on a focal plane populated with an array of bolometers. With the addition of a polarizing grid and a rotating Half-Wave-Plate before the pupil array, the instrument becomes a polarimeter, such as QUBIC.

Each pair of pupils, called a baseline, contributes with an interference fringe in the focal plane. The whole set of pupils produces a complex interference pattern that we call the synthesized image (or dirty image) of the source. For a multiplying interferometer, the observables are the visibilities associated with each baseline. In bolometric interferometry, the observable is the synthesized image, which can be seen as the image of the Inverse Fourier Transform of the visibilities.

The left panel of figure 1 shows the array of back-to-back feedhorns of the Full Instrument (FI). The back-to-back feedhorn array makes 400 pupils arranged on a rectangular grid within a circle (see [START_REF] Cavaliere | QUBIC VII: The feedhorn-switch system of the technological demonstrator[END_REF][START_REF] O'sullivan | QUBIC VIII: Optical design and performance[END_REF] for more details). The beam looking at the sky is called the primary beam, and secondary beam is the one looking toward the focal plane. In the right panel, we can see the synthesized image obtained on the focal plane, composed by 992 bolometers, when the instrument is looking at a point source located on the optical axis in the far field.

Synthesized imaging

A bolometer is a total power detector. The signal S η (r, λ) on a point r = (x, y) of the focal plane1 is the square modulus of the electric field E η (t, n, λ) averaged over time t and integrated over all sky directions n. The signal with polarization η from each direction is reemitted by each of the pupils resulting in a path difference in the optical combiner. The signal on the focal plane depends on the location of each pupil h j , its primary beam B prim (n, λ), the focal length of the combiner f , the secondary beam of the pupil on the focal plane B sec (r, λ), and the wavelength λ:

S η (r, λ) = B prim (n, λ)B sec (r, λ) j E η (t, n, λ) × exp i2π h j λ • r f 2 + r 2 -n 2 dn, (2.1)
where r is the norm of r. We use the concept of the Point Spread Function (PSF) of the synthesized beam,

P SF (n, r, λ) = B prim (n, λ)B sec (r, λ) × j exp i2π h j λ • r f 2 + r 2 -n 2 , (2.2)
to rewrite eq. 2.1 as

S η (r, λ) = |E η (t, n, λ)| 2 × P SF (n, r, λ) dn. (2.
3)

The rotating Half Wave Plate (HWP) modulates the polarized signal, with a varying angle φ, so we can write eq. 2.3 in terms of the synthesized images of the Stokes parameters:

S(r, λ) = S I (r, λ) + cos(4φ)S Q (r, λ) + sin(4φ)S U (r, λ) (2.4)
where the synthesized images are the convolution of the sky through the synthesized beam:

S X (r, λ) = X(n, λ) × P SF (n, r, λ) dn, (2.5)
X standing for the Stokes parameters I, Q or U. The signal received in the detectors with a bolometric interferometer is therefore exactly similar to that of a standard imager: the sky convolved with a beam. The only difference being that this beam is not that of the primary aperture system (telescope in the case of an ordinary imager) but is given by the geometry of the input pupil array and the beam of the pupils (see eq. 2.2). With such an instrument, one can scan the sky in the usual manner with the synthesized beam, gathering Time-Ordered-Data (TOD) for each sky direction (and orientation of the instrument) and reproject this data onto a map at the data analysis stage (see section 2.4).

Realistic case

Note that in a real detector the signal is integrated over the wavelength range defined by filters and also over the surface of the detectors [START_REF] Battistelli | The QU bolometric interferometer for cosmology[END_REF]. If one assumes that the sky signal does not vary within the wavelength range, the expression for the signal, eq. 2.4, is unchanged and one just needs to redefine the synthesized beam as

P SF (n, r d , λ k ) = P SF (n, r, λ)J λ k (λ)Θ(r -r d ) dλdr (2.6)
where J λ k (λ) is the shape of the filter for the band centered in λ k and Θ(r) represents the top-hat function for integrating over the detector whose center is at r d .

B prim × B sec (r = 0) B prim × B sec (r = 12mm) r = 0 r = 12 mm Figure 2.
Cut of the synthesized beam as a function of θ (angle between n and the optical axis) given by eq. 2.7 for a square array of 20 × 20 pupils separated by ∆h = 14 mm at 150 GHz frequency (2 mm wavelength) for a detector located at the center of the optical axis in blue and 12 mm apart in cyan. Dashed lines represent the primary beam of the pupils (here Gaussian). Resolution and peak separation depend linearly on the wavelength λ.

The monochromatic synthesized beam

If the pupil array is a regular square grid of P pupils on a side spaced by a distance ∆h, the sum in eq. 2.2 can be analytically calculated (see [START_REF] O'sullivan | QUBIC VIII: Optical design and performance[END_REF] for more detail) and the monochromatic point-like synthesized beam, assuming f is large enough to use the small-angle approximation, becomes

P SF (n, r, λ) = B prim (n, λ)B sec (r, λ)× sin 2 P π ∆h λ x f -n x sin 2 π ∆h λ x f -n x sin 2 P π ∆h λ y f -n y sin 2 π ∆h λ y f -n y . (2.7)
where n = (n x , n y ) is the off-axis angle of the source. In such a case, the synthesized beam of the monochromatic point-like detectors has the shape of a series of large peaks with ripples in between, modulated by the primary beam of the pupils. For illustrative purpose, a cut of the synthesized beam is shown in figure 2 for such a square 20x20 array of horns. The realistic synthesized beam corresponding to our circular array (see figure 1) is shown in [START_REF] O'sullivan | QUBIC VIII: Optical design and performance[END_REF] (figure 11) and shows minor differences. Figure 2 shows this approximate synthesized beam for two different detectors emphasizing the fact that the location of the peaks moves with the detector location in the focal plane. Moreover, the intensity received by the detector changes inversely with r so the two Gaussian envelopes also differ. From the expression in eq. 2.7 it is straightforward to see that the FWHM of the large peaks is roughly given by FWHM = λ (P -1)∆h while their separation is θ = λ ∆h as illustrated in figure 2. The real synthesized beam of the QUBIC Technological Demonstrator has been measured in the lab using a monochromatic calibration source and is presented in [START_REF] Torchinsky | QUBIC III: Laboratory Characterization[END_REF]. For the Technological Demonstrator, the horn array is an 8 × 8 square array. The measured synthesized beam is in overall good agreement with the approximate analytical expression in eq. 2.7 and figure 2 although excursions from this perfect case are expected due to optical defects and diffraction in the optical chain. We have shown in [START_REF] Bigot-Sazy | Self-calibration: an efficient method to control systematic effects in bolometric interferometry[END_REF] that the shape of the synthesized beam for each detector can be precisely recovered through the "self-calibration" technique that is heavily inspired from synthesis imaging techniques [START_REF] Cornwell | A new method for making maps with unstable radio interferometers[END_REF].

Monochromatic map-making

Before discussing spectro-imaging, we first describe the map-making with a Bolometric Interferometer in the monochromatic case. As shown in section 2.1, the instrument is essentially equivalent to a standard imager, scanning the sky with the synthesized beam, producing Time-Ordered-Data (TOD) that can be projected onto sky maps. The map-making will therefore be very similar to that of a standard imager. In the monochromatic case, if the sky signal is s, the TOD y can be written as:

y = H • s + n (2.8)
where n is the noise and H is an operator that describes the convolution by the synthesized beam in eq. 2.5 as well as the pointing at the different directions of the sky. H is a 2dimensional matrix: number of time samples (scaling with the number of detectors) × number of sky pixels. The noise has two contributions: photon noise and detector noise. Photon noise is the Poisson fluctuations from the temperature of the CMB (T CM B 2.7K), the atmosphere, and the internal optical components. Detector noise is given by the Noise Equivalent Power (NEP) measured in each detector. Eq. 2.8 will be generalized to the case of several frequency sub-bands in section 3.2.

For standard imagers, the H operator is such that each line (corresponding to sky pixels contributing to one time sample) only contains a single non-zero value, meaning that s is actually the sky map convolved to the instrument's resolution and that the instrument samples the convolved sky with a single peak [START_REF] Tegmark | CMB mapping experiments: a designer's guide[END_REF][START_REF] Borrill | Madcap-the microwave anisotropy dataset computational analysis package[END_REF].

In the case of a bolometric interferometer, this assumption is not valid due to the multiple peaked shape of the synthesized beam (see figure 2) which makes it impossible to use the map-making algorithms usually developed for direct imagers. We use instead an algorithm that starts from an initial guess and then simulates iterative maps s i , where i is the iterative index2 . For each of these maps, we apply the bolometric interferometer acquisition model, taking into account the scanning strategy of the sky, and we construct TOD y i that are then compared to the data TOD y using a merit function that accounts for the noise in the TOD domain. In the case of stationary and Gaussian distributed noise, the maximum likelihood solution is reached by minimizing the χ 2 :

χ 2 ( s i ) = ( y -y i ) T • N -1 • ( y -y i ) (2.9)
where N is the covariance matrix of the noise. We minimize eq. 2.9 to find the best simulated sky map, ˆ s, using a preconditioned conjugate gradient method [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF][START_REF] Shewchuk | An introduction to the conjugate gradient method without the agonizing pain[END_REF]. This is jointly done for the IQU Stokes parameters and results in unbiased estimates of the maps as shown in figure 3. 3 Spectral dependence

The polychromatic synthesized beam

As can be seen in eq. 2.7, the synthesized beam is directly dependent on wavelength and this is shown in figure 2.

The off-axis angle (given by the primary beam of the pupils), the FWHM of the peaks (hence the resolution of the maps), and the angle on the sky between two peaks all depend linearly on λ. This dependence on wavelength can be exploited to achieve spectro-imaging capabilities. Within a wide band, the synthesized beam will be the integral of the synthesized beam of all the monochromatic contributions within the band resulting in a polychromatic synthesized beam. The left panel in figure 4 shows a cross cut of the monochromatic synthesized beam for 131 and 169 GHz while the right panel shows a polychromatic synthesized beam using 9 monochromatic synthesized beams. With a bolometric interferometer operating over a large bandwidth, for each pointing towards a given direction in the sky, one gets contributions from all the multiple peaks in the synthesized beam at all frequencies. As a result, we have both spatial and spectral information in the TOD. Precise knowledge of the synthesized beam along the frequency will then allow one to reconstruct the position and amplitude of the sky in multiple frequency sub-bands.

Spectro-imaging capabilities

The synthesized beam at two different frequencies ν 1 and ν 2 will be distinguishable from one another as long as their peaks are sufficiently separated. The angular separation between the two peaks ∆θ = c∆ν ν 2 ∆h (where ∆ν = ν 2 -ν 1 and ν = √ ν 1 ν 2 ) must be large enough to unambiguously distinguish the two peaks. We apply the Rayleigh criterion [START_REF] Rayleigh | XXXI. Investigations in optics, with special reference to the spectroscope[END_REF]:

c∆ν ν 2 ∆h c ν(P -1)∆h ⇔ ∆ν ν (P -1) (3.1)
where P is the number of pupils on a side of the square-packed pupil-array. A bolometric interferometer therefore not only has a resolution on the sky FWHM θ c ν(P -1)∆h , but also in electromagnetic frequency space ∆ν ν 1 P -1 . The map-making presented in section 2.4 can be extended in order to build, simultaneously, with the same TOD, maps at a number of different frequencies as long as they comply with the frequency separation given above. The iterative TOD y i can be written as:

y i = Nrec-1 j=0 H j ˆ s ij + n (3.2)
where H j describes the acquisition (convolution+pointing) operator with the synthesized beam at frequency ν j , ˆ s ij is the sky signal estimator at iteration i for the frequency ν j and N rec is the number of reconstructed sub-bands. Similarly, as in the map at a single frequency (figure 3), one can recover the maps s j by solving eq. 3.2 using a preconditioned conjugate-gradient method (see section 4 for corresponding simulations).

The QUBIC Full Instrument (FI) has two wide-bands centered at 150 and 220 GHz with ∆ν/ν = 0.25 and a 400-feedhorn array packed on a square grid within a circular area. The grid is 22 × 22 with the corners cropped (see figure 1). This corresponds to P = 22 and FWHMν ν ∼ 0.05. It is thus possible to reconstruct approximately 5 sub-bands in each of the initial bands of QUBIC. Note that this number should just be taken as an order of magnitude for the achievable number of sub-bands. The optimal number of sub-bands is determined through optimization of the signal-to-noise ratio.

Testing spectro-imaging on simple cases

We can use the QUBIC simulation pipeline to test the spectro-imaging capabilities of bolometric interferometry in a very simple manner. Some concepts and parameters used in simulations are defined in table 1.

Extended source reconstruction

The input map used for this example to simulate TOD is composed of zeros in each pixel of its 15 input frequencies, N in , and for the three Stokes components. a Skies are generated using PySM: Python Sky Model [START_REF] Thorne | The Python Sky Model: software for simulating the Galactic microwave sky[END_REF].

b Maps are projected using HEALPix: Hierarchical Equal Area Isolatitude Pixellization of sphere [START_REF] Gorski | HEALPix: A framework for high-resolution discretization and fast analysis of data distributed on the sphere[END_REF].

Table 1. Typical parameters used in acquisition, instrument and map-making to do an end-to-end simulation. A preconditionned conjugate gradient method is used for map-making.

extended regions with a high signal-to-noise ratio are added: a square centered at 141.6 GHz and a disk centered at 156.5 GHz. Map-making is done for 5 sub-bands centered at 134.6, 141.6, 148.9, 156.5 and 164.6 GHz, with a bandwidth of ∼ 6.8, 7.1, 7.5, 7.9 and 8.3 GHz respectively. The scan is performed with 8500 points randomly placed over a 150 square degrees sky patch. Noise is included at the TOD level. The first row in figure 5 shows the input sky maps spatially convolved at the QUBIC resolution. The second row shows the reconstructed maps after map-making onto five subbands. In the first, third and fifth sub-bands, where originally the signal is zero, structures corresponding to the signals of neighboring sub-bands appear. This is due to the fact that during the map-making process, leakage occurs from the frequencies where the monochromatic signal was located towards the neighboring sub-bands due to the Frequency Point Spread Function (FPSF) that will be studied in section 4.3.

Angular resolution

As an example, we used the end-to-end pipeline to simulate the reconstruction onto 4 subbands of a point source emitting with a flat spectrum in the 150 GHz wide band. Figure 6 shows the measured (red stars) and theoretical (blue dots) values of the FWHM at the central frequency of each sub-band. Theoretical values are obtained from a quasi-optical simulation [START_REF] O'sullivan | QUBIC VIII: Optical design and performance[END_REF] on HEALPix maps and corrected by pixel size and resolution [START_REF] Gamboa Lerena | Angular resolution at map level in the QUBIC instrument[END_REF]. The difference between measured and theoretical values are up to 0.5% which makes it acceptable. The real angular resolution will be determined once QUBIC is installed on the site using far-field observations (astronomical objects and/or calibration tower).

Frequency Point Spread Function (FPSF) characterization

In section 4.1 it was shown that the reconstructed map for a sub-band has a fraction of signal coming from neighboring bands (see figure 5). In order to quantify this effect, we simulate point source reconstruction to characterize the FPSF.

If we consider I in (ν) as the monochromatic intensity of the input map, and considering ideal map-making, then the intensity of the output map I out (ν) will be given by the convolution of input signal with the FPSF:

I out (ν) = (I in ⊗ FPSF) (ν). (4.1)
Thus, for a monochromatic source at ν in with intensity I in (ν) = δ(ν -ν in ), we can measure the FPSF by investigating the reconstruction of sources placed at different frequencies. We use a high-resolution frequency grid with N in = 48 which gives us a resolution of ∼ 0.78 GHz for the 150 GHz wide band (it would be ∼ 1.15 GHz for the 220 GHz band). This grid allows to improve the map within the spectral range and thus obtain more precise information on how the signal is reconstructed at the center and edge of each sub-band.

We performed 22 independent simulations of monochromatic point sources with high signal-to-noise ratio. We kept the spatial location of the point source unchanged and we varied its frequency ν in , covering a spectral range from 133 to 162.25 GHz. We present results for map-making onto 4 sub-bands with central frequencies at ν c = 135.5, 144.3, 153.6 and 163.6 GHz. Figure 7 shows the intensity, normalized to the input one, of the central pixel of the point source (blue dots) as a function of the separation between the input frequency and the central frequency of a sub-band, i.e. d νc = (ν in -ν c )/∆ νc , where ∆ νc is the width of the sub-band. The 4 sub-bands are overlaid in figure 7. We normalize by the width of the sub-band because the frequencies are logarithmically spaced. The FPSF in each sub-band are identical up to a scale factor ν. As expected from figure 5, we observe that the FPSF extends beyond a single sub-frequency and should be accounted for in the data analysis. This means that our reconstructed sub-bands are not independent from each other and we should expect noise correlation between sub-bands. Because the FPSF is negative in the nearest band we should expect the noise correlations to be negative between neighbouring sub-bands. This will be studied in section 5.1.

Galactic dust

We demonstrate spectro-imaging capabilities by trying to recover the frequency dependence of the dust emission with simulated observations towards the Galactic center. The sky maps contain IQU Stokes parameter components and the dust model is the one provided by PySM3, named d1 [START_REF] Thorne | The Python Sky Model: software for simulating the Galactic microwave sky[END_REF]. We simulate an observation in a sky patch of 15 degree radius. The parameters of the pipeline are set in such a way that the simulated instrument has a single focal plane operating either at 150 GHz or at 220 GHz with a 25% bandwidth each. The wide band TOD are formed through the sum of a number of monochromatic TOD throughout the wide bandwidth as shown in eq. 3.2. For this simulation we have used N in (see table 1) input maps covering the ranges from 137 to 162 GHz and from 192 to 247 GHz. From these wideband TOD, we are able to reconstruct several numbers of sub-bands using spectro-imaging. We have performed simulations with N rec = 1, 2, 3, 4, 5, 8 reconstructed sub-bands. Photon noise and detector NEP (see table 1) are added as white noise for each TOD. In each case, we perform a Monte-Carlo analysis to get several independent noise realizations and also a noiseless reconstruction that will be the reference. The result of this procedure at the map level for the I and Q components, for a given realization, is shown in figures 8 and 9. The two figures display a sky reconstructed in 5 sub-bands within the 220 GHz wide band. The residual map is the difference between a reconstruction and the noiseless reference. We can see that the Galactic dust is efficiently reconstructed in the 5 sub-bands as the residuals are compatible with pure noise. Note that the noise is not white but has spatial correlations due to the deconvolution with the multi-peak synthesized beam (see paper [START_REF] Hamilton | QUBIC I: Overview and Science Program[END_REF]). We also note that the residuals are higher on the edges than in the center of the sky patch. This is due to the higher coverage of the sky in the center due to the scanning strategy.

Instead of looking at the full map, the reconstructed intensity as a function of frequency can also be studied pixel by pixel. Figure 10 shows the intensity of the input sky convolved with the instrument beam, and the reconstructed intensity for a given pixel, considering 5 sub-bands in each wide band at 150 (red) and 220 (blue) GHz. We do not display the measured points and error-bars which are not good indicators of our uncertainties due to the highly anti-correlated nature of the covariance matrix (see figure 12 in the case of 3 subbands). Instead, we have performed a Monte-Carlo-Markov-Chain (MCMC) exploration of the amplitude and spectral-index of a typical dust model (modified Black-Body, see [START_REF] Irfan | Determining thermal dust emission from Planck HFI data using a sparse, parametric technique[END_REF]) accounting for the sub-bands covariance matrix. The fit is done separately for our two physical bands at 150 and 220 GHz. The 68% C.L. is shown in light colors in each case and represents the QUBIC measurements within this band using spectro-imaging. Note that the angular resolution of the maps improves with frequency (see section 4.2) and is not accounted for. As a result this simple analysis cannot be interpreted as a measurement of the dust spectral index which would require a more detailed analysis including the beam profile to properly infer the dust property in each sky pixel (such as [START_REF] Irfan | Determining thermal dust emission from Planck HFI data using a sparse, parametric technique[END_REF]).

The dust reconstruction is also studied in the angular power spectrum using the public code NaMaster [START_REF] García-García | Disconnected pseudo-Cl covariances for projected large-scale structure data[END_REF] which computes TT, EE, BB and TE spectra where T is the temperature and E, B are the two polarization modes. Spectra are computed from a multipole moment l = 40 to l = 2×N side -1 with N side = 256 the pixel resolution parameter for HEALPix maps. We compute Inter-Band Cross Spectra (IBCS), meaning that from N rec sub-band maps, one can compute N rec (N rec +1)/2 IBCS. Having independent noise realizations allows us to make IBCS crossing two realizations, so we eliminate the noise bias. BB IBCS for 3 sub-bands in each of the 150 and 220 GHz wide bands are shown in figure 11. We plot being the B-mode angular power spectrum. In this figure, the input theoretical dust spectra coming from the PySM model d1 are superposed to the reconstructed ones.

D l = l(l+1) 2π C l , C l

Noise characterization

For the map-making described in section 2.4, we added noise to the TOD which was composed of the detector Noise Equivalent Power (4.7 × 10 -17 W/ √ Hz) and photon noise (see table 1). The goal here is to study how close to optimal (in the statistical sense) is our spectro-imaging map-making. We will study the noise behaviour as a function of the number of reconstructed sub-bands N rec . This is done at 3 different levels: on the reconstructed maps, on the power spectra computed from the maps, and on a likelihood estimation for the tensor-to-scalar ratio r.

Noise behaviour in the sub-bands at map level

We performed simulations with 40 independent noise realizations and a noiseless simulation as a reference. After map-making, residual maps are computed by taking the difference between each simulation and the noiseless reference. For each pixel, we can compute the covariance matrix, over all the noise realizations, between the sub-bands and the Stokes parameters.

The reason why we treat the pixels separately, not computing covariance over them, is that the noise level varies with the position in the sky. This is due to the coverage of the sky by the instrument beam which is not uniform. Note that the QUBIC coverage is not trivial because of the multi-peak synthesized beam.

A correlation matrix for a given pixel, considering 3 sub-bands, is shown in figure 12 and we also show the average over pixels. It can be seen that for each Stokes parameter, residual sub-bands next to one another are anticorrelated and this is seen on every pixel. However, cross-correlations between Stokes parameters are negligible. This is why we can consider the 3 correlation matrices C Ip , C Qp and C U p separately.

I 0 I 1 I 2 Q 0 Q 1 Q 2 U 0 U 1 U 2 I 0 I 1 I 2 Q 0 Q 1 Q 2 U 0 U 1 U 2 Pixel 7
1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00

I 0 I 1 I 2 Q 0 Q 1 Q 2 U 0 U 1 U 2 I 0 I 1 I 2 Q 0 Q 1 Q 2 U 0 U 1 U 2
Averaged over pixels 1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00

Noise analysis using the power spectrum

In section 5.1 it was shown that a polychromatic interferometer has anti-correlations in neighbouring bands for each Stokes parameter. Furthermore, the spatial structure of noise is studied in detail in [START_REF] Hamilton | QUBIC I: Overview and Science Program[END_REF] using end-to-end simulations. This allows us to build a fast noise simulator that reproduces efficiently the noise behaviour in the reconstructed maps. In the following, the fast noise simulator will be used in parallel with end-to-end simulations as it allows us to improve the statistics while consuming much less computing time.

We characterize the noise behaviour of spectro-imaging using the power spectrum. As shown in section 4.4, from the maps we can compute power spectra using the public code NaMaster. From N rec bands, we compute the IBCS for each TT, EE, BB and TE power spectra. As we are interested in the noise, we compute the power spectra of the residual maps containing only noise. Figure 13 shows the IBCS computed for each noise realization in the case of 3 sub-bands at 150 GHz. As we plot D l the noise bias goes as l(l + 1). We find that the IBCS within the same band (ν 0 ν 0 , ν 1 ν 1 and ν 2 ν 2 ) are positively correlated. However the IBCS crossing 2 different bands (ν 0 ν 1 , ν 0 ν 2 and ν 1 ν 2 ) are anti-correlated.

The correlations are observed in greater detail by computing the correlation matrices. In figure 14, we show the correlation matrix between -bins and IBCS for BB angular power spectrum considering N rec = 3 sub-bands in the 150 GHz wide band. In this matrix, we see that anti-correlations, in blue in the matrix, only appear between the IBCS crossing 2 different bands (ν 0 ν 1 , ν 0 ν 2 and ν 1 ν 2 in the case of 3 sub-bands) and that the correlations between bins are negligible. The same behaviour is observed for TT, EE and TE spectra.

Nearly optimal performance of spectro-imaging

In order to assess, in a manner easy to interpret, how far from optimal is spectro-imaging, we study how the tensor-to-scalar ratio is constrained as a function of the number of sub-bands N rec . The division of the wide band into a number of sub-bands for spectro-imaging could have a detrimental effect on the estimate of the tensor-to-scalar ratio r. On the one hand, we would certainly like to make as many sub-bands as possible in order to constrain the foreground spectra in a very precise manner. However, on the other hand, there is an upperlimit to the achievable number of sub-bands, when the angular distance between peaks in the synthesized beams at different frequencies becomes smaller than the peak width (angular resolution), as explained in section 3.2. We therefore expect the performance of spectroimaging to degrade when projecting data onto too many sub-bands. In fact, even for a small number of sub-bands, spectro-imaging cannot be strictly optimal because the synthesized beams at different sub-frequencies do not form an orthogonal basis. We therefore expect a certain loss in signal-to-noise ratio when performing spectro-imaging. The higher the number of reconstructed sub-bands, the more overlap there is between the synthesized beam at each sub-frequency. This results in stronger degeneracy between sub-bands, hence a higher noise in the reconstruction. This is the price to pay for improved spectral resolution. As a result, one needs to find the best balance between performance and spectral resolution for a given scientific objective.

For that purpose, we project simulated data onto an increasing number of sub-bands N rec , calculate the corresponding IBCS and in each case we compute a likelihood to estimate the tensor-to-scalar ratio r combining all sub-bands accounting for their cross-correlations. The sky model is a pure CMB sky (including lensing but no Galactic foregrounds) with r = 0. From this method, we get the error on r at 68% confidence level for each number of sub-bands. This is presented in figure 15. We normalize by the case of "spectro-imaging" with just one band. Error bars are obtained from a Monte-Carlo analysis, varying the data in the likelihood according to their diagonal uncertainties. As expected, we observe a moderate degradation due to spectro-imaging in the sense that the constraints on r become less stringent when the number of sub-bands is greater than one. This degradation slowly evolves from 25% to 40% at 150 GHz and from 10% to 20% at 220 GHz when the number of sub-bands evolves from 2 to 5. The better performance at 220 GHz is not a surprise as our horns are slightly multimoded at 220 GHz (see [START_REF] O'sullivan | QUBIC VIII: Optical design and performance[END_REF] for details) resulting in a flatter primary beam, which, in turn, favours spectro-imaging because multiple peaks of the synthesized are higher in amplitude. It is possible to project onto as many as 8 sub-bands with a corresponding performance reduction due to the fact that synthesized beam peaks become too close with respect to their width, as explained in section 3.2.

This study demonstrates that, although not optimal from the noise point of view, spectro-imaging performance remains close to optimal for up to 5 bands, providing extra . Uncertainties (68% C.L. upper limits) on the tensor-to-scalar ratio r obtained by combining an increasing number of sub-bands, normalized to that of one band, for a pure r = 0 CMB (with lensing). The slow increase of the uncertainty on r with the number of sub-bands illustrates the moderate sub-optimality of spectro-imaging and shows that we can use up to 5 sub-bands with only 40% degradation at 150 GHz (and only 20% degradation at 220 GHz). It is possible to achieve 8 sub-bands but with more significant degradation. spectral resolution that can be key for constraining foreground contamination with realistic models for which the spectrum might not be a simple power law, work being in progress on this. The appropriate balance between spectral resolution and noise performance can be adjusted for each specific analysis thanks to the fact that spectro-imaging is done entirely in post-processing.

Conclusion

In this article, we have shown how the new technique of Bolometric Interferometry offers the possibility to also perform spectro-imaging. This makes it possible to split, in post processing, the wide-band observations into multiple sub-bands achieving spectral resolution. To illustrate this method, we apply it to the case of the QUBIC instrument soon to be installed at its observation site in Argentina.

After having presented the design, concept, and mathematical aspects of the instrument and the spectro-imaging technique we have illustrated it on simple cases: monochromatic point sources, spatially extended sources, and sky maps with frequency-dependent emission such as Galactic dust. We have shown our ability to have increased spectral resolution with respect to the physical bandwidth, considering a full sky patch but also at the level of individual pixels. We studied the signal and noise behavior using Monte-Carlo simulations for an instrument like QUBIC which shows spatial and spectral correlations. We have quantified the loss of statistical performance for the measurement of the tensor-to-scalar ratio when increasing the number of sub-bands and have shown it to be moderate up to 5 sub-bands.

The precise measurement of foreground contaminants is essential for the detection of primordial B-modes. Foregrounds have spectral properties distinct from the CMB which leads to the conclusion that only a multichroic approach enables the measurement and subtraction of foreground contamination. This is usually done in classical imagers through detectors operating at distinct frequencies, each of them being wide-band in order to maximize signal-to-noise ratio. However, constraining foregrounds with such data relies on extrapolation between distant frequency bands, which may miss non-trivial variations of the spectral behaviour of complex foregrounds such as multiple dust clouds in the line of sight. In particular, scenarios where dust exhibits a certain level of decorrelation between widely separated bands, or with non constant spectral indices would be impossible to be identified with a usual wide-band analysis. Spectro-imaging could put significant constraints on such scenarios. This is being studied in detail by the QUBIC Collaboration and will be presented in the near future.

In summary, spectro-imaging improves spectral resolution within a wide physical band, while nearly preserving the optimal performance of the analysis. It may therefore become a key technique for detecting the elusive B-mode polarization of the CMB.

Figure 1 .

 1 Figure 1. Left: Feedhorn array of the Full Instrument(FI). Right: Intensity pattern on the focal plane composed by 992 bolometers for an on-axis point source in the far field emitting at 150 GHz with all horns open. The color scale in intensity is arbitrary.

KFigure 3 .

 3 Figure3. Result of the map-making for IQU Stokes parameters for a bolometric interferometer pointing in a 15 degree radius sky patch containing only CMB. The first column is the input sky convolved at the resolution of the instrument using a Gaussian with a FWHM equal to 0.4 degrees. The second column is the sky reconstructed by map-making. The last column is the difference between both. This simulation was obtained with the QUBIC full pipeline. The noise was scaled to 4 years of observations.

Figure 4 .

 4 Figure 4. Left: Monochromatic synthesized beam (MonoSB ) for 131 and 169 GHz. Each synthesized beam is modeled according to eq. 2.7 for a square array of 20 × 20 pupils separated 14 mm apart at 131 and 169 GHz. The primary beam at each frequency is shown by a dashed line. Right: Polychromatic beam (PolySB, black line) as result of the addition of 9 monochromatic synthesized beams (5 of them are shown in colored lines) spanning our 150 GHz band (131 to 169 GHz). We sample the continuous frequency band with discrete frequencies.

Figure 5 .

 5 Figure 5. Map-making for a sky full of zeros with two extended, monochromatic regions centered at 141.6 (square) and 156.5 GHz (disk). Each column corresponds to one sub-band. The first row shows the input sky maps spatially convolved at the QUBIC resolution. The reconstructed maps using the QUBIC pipeline are shown in the second row. The units are arbitrary.

Figure 6 .

 6 Figure 6. Red stars represent the angular resolution measured for a single point source emitting in the broad band after map-making onto 4 sub-bands. Blue dots are the theoretical values expected for each central frequency of the sub-bands.

Figure 7 .

 7 Figure 7. Blue dots are the intensity of the central pixel of the point source for each sub-band and each simulation as a function of d νc . Red line is the fit model using a polynomial of degree 12 (in light-red 3-σ uncertainty band). Strong-grey band represent the limits of sub-band centered in frequency ν c and soft-grey represents the limits of the adjacent sub-bands. Fit and plot were done using lmfit package.

Figure 10 .

 10 Figure 10. Intensity as a function of the frequency for N rec = 5 sub-bands in each wide band at 150 (red) and 220 (blue) GHz for a given pixel. The grey regions correspond to the unobserved frequencies outside our physical bands. Two sky pixels are shown as red stars, one in a patch centered at the Galactic center and one in the patch that QUBIC plans to observe centered in [0, -57 deg]. Red and blue dots: Input sky convolved with the instrument beam. In both cases are shown in light color the 68% C.L. regions for a modified black-body spectrum reconstructed with a MCMC from our simulated measurements and sub-band covariance matrices (see figure 12 for the case of 3 sub-bands). Maps are in µK CMB and N side = 32.

Figure 11 .

 11 Figure 11. BB Inter-Band Cross Spectra (IBCS) for 150 GHz (left) and 220 GHz (right) computed from reconstructed maps obtained with end-to-end simulations. For each IBCS, we cross-correlate 2 sub-bands with central frequencies in GHz shown in the legend of each plot. Dashed lines are the IBCS of the input sky that contains only Galactic dust. The dots with error bars show the mean and the standard deviation over 20 IBCS. Each IBCS is made with 2 maps with independent noise realizations to eliminate the noise bias.

Figure 12 .

 12 Figure 12. Correlation matrices between frequency bands ν 0 = 137 GHz, ν 1 = 149 GHz and ν 2 = 162 GHz and I, Q, U Stokes components obtained from 40 end-to-end simulations. Left: example for a given pixel. Right: The average over pixels. Blue means anti-correlations while red is for positive correlations.

Figure 13 .

 13 Figure 13. BB Inter Band Cross Spectra on the residual maps containing only noise for 3 sub-bands in the wide 150 GHz band, centered at 137, 149 and 162 GHz. Dots and error bars show average and standard deviation over 1000 independent noise realisation IBCS computed with the fast noise simulator.

Figure 14 .

 14 Figure 14. Correlation matrix between -bins and IBCS for BB angular power spectrum considering 3 sub-bands at ν 0 = 137, ν 1 = 149 and ν 2 = 162 GHz. For example, ν 0 ν 1 is the IBCS between frequencies ν 0 and ν 1 . In each black square we show the correlations between the 16 -bins used to compute the IBCS as in figure13.

Figure 15

 15 Figure[START_REF] Suzuki | The litebird satellite mission: Sub-kelvin instrument[END_REF]. Uncertainties (68% C.L. upper limits) on the tensor-to-scalar ratio r obtained by combining an increasing number of sub-bands, normalized to that of one band, for a pure r = 0 CMB (with lensing). The slow increase of the uncertainty on r with the number of sub-bands illustrates the moderate sub-optimality of spectro-imaging and shows that we can use up to 5 sub-bands with only 40% degradation at 150 GHz (and only 20% degradation at 220 GHz). It is possible to achieve 8 sub-bands but with more significant degradation.

  Two monochromatic Parameter Details N in Number of Input or True maps (in µK) used to simulate a broadband observation (TOD). Each map represents a sky at a specific frequency ν j for the IQU Stokes parameters a b . Possible values: 15, 16 or 48. N rec Number of sub-bands reconstructed from a single broadband observation. In all simulations N rec is a divisor of N in . N conv Number of convolved maps equal to N rec . Each of these maps is obtained convolving the N in input maps at the QUBIC spatial resolution corresponding to that input frequency and then averaging within the reconstructed sub-band. NEP det Detector Noise Equivalent Power added as white noise. Value: 4.7 × 10 -17 W/ √ Hz. p γ Photon noise added as white noise in time-domain, calculated from the atmospheric emissivity measured in our site, as well as emissivities from all components in the optical chain. The value is different for each detector because of their different illumination by the secondary beam B sec . The average value at 150 GHz is 4.55 × 10 -17 W/ √ Hz and 1.72 × 10 -16 W/ √ Hz at 220 GHz. θ Radius of sky patch observed in simulations. Value: 15 degrees. pointings Number of times that the instrument observes in a given sky direction aligned with the optical axis. Values > 10 4 .

  at 150 GHz and scaled proportionally to frequency. Measurements were made

		0.42 0.43 0.44					FWHM measured FWHM theoretical
	FWHM [deg]	0.38 0.39 0.40 0.41
		0.37
		0.36	135	140	145	150 [GHz]	155	160	165

r = 0 at the optical center of the focal plane.

The software uses the massively parallel libraries[START_REF] Chanial | PyOperators: Operators and solvers for high-performance computing[END_REF] developed by P. Chanial pyoperators (https://pchanial.github.io/pyoperators/) and pysimulators (https://pchanial.github.io/pysimulators/).

Acknowledgements QUBIC is funded by the following agencies. France: ANR (Agence Nationale de la Recherche) 2012 and 2014, DIM-ACAV (Domaine d'Interet Majeur-Astronomie et Conditions d'Apparition de la Vie), Labex UnivEarthS (Université de Paris), CNRS/IN2P3 (Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules), CNRS/INSU (Centre National de la Recherche Scientifique/Institut National des Sciences de l'Univers). Italy: CNR/PNRA (Consiglio Nazionale delle Ricerche/Programma Nazionale Ricerche in Antartide) until 2016, INFN (Istituto Nazionale di Fisica Nucleare) since 2017. Argentina: MINCyT (Ministerio de Ciencia, Tecnología e Innovación), CNEA (Comisión Nacional de Energía Atómica), CONICET (Consejo Nacional de Investigaciones Científicas y Técnicas).