J-Ph Bernard

A Bernard

H Roussel

I Choubani

D Alina

J Aumont

A Hughes

I Ristorcelli

S Stever

T Matsumura

S Sugiyama

K Komatsu

G De Gasperis

K Ferrière

V Guillet

N Ysard

P Ade

P De Bernardis

N Bray

B Crane

J-P P Dubois

M Griffin

P Hargrave

Y Longval

S Louvel

B Maffei

S Masi

B Mot

J Montel

F Pajot

E Pérot

N Ponthieu

L Rodriguez

V Sauvage

G Savini

C Tucker

F Vacher

I Chubani

I Ris- Torcelli

G De

Performance of the polarization leakage correction in the PILOT data

Keywords: PILOT, Interstellar Dust, Polarization, Far Infrared, systematic effects

. In this paper, we present data processing steps that were not presented in [1] and that we have recently implemented to correct for several remaining instrumental effects. The additional data processing concerns corrections related to detector cross-talk and readout circuit memory effects, and leakage from total intensity to polarization. We illustrate the above effects and the performance of our corrections using data obtained during the third flight of PILOT, but the methods used to assess the impact of these effects on the final science-ready data, and our strategies for correcting them will be applied to all PILOT data.

Introduction

Interstellar dust grains account for ≃1% of the mass of the interstellar medium (ISM). They are involved in different important processes such as photo-electric heating of the neutral interstellar gas, cooling in dense star-forming regions and the formation of molecules, including H 2 , on their surfaces. Dust emission can be used to trace the structure of the interstellar medium (ISM) in the Milky Way and in the local Universe (e.g., [2,3,4]). The thermal dust emission can be modeled using a modified blackbody spectrum in the infrared to submillimeter wavelength range, but physically motivated models remain a subject of debate, since the exact nature of the dust grains is still largely unknown. Understanding dust emission polarization is also important to devise foreground subtraction strategies for CMB experiments.

ISM dust grains absorb starlight in the visible and ultraviolet, which heats them to temperatures of ≃17 K in the diffuse ISM in the solar neighborhood in our Galaxy [5]. The polarization of thermal dust emission is believed to arise from the irregular shape of dust grains, and their alignment. The global alignment is believed to be the result of fast grain rotation and relaxation processes slowly bringing the grain minor axis onto the local magnetic field direction (e.g., [6,7]). The grain thermal emission being stronger along the long axis of the grain, the global partial alignment causes a fraction of the thermal emission to be linearly polarized in a direction orthogonal to the magnetic field direction as projected on the sky. For the same reason, non-polarized starlight passing through the ISM with aligned dust grains also becomes polarized, with preferential absorption along the long axis of the grains leading to extinction in the visible and the near-infrared (NIR) being polarized parallel to the magnetic field lines.

First measurements of the polarized extinction in the visible and NIR date from the 1960s (see large catalogs such as in [8]). These studies allowed accurate measurements of the spectral shape of the polarized extinction curve, also known as the Serkowski law ( [9]), which is an efficient way of constraining the size distribution of dust grains. Measurements of the thermal dust emission in polarization are more recent. The balloon experiment Archeops ( [10]) mapped the polarized dust emission at 353 GHz with ∼ 13 ′ resolution over ∼ 20% of the sky. These measurements indicated high po-larization levels (up to 15%) in the diffuse ISM. More recently, the Planck satellite mapped the polarized emission over the whole sky in 7 spectral bands in the wavelength range 850 µm (353 GHz) to 1.0 cm (30 GHz) [11]. At the highest frequencies, thermal dust dominates the polarization signal, while low frequencies are typically dominated by polarized synchrotron emission. Analysis of the polarized thermal dust emission at 353 GHz ( [11]) indicated a good correlation with polarized extinction. [12] showed that the overall thermal dust polarization fraction is only a few percent of the total dust emission over most of the sky, but confirmed the existence of highly polarized regions at high galactic latitudes with polarization fractions reaching up to 22%. These studies also demonstrated that the thermal dust polarization fraction varies by large factors on small scales. These variations appear linked to the total gas column density, with dense regions exhibiting lower polarization, and to the structure of the magnetic field, with regions showing the most B-field rotation on the plane of the sky also being the least polarized. This latter behavior was shown to be consistent with predictions of MHD models of the ISM (see [13]). As a consequence of the above studies, the polarized dust thermal emission is now recognized as a dominant foreground contaminant to the observation of the Cosmic Microwave Background (CMB) polarization (see [14]).

Several other facilities allow observations of the thermal dust emission from airborne and ground-based telescopes. The HAWC+ instrument on SOFIA has polarization capabilities in 4 bands from 53 µm to 214 µm ( [15]). The SCUBApol instrument on JCMT [16] can also map thermal dust polarization at 850 µm. The NIKA2 instrument on the IRAM 30m telescope [17] can be used to measure polarization at 260 GHz (1.1 mm) with angular resolution of 10". Finally, the ALMA interferometer allows polarization measurements in band 7 (350 GHz) with very high angular resolution ( [18]). The BLASTPol instrument ( [START_REF] Fissel | Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy V[END_REF]) measures polarized dust emission in 3 bands from 250 µm to 500 µm. In most cases, these facilities are limited in sensitivity to observations of very bright regions and/or suffer from atmospheric absorption or emission fluctuations. Because they can only map fields of view that are limited in size, at much better angular resolution than Planck, a comparison of their results with those of Planck for the same region is at best very difficult, sometimes impossible.

Measuring the spectral and spatial variations of polarized dust emission provides a potentially powerful constraint on the physics of dust grains (see for instance [20]), and is crucial to accurately separate the contribution of the polarized Galactic foreground from the CMB signal. To date, spectral variations of dust polarization have been only poorly constrained by observations. [21] established the first reliable measurement of the spectral variations over the Planck frequency range, using the average dust emission over a care-fully selected fraction of the sky. This study concluded that the polarization fraction is roughly constant across 353 GHz to 100 GHz, with some indication (at the 3σ level) that the polarization fraction decreases with decreasing frequency. This measurement of the spectral shape of the dust polarization fraction is extremely challenging due to the decreasing brightness of dust emission at low frequencies and the increasing contribution of polarized synchrotron emission and unpolarized sources such as spinning dust emission and free-free. At frequencies above 353 GHz, most existing measurements have been obtained by large ground-based or airborne telescopes, which can only map very restricted regions around bright sources. Differences in resolution and the differential scale filtering necessary to subtract atmospheric emission complicate an co-analysis of these measurements and the Planck data. As a consequence, there is so far very little information available about the polarized SED of thermal dust emission. A key objective of the PI-LOT mission is to improve our understanding of the thermal dust polarization signal, by measuring it at higher frequencies than Planck in the far-infrared, at an angular resolution and spatial coverage that enables a robust co-analysis with the Planck data.

Measurements of astrophysical polarization are difficult because the signal is extremely weak. Most, if not all, of the instruments mentioned above have encountered difficulties in accurately measuring polarization at low intensities due to systematic instrumental effects. Some of these effects result from well-understood phenomena, such as imperfect inter-calibration of detectors, inaccurate correction for time constants of detectors and for electronic cross-talk, ADC conversion, unmasked glitches, etc. Other systematic effects have been discovered during data processing, such as the spurious contributions from molecular gas spectral lines in the signal [22] and bandpass mismatch between detectors [23,24], both of which were encountered in the Planck data and required dedicated complex treatment [25]. Another example is the effect of the Gore-Tex membrane in front of the JCMT which requires special treatment ( [START_REF] Friberg | Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy IX[END_REF]). Recently, a leakage from intensity to polarization has been identified by several experiments including NIKA2 ( [START_REF] Ritacco | [END_REF][START_REF] Ajeddig | European Physical Journal Web of Conferences[END_REF][START_REF] Ajeddig | [END_REF]), and HAWC+ ( [30]) as a clear limitation to the accuracy of polarization measurement. This effect appears to originate from imperfections of the optical systems that lead to asymmetries in the optical ray propagation through the instrument, producing artificial polarization signal on unpolarized sources. The exact origin is not fully understood and may be instrument dependent.

In this paper, we present the method used to correct for the polarization leakage in the PILOT data and evaluate its performance. We describe two other systematic effects that have an electronic origin, -detector cross-talk and a readout electronics memory effect -that affect the PILOT point spread function (PSF) and must be addressed prior to the leakage correction. We give a short description of the instrument and the flights in Sect. 2 and Sect. 3 respectively. In Sect. 4 we present observations of Jupiter obtained during flight 3, which show the effect of read-out latency, cross-talk and leakage. We use the Jupiter data to characterize and correct for the above systematic effects. We show residual maps to assess the uncertainties associated with residual systematic effects after correction and measure the performance of the leakage correction, in Sect. 4.4. We summarize our conclusions in Sect 5. A complete description of the PILOT instrument is available in [START_REF] Bernard | the PILOT team[END_REF]. Here, we only give a brief description for completeness. Table 1 summarizes the main characteristics of the instrument.

The PILOT instrument

The telescope optics comprises an off-axis paraboloid primary mirror (M1) with diameter of 0.83 m and an off-axis ellipsoid secondary mirror (M2). The combination respects the Mizuguchi-Dragone condition to minimize depolarization effects (see [START_REF] Bernard | the PILOT team[END_REF]32]). All optics following M1, including M2, are cooled to a cryogenic temperature of 2 K.

Following the Gregorian telescope, the beam is folded using a flat mirror (M3) towards a re-imager and a polarimeter. Two lenses (L1 and L2) are used to re-image the focus of the telescope on the detectors. A Lyot-stop is placed between the lenses at a pupil plane that is a conjugate of the primary mirror. A rotating Half-Wave Plate (HWP), made of sapphire, is located next to the Lyot-stop. The bi-refringent material of the HWP introduces a phase shift between the two orthogonal polarization components of the incident light. A polarization analyzer consisting of parallel metallic wires is placed at a 45 o angle in front of the detectors, in order to transmit one polarization to the transmission (TRANS) focal plane and reflect the other polarization to the reflection (RE-FLEX) focal plane. Observations at two or greater different HWP angles allow us to reconstruct the Stokes parameters I, Q and U as described in Sect. 2.1. Each of the TRANS and REFLEX focal planes includes 1024 bolometers (4 arrays of 16 X 16 pixels). They are cooled to 300 mK by a closed cycle 3 He fridge. The detectors were developed by CEA/LETI for the PACS instrument on board the Herschel satellite.

In order to reconstruct the pointing of the instrument, we use the Estadius stellar sensor developed by CNES for stratospheric applications and described in [START_REF] Montel | 22nd ESA Symposium on European Rocket and Balloon Programmes and Related Research[END_REF]. This system provides an angular resolution of a few arcseconds, which is required to optimally combine observations of the same part of the sky obtained with various polarization analysis angles. A key feature of Estadius is that it remains accurate even with fast scan speeds (up to 1 o /s). An internal calibration source (ICS) is used inflight to calibrate time variations of the detector responses. This device is described in [START_REF] Hargrave | Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series[END_REF][START_REF] Hargrave | IR Space Telescopes and Instruments[END_REF]. The source is located behind mirror M3 and illuminates all detectors simultaneously. It is driven using a square modulated current. The current and voltage of the source are measured continuously during flight, in order to monitor the power dissipated in the source.

Polarization measurements

Assuming a perfect HWP, the PILOT measurements m are related to the input Stokes parameters I, Q, U of partially linearly polarized light through

m = R xy T xy × [I ± Q inst cos 4ω ± U inst sin 4ω] + O xy , (1)
where R xy and T xy are the detectors response and optical transmission respectively, and O xy is an arbitrary electronics offset. For the configuration of the HWP and polarizer in the instrument, ω is the angle between the HWP fast axis direction and the horizontal direction measured counterclockwise as seen from the instrument. The ± sign is + andfor the REFLEX and TRANS arrays respectively (see [START_REF] Bernard | the PILOT team[END_REF]). Note that, with the above conventions, Q inst and U inst are defined with respect to instrument coordinates in the IAU convention, with Q inst =0 for vertical polarization. For PILOT, ω is related to a mechanical HWP position called HW P pos , which can be varied continuously over the range 1 ≤ HW P pos ≤ 8 as

ω = 87.25 o -(HW P pos -5) × 11.25 o , (2)
allowing the HWP fast axis to vary by approx. ±45 o around the vertical direction. When referring to the sky polarization Q and U, Eq. 1 becomes

m = R xy T xy × [I ± Q cos(2θ ) ± U sin(2θ )] + O xy , (3)
where θ = 2 × ω + φ is the analysis angle, φ is the time varying parallactic angle measured counterclockwise from equatorial north to Zenith for the time and direction of the current observation, and Q and U are in the IAU convention with respect to equatorial coordinates. In practice, maps of Q and U are derived from observing the same patch of sky with at least two values of the analysis angle taken at different times in general. Inversion to derive sky maps of I, Q and U can be done through polarization map-making algorithms (see for instance [36]). The light polarization fraction p and polarization direction ψ are then defined as:

p = Q 2 + U 2 I (4)
and

ψ = 0.5 × arctan(U/Q). (5)
3 The PILOT flights and observations PILOT is carried to the stratosphere by a generic gondola suspended under an open stratospheric balloon through a flight chain, with a helium gas volume of ∼ 800 000 m 3 at ceiling altitude. The flights are operated by the French National Space Agency (CNES) with launch campaigns involving several international balloon experiments (up to six per campaign). At ceiling altitude, the instrument can be pointed towards a given sky direction using the gondola rotation around the flight chain and rotation of the instrument around the elevation axis (see [START_REF] Bernard | the PILOT team[END_REF]). Scientific observations are organized into individual observing tiles (also called observations for short) during which a given rectangular region of the sky is scanned by combining the azimuth and elevation rotations. The flight plan is built taking into account the various observational constraints such as the visibility of astronomical sources, the minimum allowed angular distance between the instrument optical axis and bright sources such as the Sun or the Moon, elevation limits due to the presence of the Earth at low elevations and the balloon at high elevations. The expected performance of the instrument is taken into account when establishing the flight plan, in order to distribute the observing time according to the science objectives, and to evenly distribute both the polarization analysis directions (angle θ in Eq. 3) and the scanning directions for any given astronomical target.

The first two flights of the PILOT experiment took place from the launch-base facilities at the airports of Timmins (Ontario, Canada) in 2015, and Alice Springs (Australia) in 2017, respectively. A detailed description of the characteristics of these flights and the corresponding observations are presented in [37]. In this paper, we focus on the instrument performance during the third PILOT flight.

Performance during flight#3

The third flight of the PILOT instrument took place from Timmins on September 24 2019, as part of a balloon experiment launch campaign led by CNES and the Canadian Space Agency (CSA).

The flight lasted approximately 26 hr, during which 21 hr of scientific observations were obtained. The launch took place at 5:36 AM local time. The experiment reached ceiling altitude about 2.3 hr after launch. The instrument reached an altitude of 39 km, slowly decreasing to 37 km during the first day of the flight. The altitude decreased to 34 km during the night due to the lower buoyancy force of the balloon. During the second day, the altitude rose again, reaching 37.5 km just before the gondola was dropped in Quebec. The temperatures of the focal planes evolved slightly with altitude during the ceiling period and remained in the range 296.5 to 297 mK and 300 to 301 mK for the TRANS and REFLEX focal planes respectively during the day, and ≃ 297.5 mK and ≃ 301.5 mK during the night. The higher nocturnal temperatures are due to less efficient pumping on the He bath. Out of the eight bolometer arrays, array #1 (TRANS), array #5 (TRANS) and array #6 (TRANS) were not operational during flight#3. The footprint of the available arrays on the sky is shown in Fig. 6.

The balloon followed a trajectory towards the north-east during most of the flight. We successfully used the two telemetry antennas located in Timmins and Chibougamau. The gondola was recovered about 900 km north-east of the launch site, north of Saguenay, Quebec. The gondola was brought back to the Timmins base using a helicopter and a truck.

The gondola and the instrument suffered no major damage from landing or recovery, which was later confirmed by a thorough inspection following the return of the instrument to France. The astronomical sources targeted during flight#3 are listed in Table 2. In the following, we concentrate on the analysis of the data obtained on Jupiter, which we use to characterize systematic effects.

Systematic effects

In this section we describe three instrumental systematic effects of the PILOT data, not addressed in [1]. Two of these effects are related to the readout electronics of the PILOT detectors, which we refer to as cross-talk and read-out latency. The third effect is produced by the optics of the instrument, which we refer to as leakage (see Sect. 4.3). Here, we describe the manifestation of each of these effects on the instrument Point Spread Function (PSF) as observed on Jupiter during flight#3, how we measured the parameters used in the correction methods, how the corrections were performed and the overall performance of the corrections, as measured on the Jupiter observations.

During each flight of the instrument, we observed planets, which can be considered point sources at the resolution of PILOT. These observations can be used to assess the optical quality through a measurement of the PSF. During Flight#3, we observed Jupiter at its maximum elevation of ≃ 17 • during about 30 min at the start of the flight. We obtained eight individual maps using eight positions of the HWP, sampling the available analysis range uniformly. The maps were obtained at two different scanning angles to enable low frequency noise removal. The data were corrected for the responses calculated on the residual atmospheric signal and the ICS calibration signals, and corrected for the effect of the detector time constants through deconvolution, as described in [1]. The signal was then processed using the Scanamorphos map-making software described in [38] and as used in its polarization version in [37] to produce maps of the Stokes parameters I, Q and U and the corresponding variances and co-variances. Note that these maps can also be obtained in instrument coordinates, also referred to as crosselevation and elevation, by setting the parallactic angle to zero in Eq. 3. This representation is optimal to reveal and characterize effects that project in the focal plane of the instrument, since it avoids blurring through sky rotation.

In order to produce PSF maps that are sufficiently accurate to be used for leakage correction, we constructed Jupiter maps with a pixel size of 6 ′′ . We generated individual maps for each detector array and for each individual observation of the planet. We use these maps for assessing potential temporal or focal plane variations of the systematic effects.

Figure 1 shows the maps of Jupiter in instrumental coordinates where all the data from flight#3 have been used.

Below, we use these maps to investigate systematic effects affecting the PSF in polarization, and to measure the parameters used in the correction method.

Cross-talk

Figure 1 shows the total intensity map of Jupiter obtained during flight#3 in instrument coordinates. A blurred linear extension is clearly visible across the PSF from the lowerleft to upper-right of the image. As illustrated in Fig. 2, this direction corresponds to the orientation of the individual pixel lines on the arrays, which are rotated 45 • with respect to instrument coordinates. The readout electronics of the PILOT detectors is such that the signals from bolometers along each line are transferred simultaneously to a buffer unit (BU) with 16 registers for amplification. We interpret the observed effect as cross-talk between pixels along a given detector line, which could happen within the BU. As the strong signal from the peak of the optical PSF falls on a given pixel of the array, part of its intensity is transferred through cross-talk to other pixels along the corresponding detector line, producing the observed pattern. This is illustrated in Fig. 2 which shows the projection of the array footprint on the Jupiter map.

We model the cross-talk as the transfer of a fraction f CT (i, k) of the signal from pixel i to pixel k along line j. As the transfer is likely to occur in the buffer unit which is common to all lines, we further assume that f CT (i, k) is independent of line number j. We therefore subtract cross-talk following

d ′ i j = d i j -∑ k =i f CT (i, k) × d k j + ∑ k =i f CT (i, k) × d i j , (6)
where d i j and d ′ i j are the signal before and after cross-talk subtraction respectively, and the summations are carried out over all pixels along the line under consideration. The two terms correspond to the signal received and given by the considered pixel. We also assume that, cross-talk being an induction effect, it is symmetric with f CT (i, k) = f CT (k, i).

We searched for possible variations of the cross-talk coefficient along detector lines. For this, we analyzed jointly the recordings of pixels receiving a strong glitch (normally masked out during the processing) and the stacked signal of simultaneous recordings of other pixels of the same line. We correlated the stacked signals from the main glitched pixel and the cross-talk pixel. This analysis did not produce convincing evidence for variations of f CT (i, k) along the pixel lines. We therefore assume that f CT (i, k) does not vary across a given array and we search for a single value of the pixel-to-pixel cross-talk coefficient f CT for each array.

In order to determine f CT for each array, we defined a cross-talk region and a reference region in each image of Jupiter (see Fig. 2). Both regions share the same average distance from the planet so that they would have the same brightness in the absence of cross-talk but are centered on a regions of high and low cross-talk signals respectively. The value of f CT for each array was found through a χ 2 minimization of the signal difference between the cross-talk and the reference regions in images of Jupiter obtained with each array separately. At each iteration of the minimization, the cross-talk signal was subtracted from the timeline using Eq. 6 and a new image was produced. The resulting values of the cross-talk parameters f CT are given in Tab. 3. The derived values appear consistent between arrays and at a level just below 10 -3 . The values are comparable between arrays and there appear to be no particular similarities between parameter values for arrays sharing the same BU.

The Jupiter map after correction of the cross-talk using the parameters given in Tab. 3 and Eq. 6 is shown in Fig. 4. Following the correction for cross-talk, the amplitude of the effect as measured in the regions defined above is about 0.9% of the PSF peak value, significantly smaller than the initial value of 4.8%.

Read-out latency

Figure 2 shows a map of Jupiter obtained during flight#3 where the Scanamorphos glitch detection was inhibited during processing. The map is overlaid with the footprint of one of the PILOT bolometer arrays. The image clearly shows some positive signal located precisely one array away from Jupiter along the column direction of the array. This fake source appears in the data stream of bolometers located on line 1 of each array, only when a bright source is present on the same column on line 16. This effect had already been seen clearly in calibration data, when a bright source was moved over all pixels of all arrays ( [START_REF] Misawa | the PILOT team[END_REF]). It is attributed to latency in the time-multiplexed readout electronics, which we refer to as read-out latency. This effect transfers some of the signal from one readout to the next along the readout direction (see Fig. 2), including when the readout goes from line 16 back to line 1, which creates the fake positive source in the map. The effect is also seen as a negative shadow of the cross-talk signal described in Sect. 4.1, which indicates that the read-out latency effect is mostly negative during transfer across the array and positive when readout is reset to line 1. The fact that we see the effect of the read-out latency on the cross-talk signal also shows that the read-out latency happens after the cross-talk in the detection chain, and as a consequence it needs to be corrected before cross-talk in the data processing.

We model the read-out latency effect as the transfer of a fraction f RL (i, j) of the signal from readout j to readout Fig. 3: Parameters for the read-out latency derived for each array, plotted as a function of the readout column. The black curve shows the parameter values between adjacent pixels along the column. The blue curve shows the parameter values when readout returns from line 16 to line 1. j + 1 along a given column i. We correct for the effect in the timelines using

d ′ i j = d i j -f RL (i, j -1) × d i j-1 + f RL (i, j + 1) × d i j+1 , (7)
where all readouts have been ordered with time of acquisition. The first term corresponds to the signal received by the considered pixel from the previous readout and the second term corresponds to the signal given to the next readout. We further assume the same value for f RL (i, j) between all consecutive readouts, except for multiples of 16 ( j = n × 16) with value f RL (i, 16) in Eq. 7. In order to measure the parameter f RL (i, j) for each readout column i and array for j = n × 16 , we constructed images of Jupiter in instrument coordinates using only the signal from that individual column of that individual array, using a timeline corrected according to Eq. 7. We optimized the value of f RL (i, j) in order to minimize the difference between the average signal in two rectangular boxes located on both sides of the crosstalk extension around Jupiter, as shown in Fig. 2. In order to measure f RL (i, 16), we performed a similar minimization but minimizing the signal in a square box centered on the fake source in the images as shown in Fig. 2. In both cases, the minimization was performed using the χ 2 minimization algorithm implemented in the IDL routine mpfitfun.

The values derived for f RL (i, j) are shown in Fig. 3 for each column of each array. The values are generally negative, while f RL (i, 16) is generally much smaller in absolute value but mostly positive.

The Jupiter map after correction of the read-out latency using the parameters shown in Fig. 3 and Eq. 7 is shown on Fig. 4. Note that this correction produces a significant shift of the planet peak position. Since we compute sky coordinates of each data sample using the data from the Estadius stellar sensor and correcting for the offset between the sensor and the PILOT instrument optical axis using the position of observed bright sources, we recompute coordinates following that correction, which we use for the rest of the data processing and analysis.

Leakage

Figure 1 shows the I, Q and U maps obtained on Jupiter during flight#3 projected in instrument coordinates (elevation and cross-elevation). While we expect Jupiter to show no polarization, the maps clearly show some residual Q and U structures at the location of the planet.

Figure 5 shows the cumulative profiles of the I, Q and U as a function of the integration radius in the maps of Fig. 1, as well as the corresponding profiles for polarization fraction p and polarization angle ψ. It is clear that even when averaged over a large area, the Q and U cumulative profiles do not converge to zero. The corresponding values of the leakage are summarized in Table 4. Before leakage correction, the leakage polarization is of the order of p=2.4% when integrated over the whole PSF and p=4.7% when integrated over the PSF down to its FWHM. This is non-negligible compared to the typical polarization of thermal dust in the ISM, with a most likely value of around 3%, as measured by Planck at 353 GHz ( [11]). Given the amplitude of this leak- age effect compared to the expected astrophysical signal, it clearly needs to be subtracted accurately from the data. This effect is known as instrumental polarization and is also often referred to as leakage from total intensity to polarization. It has been observed in the polarization data of many, if not all, instruments measuring polarization. Some procedures have been proposed to subtract the effect from the data, for instruments such as NIKA2 ( [START_REF] Ritacco | [END_REF]), ACTPol ( [40]) or HAWC+ ( [30]). The origin of the effect is unclear, but observations suggest that it is due to the propagation of light in the instrument. As seen in Fig. 1, the leakage pattern of the PILOT instrument does not show any distinctive structure. This is unlike the pattern observed for the NIKA2 instrument for instance. This difference may be due to the fact that we use an off-axis telescope that has no occultation by the support for the secondary mirror. The leakage is also seen to produce mostly U in instrument coordinates, which indicates that the instrumental polarization is mostly horizontal in those coordinates, corresponding to a polarization vector roughly vertical. This is compatible with the leakage being due to asymmetries in the optical system, which is mostly symmetrical with respect to the vertical direction for the PILOT instrument (see for instance [START_REF] Bernard | the PILOT team[END_REF]).

In order to investigate possible variations of the leakage with the position in the focal plane, we constructed separate maps of Jupiter with the 5 operational arrays available during flight#3. For each image, we computed the integrated leakage polarization fraction and angle integrated over the FWHM of the total intensity beam. Figure 6 shows the distribution of the polarization fraction and orientations in the PILOT focal plane. The direction is roughly vertical across the focal plane and the fraction also varies slightly. Note that the directions are consistent between arrays 2 and 8 which are optical conjugates on the sky.

Leakage correction performances

The scheme we use to subtract the leakage is illustrated in Fig. 7. We adopt the description proposed by [START_REF] Ritacco | [END_REF] for the NIKA2 data, in which the leakage can be computed as the convolution of the total intensity map of the sky with the leakage PSF measured on a planet. The method implemented in the PILOT pipeline involves rotating the Q and U leakage PSFs by the parallactic angle to bring them to the correct sky orientation. The rotated leakage PSF maps are then convolved with the total intensity map of the sky to produce some leakage Q and U maps. Note that this has to be done at each time sample of the time-line to account for the continuous sky rotation. The leakage Q and U maps are then interpolated at the sky coordinates corresponding to each data sample, in order to predict a leakage signal through Eq. 1. That signal is then subtracted from the original timeline and the corrected timeline is used to produce a map corrected for the leakage, using the Scanamorphos map-making algorithm. Note that, unlike in [START_REF] Ritacco | [END_REF], we do not subtract any fixed polarization contribution other than the above leakage. The performances quoted therefore reflect the correction of the Fig. 5: Cumulative intensity profiles as a function of integration radius from the peak of the total intensity in Jupiter maps. The black curve shows the cumulative profile of the total intensity. The red and dark-blue curves show the Q and U profiles multiplied by 10. The light blue and yellow curves show the corresponding average polarization leakage fraction in units of 10%, and the corresponding polarization angle divided by 90 o . The top panel shows the initial leakage characteristics before correction. The bottom panel shows the leakage characteristics after correction for readout latency, cross-talk and leakage. The erratic behavior of the polarization angle on this panel is due to the very low magnitude of Q and U after correction. intensity to polarization leakage through the algorithm described here.

For the leakage PSF maps, we use the Jupiter maps shown in Fig. 1 computed in instrument coordinates and with Q and U also in reference to instrumental coordinates. As a consequence, we set the parallactic angle φ to zero in Eq. 1 when computing the leakage contribution to be subtracted from the original timeline. For the total intensity map, we use the Herschel map of the astronomical object extracted from the ESA Herschel Science Archive 1 , which we reproject into the appropriate equatorial coordinates. This approach is pre- ferred to deconvolving our own intensity map from our total intensity PSF, given the accuracy required for a proper subtraction of the leakage signal. In the case of Jupiter, we use a fake source map where the total intensity map of the planet is computed as a Gaussian with FWHM of 18.1 ′′ , i.e. the FWHM of the 250 µm Herschel band. In both cases, we correlate the observed total intensity PILOT map of the object with the Herschel map and use the linear scale factor between the two images to rescale the Herschel map prior to using it to predict the leakage. In practice, the calculation of the leakage maps is performed at discrete parallactic angle values covering the range for each PILOT observation tile, with a discretization step of ≃ 1 • and the de-projected timelines are interpolated at the actual parallactic angle value of each data sample from those maps, using linear interpolation. The above processing is computed independently for each detector array and we use the PSF of each array as measured in maps of Jupiter computed for that array only. The data processing steps proceed clockwise from top-left. The I, Q and U PSF as measured on Jupiter are rotated to follow the time-dependent sky rotation. They are then convolved to the high resolution Herschel satellite sky map of the considered object (here the Orion A region) at 250 µm sky to produce instantaneous maps of the Q and U leakage. The instantaneous maps are re-observed to produce a correction time-line which is subtracted from the original data. Our final leakage-corrected, science-ready maps are constructed by applying our map-making algorithm to this corrected timeline.

To first order, the processing therefore takes into account the focal-plane variations of the leakage shown in Fig. 6.

The residual polarization after leakage correction on the Jupiter data is shown in the bottom panel of Fig. 5. It can be seen that the polarization leakage is strongly reduced compared to the profiles prior to correction shown in the upper panel. The residual polarization leakage as measured at the FWHM radius of the PSF and integrated over the whole PSF are given in tab. 4 and are respectively 0.17% and 0.38%. The spatial distribution of those residuals in the image shown on the radial profiles is likely due to uncertainties in the pointing reconstruction or in the PSF shape, discretization steps used for map rotations, etc. . . We also stress that a high level of accuracy must be preserved at each step of the leakage correction process, e.g. during the map rotations, convolution and timeline de-projection. Most common reprojection routines do not guarantee sufficient accuracy, and we were obliged to use drizzling methods [41] at each step. It also required working with map pixels of 6 ′′ consistent with the Herschel resolution, which are ≃ 22 times smaller than the PILOT beam, to produce acceptable accuracy of the correction. We note that, due to filtering of large scale emission inherent to measurements with bolometers, the leakage on extended sources should be lower than measured here on a point-like source. We consider that the range of performances given in Tab. 4 reflect those attain-able on astrophysical sources with the current leakage subtraction method.

Conclusions

In this paper, we have presented the methods used to correct for residual systematic effects in the PILOT data, in addition to those already described in [1]. In particular, we describe the cross-talk and read-out latency systematic effects that affect the shape of the instrument PSF in total intensity. We also described how we measure and correct for the intensity to polarization leakage and the method we use to subtract it from the data. The cross-talk and read-out latency effects are observed in the total intensity maps of Jupiter obtained during the third flight of the PILOT instrument as distortions of the instrument PSF. We measured the parameters characterizing those effects by using a simple pixel-to-pixel transfer model and derived the transfer coefficients by minimizing the PSF defects in the images of Jupiter. Our analysis showed that the read-out latency effect is observed on the cross-talk signal, indicating that it arises after cross-talk in the detection chain and therefore needs to be corrected first.

Following the above correction, images of Jupiter in polarization show polarized signal at the level of ≃ 3%, which we interpret as leakage from total intensity to polarization, also known as instrumental polarization. Polarization leakage likely affects most instruments measuring polarization in the FIR/submm at a similar level. Using images of Jupiter obtained separately for the five bolometer arrays that were operational during flight#3, we showed that the leakage is mostly oriented parallel to the symmetry axis of the instrument, which points towards residual asymmetries of the optics as the origin of the leakage. We also showed that the polarization direction and fraction vary across the focal plane, an effect that we also take into account in the correction. We correct for the leakage in the PILOT pipeline following the method initially proposed by [START_REF] Ritacco | [END_REF] to correct for the leakage in the NIKA2 instrument. We use the I, Q and U PSFs measured on Jupiter, rotated to follow sky rotation and convolved with a scaled intensity sky map obtained by the Herschel satellite at 250 µm to predict maps of the leakage in sky coordinates. Those convolved leakage maps are deprojected onto the timeline of each detector to infer the correction for each detector, taking into account at first order the observed spatial variations of the leakage. We emphasize that accuracy must be preserved at each step of the process in the map rotations, convolution and timeline de-projection. This is not guaranteed by most reprojection routines. Applying our correction strategy to the Jupiter data and using a simple synthetic PSF model for the planet yields a residual polarization fraction lower than ≃ 0.4%, which we regard as the accuracy of our leakage correction method.

Fig. 1 :

 1 Fig. 1: Total intensity I (left), Q (center) and U (right) beam maps obtained on Jupiter during flight#3 before correction of systematic effects discussed here. The images are shown in instrument coordinates with elevation increasing upward and cross-elevation (azimuth) increasing to the right. The images are shown in arbitrary units in logarithmic scale for I and linear scale for Q and U. The elongation across the first diagonal in the total intensity image is due to cross-talk. The negative shadow of the cross-talk signal, the PSF distortion visible above the lower left-upper right diagonal and the faint residual source appearing below the planet along the other diagonal are due to read-out latency. The non-zero Q and U originate from intensity to polarization leakage. These effects are described in Sect. 4

Fig. 2 :

 2 Fig. 2: From left to right: (a): Schematic view of the array showing the readout direction and the effect of cross-talk along lines and read-out latency effects. (b): Image of Jupiter obtained by leaving in glitch detected samples, overlaid with the footprint of the detectors, showing the effect of the cross-talk along the lines and the presence of a ghost opposite to the planet along readout columns due to read-out latency. (c): The regions overlaid on the image of Jupiter delineate the zones used to measure the read-out latency parameters (red) and the corresponding reference regions (black). (d): The regions overlaid on the image of Jupiter delineate the zones used to measure the cross-talk effect parameters (red) and the corresponding reference regions (black).

Fig. 4 :

 4 Fig. 4: Images of Jupiter constructed from uncorrected data (left), corrected for read-out latency (middle), corrected for read-out latency and cross-talk (right). All images are shown in the same arbitrary units in logarithmic scale.

Fig. 6 :

 6 Fig.6: Distribution of the leakage polarization direction and intensity across the PILOT focal plane. The focal plane is shown in instrument coordinates, with elevation corresponding to the vertical direction, in offsets from the focal plane center. The circles show the positions of individual pixels of the TRANS arrays. The labels give the name of each operational array in red and blue for the REFLEX and TRANS focal plane respectively. The lines show the direction and amplitude of the intensity to the polarization leakage as measured on Jupiter for each array independently, in red and blue for the REFLEX and TRANS arrays respectively. The scale on the bottom left shows 5%.

Fig. 7 :

 7 Fig. 7: Schematic representation of the method used to correct for the intensity to polarization leakage in the PILOT data.The data processing steps proceed clockwise from top-left. The I, Q and U PSF as measured on Jupiter are rotated to follow the time-dependent sky rotation. They are then convolved to the high resolution Herschel satellite sky map of the considered object (here the Orion A region) at 250 µm sky to produce instantaneous maps of the Q and U leakage. The instantaneous maps are re-observed to produce a correction time-line which is subtracted from the original data. Our final leakage-corrected, science-ready maps are constructed by applying our map-making algorithm to this corrected timeline.

Table 1 :

 1 Key optical characteristics of the PILOT instrument.

	Telescope type	Gregorian
	Equivalent focal length [mm]	1790
	Numerical aperture	F/2.5
	FOV [ o ]	1.0 × 0.8
	Ceiling altitude	∼3 hPa
	Pointing reconstruction	translation= 1 ′′ , 1σ
		rotation= 6 ′′ , 1σ
	Gondola mass	∼1100 kg
	Primary mirror (M1)	Off-axis parabolic
	M1 diameter [mm]	930 × 830
	M1 used diameter [mm]	730
	Focal length [mm]	750
	Detector type	multiplexed
		bolometer arrays
	Total number of detectors	2048
	Detectors temperature [mK]	300
	Sampling rate [Hz]	40
	Photometric channels
	λ 0 [ µm]	240
	ν 0 [GHz]	1250
	∆ ν/ν	0.27
	beam FWHM [ ′ ]	1.9
	Minimum Strehl ratio	0.95

Table 2 :

 2 Observations obtained during flight#3.

	Source	Observation Time [mn]	Map size [deg x deg]	Total depth [deg 2 /h]
	Aquila Rift	128.5	7 x 2	6.5
	Crab nebula 100.	1.5 x 1.2	1.1
	Fan	118.5	5 x 3.2	0.8
	Jupiter	33.	2 x 1	3.6
	M31	301.6	4 x 1.8	1.4
	MW L133	101.58	3 x 2.8	5.0
	Orion	140.1	5 x 2.5	5.3
	Tau B211	50.1	2 x 1.8	4.3
	Tau L1506	160	2 x 1.9	1.4
	SkydipM31	20.	32 x 2	n/a
	SkydipPol	33.1	44 x 2	n/a

Table 3 :

 3 Cross-talk parameters derived for each array. Column 1 gives the name of the corresponding focal plane, column 2 gives the array number, column 3 gives the Buffer Unit (BU) number associated to reading the array, column 4 gives the value for the cross-talk coefficient derived.

	Focal Plane Array BU	f CT
	TRANS	2	2	9.0810 -4
	REFLEX	3	3	8.3110 -4
	REFLEX	4	3	9.0810 -4
	REFLEX	7	4	9.6910 -4
	REFLEX	8	4	9.3310 -4

Table 4 :

 4 Leakage correction performance. Columns 1 to 3 indicate which corrections are applied. Column 4 gives the polarization fraction of the leakage as integrated over the whole PSF. Column 5 gives polarization fraction of the leakage as integrated up to the FWHM of the PSF.

	RL CT leakage p f whm [%] p tot [%]
	-	-	-	5.36	2.19
	x	-	-	5.38	2.62
	x	x	-	4.74	2.41
	x	x	x	0.17	0.28

Acknowledgements PILOT is an international project that involves several European institutes. The institutes that have contributed to hardware developments for PILOT are IRAP and CNES in Toulouse (France), IAS in Orsay (France), CEA in Saclay (France), Sapienza University in Rome (Italy), Cardiff University (UK) and the Scientific Support Office at ESTEC (NL). This work was supported by the CNES. It is based on the PILOT data obtained during three flight campaigns operated by CNES, under the agreement between CNES and CNRS / INSU. It benefited from the active participation of CNES, IRAP, IAS, CEA and ESTEC to the flight campaigns. This work was supported by the Programme National "Physique et Chimie du Milieu Interstellaire" (PCMI) of CNRS / INSU with INC/INP co-funded by CEA and CNES. DA acknowledges SC MES RK grant No. AP08855858 and Nazarbayev University Grant Programme No. 110119FD4503. This work was supported by World Premier International Research Center Initiative (WPI), MEXT, Japan.