text:
# The effect of signal digitisation in CMB experiments

## 1 Introduction

Since the successful detection of CMB anisotropy at angular scales $`\theta \stackrel{_>}{_{}}7^{}`$ by the $`COBE`$-DMR experiment (Smoot et al. 1992; Bennet et al.1996; Górski et al. 1996), a large number of detections on smaller angular scales with increased sensitivity by both ground-based and balloon-borne experiments have followed (e.g. Lasenby et al. 1998, De Bernardis & Masi 1998, Bersanelli & Maino & Mennella 2002 for reviews). Recent exciting results from BOOMERanG (De Bernardis et al. 2000), MAXIMA-1 (Balbi et al. 2000), DASI (Pryke et al. 2001) and CBI experiments (Mason et al. 2001) have provided strong support to the inflationary scenario for structure formation and a universe with $`\mathrm{\Omega }_01`$. However to fully extract the wealth of cosmological information encoded in the CMB angular power spectrum, extreme mapping sensitivity is required. Only space missions can provide the proper combination of full-sky view and low level of systematic effects that is necessary for “precision” CMB measurements (Danese et al. 1996), as demonstrated by the excellent results from the NASA satellite MAP 111http://map.gsfc.nasa.gov/ (Bennett et al. 2003). The ESA satellite Planck 222http://planck.esa.nl (Mandolesi et al. 1998; Puget et al. 1998) is designed to measure the CMB anisotropy with an accuracy mainly set by astrophysical limits on a wide frequency range (30–900 GHz) and with unprecedent angular resolution and sensitivity.

While space experiments benefit from unique environmental conditions, numerous technical issues need to be addressed. One of them is the limited bandwidth for downloading of data to the Earth which often imposes strong limits on the available telemetry rate. For this reason it is often necessary to use loss less compression methods (Maris et al. 2000a) in order to reduce the data rate. Because the efficiency compression algorithms depends on the quantization step applied to the observed signal by the instrument acquisition chain on board the satellite, it is necessary to find the optimal trade-off between the required compression rate and the quantisation step, which needs not to be too high to preserve the scientific quality of the measured data.

If the quantization step, $`q`$, is significantly smaller than the signal root mean square (RMS), $`\sigma _\mathrm{T}`$, and no saturation of the quantization levels occurs then the distortion induced by the signal quantization can be well approximated by an additive source of non-Gaussian white noise with variance $`q^2/12`$ (Kollar, 1994). Corrections to this simplified model are required when $`\sigma _\mathrm{T}/q\stackrel{_<}{_{}}1`$, which depend on the kind of the total signal in input to the quantization stage of the on-board acquisition chain (Kollar 1994). For example, in the case of a pure white noise signal they are exponentially damped by a factor $`\mathrm{exp}[2(\pi \sigma _\mathrm{T}/q)^2]`$. On the contrary, in the case of a sinusoidal signal (e.g. like the CMB dipole), corrections may by quite large and no simple analytical expressions can be written.

It is worth noting that, in general, the observing strategy implies repeated observations of the same region of the sky on different time scales with detectors operating at different frequencies. Therefore, considering that a sample is the result of the average of $`N`$ repeated observations, we have that for a normally distributed signal with $`\sigma _\mathrm{T}/q1`$ the RMS of the quantization noise, $`\sigma _q`$ is:

$$\sigma _q\frac{q}{\sqrt{12N}}.$$
(1)
In CMB experiments, however, the approximations of normal distributed signal and low $`\sigma /q`$ are often not strictly applicable. In fact the astrophysical signal is in general not normally distributed, and it is often required to use large quantization steps (i.e. $`\sigma _\mathrm{T}/q1`$). The instrumental noise power spectrum usually shows a $`1/f^\alpha `$ spectral behaviour, with $`\alpha 12`$ (see Burigana et al. 1997a and Seiffert et al. 2002), while the number of samples is not uniformly distributed when projected onto a sky map. Therefore it is apparent that the validity of eq. (1) needs to be assessed in the context of space CMB experiments in order to evaluate possible deviations from this simple formula.

In this paper we evaluate the impact of the quantization error on CMB measurements from space and provide tools and methods to optimize the quantization step to reach the required compression rate with minimal impact on the scientific output.

We show examples relative to the Planck/LFI 30 and 100 GHz channels. However, the formalism, tools and methods developed here are general and applicable to CMB anisotropy experiment based on differential receivers, as well as to other kinds of experiments with large quantization steps.

Using Monte Carlo simulations we estimate the effect of software quantization on the data considering various quantization functions and for different values of the quantization step. We also discuss the impact of different experimental conditions like the total number of repeated measures for sample, $`N_s`$, the white noise RMS, $`\sigma `$, the $`1/f^\alpha `$ noise and the scanning strategy. The impact of quantization is evaluated at three different levels: on time ordered data, on sky maps and their statistical moments, and on the angular power spectrum.

In Sect. 2 we describe the analytical model of signal quantization and discuss the adopted numerical approach. Sect. 3 describes our results in terms of effects on TOD and maps including their statistical properties, and on the angular power spectrum. In Sect. 4 we present the possibility to reduce the effect of quantization in the recovered power spectrum. Finally, we discuss our main results and draw the conclusions in section 5. To improve readability, some mathematical details are reported in Appendix A while in Appendix B we briefly verify that the extra-quantitazion step introduced via hardware by the analog to digital converter (ADC) quantization represents a minor effect with respect to the subsequent software quantization extensively discussed in the paper. In Appendix C we report on statistical moments of the map adopted in the simulations and of several foregrounds, for comparison.

## 2 Analytical model and numerical simulations

### 2.1 Analytical model

During a typical CMB experiment from space, the measured signal is first digitized and compressed on board and subsequently reconstructed on ground. Furthermore the scanning strategy usually involves multiple measurements from the same pointing direction that are subsequently averaged during ground data analysis. Therefore in order to assess the impact of the quantization error we need to consider the complete process, from quantization, to reconstruction and averaging. In the following of this paper this process will be referred to as “QRA process”.

Let us neglect the data compression step (that is out of the scope of our paper) and consider first the on-board quantization. In general the instrument output will be constituted by a Time Ordered Data stream (TOD) representing the sky temperature anisotropy; in the case of differential receivers, for example, the TOD will contain consecutive samples of temperature differences (sky minus reference signal) so that the $`i^{th}`$ sample, $`\mathrm{\Delta }T_i`$, will be quantized as (Maris et al. 2000a ):

$$\mathrm{\Delta }T_i^\mathrm{Q}=𝒬\left[\frac{\mathrm{\Delta }T_i\mathrm{\Delta }T_{}}{q}\right],$$
(2)
where $`𝒬[x]`$ and $`\mathrm{\Delta }T_{}`$ represent the quantization operator and a zero level.

The second step is the data reconstruction, i.e. the inverse of eq. (2). This second step is necessary because both the quantization step $`q`$ and the offset $`\mathrm{\Delta }T_{}`$ will not necessarily be constant during the mission; this implies that eq. (2) needs to be inverted in order to be able to compare measurements obtained in different periods of observation. The data reconstruction step can be defined as:

$$\mathrm{\Delta }T_i^\mathrm{R}=\left(\mathrm{\Delta }T_i^\mathrm{Q}+𝒪_{}\right)q+\mathrm{\Delta }T_{}$$
(3)
where $`𝒪_{}`$ is an offset that depends on the considered quantization operator, and is defined by:

$$𝒪_{}=\mathrm{E}[\mathrm{\Delta }T_i]/q\mathrm{E}[\mathrm{\Delta }T_i^\mathrm{Q}]$$
(4)
where $`\mathrm{E}[X]`$ represents the expectation of the random variable $`X`$.

The value of $`𝒪_{}`$ is dependent on the quantization operator and assumes the values of $`0`$, $`+1/2`$, $`1/2`$ for the $`\mathrm{floor}(x)`$, $`\mathrm{ceil}(x)`$ or $`\mathrm{round}(x)`$ operators, respectively. In the case of the $`\mathrm{trunc}(x)`$ operator the offset depends on the statistics of the random variable $`X`$; if $`X>0`$ then $`\mathrm{trunc}(x)\mathrm{floor}(x)`$, if $`X<0`$, $`\mathrm{trunc}(x)\mathrm{ceil}(x)`$.

In our analytical treatment we will neglect the presence of an offset; its effect will be discussed in the next section. Other simplifying assumptions are that the signal is quantized in equally spaced steps and that $`\mathrm{\Delta }T_{}`$ is always $`>0`$. With simple considerations we can estimate the effect of the QRA process on the reconstructed maps and power spectra. We consider the case of typical CMB anisotropy experiments consisting of collections of many TODs each one based on the averaging over $`N_\mathrm{s}`$ repetitions of the observation of the same stripe in the sky 333For example, Planck TODs consist of about $`1.1\times 10^4`$ scan circles (for about 14 months of observations) each of them observed 60 times. See also section 2.2..

Assuming that the QRA acts mainly as an added white noise source, the effect of quantization on the power spectrum can be easily estimated. Let $`C_{\mathrm{},N_\mathrm{s}}^{wn}`$ the white-noise induced excess of power and $`C_{\mathrm{},N_\mathrm{s}}^{\mathrm{QRA}}`$ the equivalent excess of power induced by the QRA process; then we have that $`C_{\mathrm{},N_\mathrm{s}}^{wn}=K_{\mathrm{wn}}^2\sigma ^2/N_\mathrm{s}`$ and $`C_{\mathrm{},N_\mathrm{s}}^{\mathrm{QRA}}=K_{\mathrm{qra}}^2q^2/(12N_\mathrm{s})`$, where $`N_\mathrm{s}`$ represents the number of times that each pixel in the sky is measured during a single stripe and the pixel dependent normalization constants $`K_{\mathrm{wn}}^2`$ and $`K_{\mathrm{qra}}^2`$ are determined by the scanning strategy (and the processing of the data, assumed to be a linear process). Since both quantized and not-quantized signals are sampled and processed in the same way, in principle we expect $`K_{\mathrm{wn}}^2=K_{\mathrm{qra}}^2`$. Therefore, for any $`N_\mathrm{s}`$ we have:

$$\frac{C_{\mathrm{},N}^{wn}}{C_{\mathrm{},N}^{\mathrm{QRA}}}=\left(\sqrt{12}\frac{\sigma }{q}\right)^2.$$
(5)
The total noise power spectrum after quantization can be well approximated by:

$$C_{\mathrm{},N}^{total}=C_{\mathrm{},N}^{wn}\left[1+\left(\sqrt{12}\frac{\sigma }{q}\right)^2\right].$$
(6)
For $`\sigma /q2`$ (as in the case of the Planck baseline) the noise power excess will be increased by $`2`$%. A similar relation holds for the ratio between the variance of the quantization noise and the variance of the white noise in a map:

$$\sigma _{\mathrm{map},\mathrm{total}}^2=\sigma _{\mathrm{map},\mathrm{wn}}^2\left[1+\left(\sqrt{12}\frac{\sigma }{q}\right)^2\right],$$
(7)
where $`\sigma _{\mathrm{map},\mathrm{wn}}`$ is the RMS white noise in the map introduced by the instrument. This formula holds both for the noise average variance on the map and for the noise variance of a given pixel in the map.

### 2.2 Numerical simulations

The analytical model described in the previous section allows a simple, first-order estimation of the QRA effect. We have performed numerical simulations to verify and quantify the validity and the accuracy of the analytical model and to study effects (like the presence of the skysignal and $`1/f^\alpha `$ noise) not accounted for in the analytical approach. The simulations are representative of the Planck mission. In particular we consider here the 30 GHz and 100 GHz channels of Planck/LFI assuming angular resolutions of $`33^{}`$ and $`10^{}`$, respectively.

#### 2.2.1 Simulation tools

To evaluate the effect of QRA process through numerical simulations we have $`a)`$ simulated the data acquisition during the space mission (with and without signal QRA processing) producing TODs; $`b)`$ produced reconstructed sky maps and power spectra from simulated TODs; $`c)`$ compared the results when including or not the QRA process.

The simulation of the data acquisition during the mission has been performed using the Planck Flight Simulator developed by Burigana et al. (1997b ). In this code we have included numerical tools to simulate data acquisition, signal quantization and signal reconstruction to allow the parallel generation of a quantized data stream plus a data stream of reconstructed data and to perform the statistical evaluation of the quantization error at TOD level (Maris et al. 2000b ).

The Flight Simulator is a code that simulates the CMB measurement according to the Planck satellite scanning strategy. Planck is a spinner (with the telescope observing axis nearly perpendicular to its spin axis) that will map the sky through a large set of near great circles at a spin rate of 1 r.p.m. The satellite will orbit around the Lagrangian point L2 of the Sun-Earth system maintaining the spin-axis always in the anti-Sun direction. The measurement strategy is to acquire data from the same sky circle for about one hour before repointing the spin axis by approximately $`2.5^{}`$, so that the same sky pixel is measured $`60`$ times in each scan circle.

The main simplifying assumptions concerning the scanning strategy that we have considered in our simulations are: (i) no drift or variation in the spin rate (equal to 1 r.p.m. for all the simulation), (ii) instantaneous repointing, (iii) perfect pointing with $`85^{}`$ observing angle with respect to the spin axis.

The only astrophysical signal considered in the simulation has been a map describing the CMB fluctuations, therefore neglecting the effect of the dipole and of the galactic and extragalactic contributions. To the astrophysical signal we have also added white noise representative of the 30 GHz and 100 GHz radiometers (with an RMS of 1.056 mK and 3.325 mK, respectively, corresponding to integration times of 0.03 s and 0.009 s). In some simulations we also added a 1/$`f`$ component to study the effect of correlated noise on quantization; to this aim we considered two cases, with a 1/$`f`$ component having a knee frequency of 0.1 and 0.5 Hz, respectively. The first is representative of a typical worst case expected in LFI radiometers (for which we require a knee frequency $`0.05`$ Hz) while the second has been studied to enhance any possible effect connected to 1/$`f`$ noise but is clearly non representative of expected performances.

The code developed by Burigana et al. (1997b ) and Maino et al. (Maino:destriping (1999)) has been used to generate maps from TODs. With this code it is also possible to significantly reduce correlated systematic effects as the $`1/f^\alpha `$ noise stripes from the data before generating the final maps. Power spectra have been obtained from the generated maps using the anafast code of the HEALPix package from Górsky et al. 1998, adopted here to pixelize the sky.

#### 2.2.2 Numerical procedures

Let us now consider the time ordered data stream produced by the Flight Simulator as an ordered sequence of measured temperature differences $`\mathrm{\Delta }T_i`$, as in eq. (2). If we now take into account that for each pointing direction in the sky, indicated by the subscript $`p`$, on a given stripe (i.e., a scan circle for Planck) there will be $`N_\mathrm{s}`$ repeated measurements then we can reorder the TOD according to two indices, $`(s,p)`$, the first one ($`s=1,\mathrm{},N_\mathrm{s}`$) representing the repeated pixel measurements and the second one representing the pointing direction. Therefore we will have that $`\mathrm{\Delta }T_i\mathrm{\Delta }T_{s,p}`$ ($`i(s,p)`$).

Now for a single sample $`(s,p)`$ the quantization error introduced by a quantization step, $`q`$, can be defined as:

$$\mathrm{\Delta }T_{s,p}^\mathrm{E}=\mathrm{\Delta }T_{s,p}^\mathrm{R}\mathrm{\Delta }T_{s,p},$$
(8)
where the superscript $`R`$ indicates a sample that has been quantized and then reconstructed. In a similar way we can define the following quantities averaged over $`N_\mathrm{s}`$ multiple measurements for each direction, $`p`$, in the sky:

$$\mathrm{\Delta }T_p=\frac{1}{N_\mathrm{s}}\mathrm{\Sigma }_s\mathrm{\Delta }T_{s,p},$$
(9)
$$\mathrm{\Delta }T_p^R=\frac{1}{N_\mathrm{s}}\mathrm{\Sigma }_s\mathrm{\Delta }T_{s,p}^\mathrm{R},$$
(10)
$$\mathrm{\Delta }T_p^E=\mathrm{\Delta }T_p^R\mathrm{\Delta }T_p=\frac{1}{N_\mathrm{s}}\mathrm{\Sigma }_s\mathrm{\Delta }T_{s,p}^\mathrm{E}.$$
(11)
To characterize the quantization error for the LFI signal at TOD level we study the behavior of the distribution moments of $`\mathrm{\Delta }T_p`$, $`\mathrm{\Delta }T_p^R`$ and $`\mathrm{\Delta }T_p^E`$. In particular we analyze their expectation, RMS, skewness and kurtosis as a function of the quantization step $`q`$, of the quantization function $`𝒬[]`$ and of the number of samples $`N_\mathrm{s}`$.

The time ordered data with and without quantization have then be used to produce maps and power spectra in order to study the quantization error at these other two levels.

The simulations have been carried out considering three distinct quantization steps, i.e. $`q0.5`$, 1.0, 2.0 mK for the 30 GHz channel and $`q1.2`$, 2.4, 4.9 mK for the 100 GHz channel. These three values of $`q`$ correspond, for each frequency channel, to $`\sigma /q2`$, 1 and 1/2, respectively. The case $`\sigma /q2`$ corresponds to the current Planck baseline; the case $`\sigma /q1/2`$ is representative of a case in which the noise model for the quantization distortion is expected to fail, while the intermediate case was chosen evaluate possible deviations of the dependence of the quantization error on $`\sigma /q`$ from a power law.

## 3 Results

In this section we discuss the effect of signal quantization at the three different levels of time ordered data, final maps and power spectra. The results presented in this section have been produced considering the quantization operator $`\mathrm{floor}(x)`$; this assumption can be made without any loss of generality, as all the quantization operators can be considered equivalent provided that the proper offset is applied in the reconstruction formula, eq. 3.

### 3.1 Effect of quantization at TOD level

To study the effect of quantization at the level of the time-ordered data stream we focus on the statistical properties of the data-stream by evaluating the various moments of the statistical distribution of the quantization error.

#### 3.1.1 Mean value of the quantization error distribution

The first statistical quantity that we have analyzed is the mean value of the quantization error to evaluate if the QRA process introduces any spurious bias in the data and how this varies with $`q`$. This is relevant in the context of space missions, in which possible variations in the noise properties of the detectors will require slight modifications in time of $`q`$ in order to maintain constant the $`\sigma /q`$ ratio required for optimal compression. Therefore if any bias in the data is introduced by quantization it is reasonable to expect that this bias will vary in time (because of the changes in $`q`$) leading to systematic errors in the final maps and power spectra.

We have calculated the mean value of the quantization error ($`\mathrm{\Delta }T_{60}^\mathrm{E}=\mathrm{mean}[\mathrm{\Delta }T_{60,p}^\mathrm{E}]`$, see eq. 11) for a $``$ 1 year-long data stream, considering different values of the quantization step $`q`$. Our results show that quantization introduces a bias in the data at a level $`|\mathrm{\Delta }T_{60}^\mathrm{E}|/q<10^3`$ which implies $`|\mathrm{\Delta }T_{60}^\mathrm{E}|<2`$ $`\mu `$K even for the largest values of $`q`$. This offset effect is small, and in addition it can be efficiently removed by destriping algorithms.

Therefore we can expect systematic effects in the time ordered data at a very low level and the time scales of the updating of $`q`$ on board ($`>`$ 24 hours). Furthermore these effects can be efficiently reduced by applying destriping and/or map-making methods to the time-ordered data and therefore they do not represent a concern444In fact, the current destriping codes are able to significantly reduce the impact of much larger drifts, such as those induced by the $`1/f^\alpha `$ noise or by periodic fluctuations induced by thermal fluctuations..

#### 3.1.2 RMS of the quantization error distribution

The rms of the quatization error can be easily evaluated from the theoretical eq. (1). The validity of this relationship has been checked numerically by calculating the rms of the quantization error (indicated by $`S_N`$) for several TODs generated with various values of $`q`$ and $`N`$ around the reference values $`N_0=60`$ and $`q_0=0.305`$ mK/adu, equivalent to $`\sigma /q=3.462`$ for the 30 GHz channel and $`\sigma /q=10.902`$ for the 100 GHz channel. A fit of $`S_N`$ with the relation:

$$\mathrm{ln}S_N=\mathrm{ln}\mathrm{\Sigma }_0+m_\mathrm{N}\mathrm{ln}(N_\mathrm{s}/N_0)+m_\mathrm{q}\mathrm{ln}(q/q_0).$$
(12)
yielded $`\mathrm{\Sigma }_0q_0/\sqrt{12N_0}`$, $`m_\mathrm{N}1/2`$ and $`m_\mathrm{q}1`$ with an accuracy up to the third decimal figure. This confirms that in our case eq. (1) is a very good approximation of the rms value of the quantization error distribution.

#### 3.1.3 Higher moments of the quantization error distribution

In the ideal case of a normally distributed random signal with $`\sigma /q1`$ the QRA error will also be nearly normally distributed and the indices of skewness and kurtosis 555In what follows we will use sometimes the terms “skewness” and “kurtosis” instead of skewness and kurtosis index. These indices, combinations of the distribution statistical moments $`\mu _i`$, are defined by $`\mathrm{Skewness}\mathrm{index}=\mu _3/\mu _2^{3/2}`$ snd $`\mathrm{Kurtosis}\mathrm{index}=\mu _4/\mu _2^23`$. are expected to vanish. On the other hand if $`\sigma /q`$ approaches 1 then it is reasonable to expect a change in the values of the higher moments of the statistical distribution for the QRA error.

From theory, the effect on higher moments may be evaluated from the characteristic function for the distribution of the sampled values after QRA processing $`\mathrm{\Theta }\mathrm{\Pi }_{N_\mathrm{s}}^{\mathrm{QRA}}(\omega )=\mathrm{\Phi }_{N_\mathrm{s}}(\omega )\mathrm{\Pi }_{N_\mathrm{s}}(\omega )`$, where $`\omega `$ is a real variable, $`\mathrm{\Phi }_{N_\mathrm{s}}(\omega )`$ is the characteristic function for the QRA process (see appendix A), $`\mathrm{\Pi }_{N_\mathrm{s}}(\omega )`$ is the corresponding function for the averaged input data. The expectation for central moments of order $`k0`$ being $`\mu _k=(i)^kd^k\mathrm{\Theta }\mathrm{\Pi }_{N_\mathrm{s}}^{\mathrm{QRA}}(\omega )/d\omega ^k|_{\omega =0}`$. Since the QRA error is even distributed around zero expectations for its moments for $`k`$ odd are null. If the same is true for the distribution of the unprocessed data then the expectations for all the $`\mu _k`$ with $`k`$ odd is null. The same does not hold for the expectation of moments with $`k`$ even, in this case the QRA process may introduce a not null bias for some or all of them. In addition the QRA process randomly scatters about their expectation both odd and even moments. Then QRA process randomly changes the measured value of moments of any order $`k`$ for a given realization respect to the values which have been measured in the case of no-QRA. Higher order moments, at least up to the fourth order, are needed to full-characterize the statistical properties of the QRA error. However the complete study of these features is beyond the scope of this work and we limit our investigation to few cases illustrative of the methods which may be applied to study them and the kind of perturbation which may be expected on real data. We may attempt a deeper investigation of this effect as numerical tools will improve enough to allow the generation and the analysis of a large data set of simulated missions.

To evaluate the departure of the averaged quantization error from a normal distribution we have calculated the skewness and the kurtosis for the distribution of $`\mathrm{\Delta }T_p`$, $`\mathrm{\Delta }T_p^R`$ and $`\mathrm{\Delta }T_p^E`$ for various values of $`q`$ and $`N_\mathrm{s}`$.

As expected, the QRA process alters the sampling skewness in the input signal, so that it has meaning to put:

$$\frac{\mathrm{\Delta }\mathrm{Skewness}}{\mathrm{Skewness}\left[\mathrm{\Delta }T_p\right]}=\frac{\left(\mathrm{Skewness}\left[\mathrm{\Delta }T_p^R\right]\mathrm{Skewness}\left[\mathrm{\Delta }T_p\right]\right)}{\mathrm{Skewness}\left[\mathrm{\Delta }T_p\right]}0.$$

In our simulation we observe values of $`|\mathrm{\Delta }\mathrm{Skewness}/\mathrm{Skewness}[\mathrm{\Delta }T_p]|`$ up to $`50\%`$ at the TOD level, comparable to the random fluctuations induced by the sampling variance. It is worth noting that this value depends on $`N_s`$ so that the variation in skewness decreases with increasing the number of samples of the same data point and that on the final map (in which each pixel will be the result of averaging a large number of measurement, especially close to the ecliptic poles) the impact on the skewness of the CMB distribution is even smaller (see Table 2) and varying with the colatitude.

In Fig. 1 we show the behavior of the kurtosis for the QRA error ($`\mathrm{E}[\mathrm{Kurtosis}[\mathrm{\Delta }T_p^\mathrm{E}]]`$) as a function of $`N_\mathrm{s}`$. From the plot it is apparent that the expectation of the kurtosis goes to zero for increasing values of $`N_\mathrm{s}`$. For an input signal with small expected kurtosis, the value of the sampling kurtosis at TOD level can be changed by the QRA process by a factor up to $`2`$, again comparable with the sampling variance for the kurtosis in our TODs. As for the skewness the impact on the kurtosis on maps is even smaller. However it is important to recall that, while skewness is changed randomly by the QLA, i.e. no bias is added by QRA to the skewness estimator, the same does not hold for the kurtosis. For the kurtosis, the bias is $`3/[1+12(\sigma /q)^2]^2`$, equivalent to $`0.001`$ for $`\sigma /q2`$. However, once the quantization step and the scanning strategy is defined, it is a simple matter to calculate this bias and to remove it.

A non-zero value of the kurtosis in the distribution of the QRA error changes the statistical significance of the confidence limits respect to the usual definition of “standard error” for which “$`1\sigma `$ corresponds to a 68.27% confidence level”. To recover the standard error definition we have scaled the variance of the quantization error by a multiplicative factor $`\alpha _{N_\mathrm{s}}\stackrel{_>}{_{}}1`$, which decreases monotonically for increasing $`N_\mathrm{s}`$, whose $`lim_{N_\mathrm{s}\mathrm{}}\alpha _{N_\mathrm{s}}=1`$, and $`\alpha _1=1.18`$. Appendix A describes the computation of $`\alpha _{N_\mathrm{s}}`$.

These results allow to safely apply the usual error propagation to combine both signal plus noise and the QRA error with an accuracy better than some $`\mu `$K. The $`1\sigma `$ confidence level for a given temperature average, $`\mu _{N_\mathrm{s}}`$ may be expressed as:

$$\mu _{N_\mathrm{s}}\pm \sigma _{N_\mathrm{s}}\pm \alpha _{N_\mathrm{s}}\sigma _{q,N_\mathrm{s}}\mu _{N_\mathrm{s}}\pm \sqrt{\sigma _{N_\mathrm{s}}^2+\alpha _{N_\mathrm{s}}^2\sigma _{q,N_\mathrm{s}}^2}.$$
(13)
The first term down the square root represents the effect of the white-noise, while the second term the effect of the QRA error. Note that since the $`\alpha _N`$ correction equations (6) and (7) are valid only for $`N_\mathrm{s}`$ sufficiently large to have $`|\alpha _{N_\mathrm{s}}1|1`$. Otherwise the $`\left(\sqrt{12}\sigma /q\right)^2`$ has to be scaled by $`\alpha _{N_\mathrm{s}}^2`$ in eq. (6) and eq. (7). In the case of Planck-LFI the measurement redundancy ($`N_\mathrm{s}60`$) is such that the distributions of the QRA error is sufficiently well approximated by a normal distribution to apply the usual formula for the error propagation.

### 3.2 Effect on the reconstructed map

In Fig. 2 (top panel) we show a map of the QRA error for $`\sigma /q2`$ for a 30 GHz detector. The map is pixelised according to the HEALPix scheme (Gòrski et al. 1998) with a pixel-size of about 13.7 arcmin corresponding to $`N_{\mathrm{side}}`$ ($`12\times N_{\mathrm{side}}^2`$ is the number of pixels in the map). To generate the map we first produced two TODs of the same simulation (including CMB, white and 1/$`f`$ noise) one with quantization and the other without. The error map was then obtained by differencing the maps generated by the two data streams. A simple “visual” inspection of the map indicates the absence of evident correlate structures666The darker narrow band visible in the map is not due to quantization but to the overlapping of a subset of scan circles. In this region the integration time is greater and, therefore, the noise level is smaller. at few $`\mu `$K level. Although the peak-to-peak variation of the quantization error on the map may appear large ($``$50 $`\mu `$K) we must underline that quantization does not act as a correlated systematic effect, but is mainly an added white noise at a few % level over the detector white noise (see eq. (7)).

In the bottom panel of Fig. 2 we show the ratio between the QRA error distribution and a normal distribution with the same mean and RMS, indicating that the QRA error is not completely normal distributed, with a number of error samples larger than $`2\sigma `$ that is less compared to a normally distributed noise. These “non-Gaussian samples”, however, are mostly contained in a symmetric band at ecliptic latitudes between $`\pm 15^{}`$, and owing to the larger number of samples, this feature of the histogram is less large for the 100 GHz channel than the 30 GHz.

In Tab. 1 we report the % noise increase determined by quantization for various values of $`\sigma /q`$. The comparison between the values determined by simulations with those calculated by eq. (7) shows the departure from the simple theoretical model described in Sect. 2.1. For $`\sigma /q2`$ eq. (7) still represents a good estimate of the main effect of quantization.

There is an interesting parallelism between what happens to moments for a TOD of pure white noise when the QRA process is introduced, and what happens for the same moments in the map obtained from the same TOD. To illustrate it, let be to denote with $`\mu _{k,map,wn,QRA}`$ and $`\mu _{k,TOD,wn,QRA}`$ the $`k`$-th moment for the map obtained from QRA processed data and its corresponding TOD from which the map is obtained, and with $`\mu _{k,map,wn}`$ and $`\mu _{k,TOD,wn}`$ the corresponding moments obtained for the map without quantization and the corresponding TOD. Looking to pixels as averages of $`N_p`$ samples, where $`N_p1`$ is the number of TOD samples entering a given pixel $`p`$ in the map, it is possible to classify pixels as a function of $`N_p`$. If $`[N_p]`$ is the class of pixels averages of $`N_p`$ samples, $`g(N_p)`$ the partition function for the pixels over the various classes it is easy to demonstrate that, for any value of $`N_\mathrm{s}`$ and $`\sigma /q`$, and $`k=2`$ the relation

$$\frac{\mu _{k,\mathrm{map},wn,\mathrm{QRA}}}{\mu _{k,\mathrm{map},wn}}=\frac{\mu _{k,\mathrm{TOD},wn,\mathrm{QRA}}}{\mu _{k,\mathrm{TOD},wn}}$$
(14)
holds exactly. From this eq. (7) may be easily derived. For higher order moments ($`k>2`$) the relation in eq. (14) holds better and better as $`N_\mathrm{s}`$ is larger and larger. Indeed, if $`N_\mathrm{s}`$ is not large enough other terms including $`g(N_p)`$ will affect the right hand side of eq. (7). In the case of Planck having $`N_\mathrm{s}=60`$ these terms are negligible at least up to $`k=4`$. In conclusion, from the same reasoning leading to eq. (7) it is possible to draw equations linking the variations of higher order moments on TODs to the corresponding variations results of higher order moments on maps.

As anticipated in the section related to the TODs, we compared the kurtosis in a map which has been QRA processed ($`\mathrm{Kurtosis}[\mathrm{QRA}]`$), against the kurtosis in the corresponding original map ($`\mathrm{Kurtosis}[\mathrm{Orig}]`$) giving the ratio of the variation ($`|\mathrm{\Delta }K/K|=|\mathrm{Kurtosis}[\mathrm{QRA}]\mathrm{Kurtosis}[\mathrm{Orig}]|/\mathrm{Kurtosis}[\mathrm{Orig}]`$). In general $`|\mathrm{\Delta }K|/K`$ is greater for the 30 GHz channel than for the 100 GHz. In particular for the 30 GHz channel considering maps which does not contain $`1/f`$ noise and consequently have not been destriped the kurtosis may be varied even of a 30% by the QRA process. However, comparing maps which does not contain $`1/f`$ noise and have not been destriped with maps originated from data containing $`1/f`$ noise and destriped it is possible to see that the relative effect on the kurtosis is much smaller than in the previous case. This is nothing else than the fact that $`1/f`$ noise and destriping increases the kurtosis of the unquantized signal given at the denominator of the $`|\mathrm{\Delta }K/K|`$ ratio, but leaves practically unchanged the $`\mathrm{\Delta }K`$ difference. From Tab. 1) it is also possible to see how $`|\mathrm{\Delta }K/K|`$ scales with $`\sigma /q`$ both for the 30 GHz and 100 GHz maps. Comparing maps which contains even the $`1/f`$ noise and have been destriped (last column of $`|\mathrm{\Delta }K/K|\stackrel{_<}{_{}}4.2\%`$ ($`1.5\%`$). For both 30 GHz and 100 GHz maps $`|\mathrm{\Delta }K/K|`$ increases up to a factor of about 7 decreasing $`\sigma /q`$ from $`2`$ down to $`1/2`$.

The square root of the sampling variance of $`K`$ is $`\sqrt{24/N_{\mathrm{pix}}}`$ (see e.g. Kendal & Stuart 1977), i.e. $`\sqrt{2}/N_{\mathrm{side}}`$ in the case of the HEALPix scheme, over a map with $`N_{\mathrm{pix}}`$ pixels obtained from a single receiver (radiometer for LFI). For maps at $`N_{\mathrm{side}}=256`$ (512) this gives $`5.5\times 10^3`$ ($`2.8\times 10^3`$), which is equivalent to a sensitivity $`|\mathrm{\Delta }K/K|3050\%`$ (20%) over a map from a single receiver in the case of a sky of pure CMB fluctuations (see Tab. C.1 in Appendix C). When $`|\mathrm{\Delta }K/K|`$ is scaled to take into account that the final map is derived from a combination of more radiometers, the sensitivity to $`K`$ improves to 0.0028 (0.00057) at 30 GHz (100 GHz), equivalent to $`|\mathrm{\Delta }K/K|1525\%`$ (4%). The $`|\mathrm{\Delta }K/K|`$ values induced by the QRA process (see Tab. 2), although not large, is not negligible in the case of maps of pure CMB anisotropy plus noise, particularly because it is not a statistical but a systematic effect. At 30 GHz, the single receiver noise is small, and the above conclusion does not change significantly including also noise contribution to the kurtosis, while at 100 GHz, where a single receiver noise is quite large for LFI, the reference value of $`K`$ changes significantly including or not the noise, and the comparison between $`|\mathrm{\Delta }K/K|`$ representing the sensitivity and that reprenting the QRA effect may depend strongly on the considered realization.

On the other hand, the microwave sky includes fluctuations also from Galactic and extragalactic foregrounds, with remarkable non Gaussian distributions. In Appendix C we report the statistical moments of the most relevant microwave components and of composite maps. As evident, even including Galactic cuts and removing bright sources and Sunyaev-Zeldovich effects (see last three lines of Tab. C.3), foregrounds change the sky kurtosis index with respect to the case of a pure CMB anisotropy sky much more than the QRA process. Therefore, the foreground removal is an aspect much more important than QRA process for the estimation of CMB anisotropy high order moments.

Finally in Fig. 3 we show the variation with the ecliptic latitude of the variance of the QRA error, which decreases as a function of $`\mathrm{sin}\theta `$, where $`\theta `$ is the pixel colatitude. The figure also shows that for $`\theta >15^{}`$ the variance can be well approximated with $`a+b\mathrm{sin}\theta /\mathrm{sin}50^{}`$, where the best-fit coefficients $`a`$ and $`b`$ are listed in the last two columns of Table 1. This relationship represents a useful method to estimate the rms QRA error at a given latitude.

### 3.3 Effect on the power spectrum

In Fig. 4 we show an example relative to a Planck-LFI 30 GHz radiometer that highlights the qualitative features of the quantization effect on the measured power spectrum. To enhance the effect we have chosen a value $`\sigma /q0.5`$, which is clearly not representative of typical conditions.

The figure shows three power spectra. The first one (lower curve) is the power spectrum of the original map. The second (middle curve) is the power spectrum of the map “observed” by the simulator without quantization; the power increase at high multipoles is due to the instrumental white noise. The third (upper curve) is the power spectrum of the “observed” map including quantization. The inset is an enlargement of the power spectrum for $`400\mathrm{}767`$.

From the inset it is evident that the QRA process introduces a power excess in the spectrum. A closer look reveals that most of the features of the original power spectrum are preserved. At multipoles larger than a critical value 777$`\mathrm{}_{\mathrm{crit}}70`$ for the 100 GHz channel and $`\mathrm{}_{\mathrm{crit}}30`$ for the 30 GHz channel. ($`\mathrm{}>\mathrm{}_{\mathrm{crit}})`$ this power excess approaches a value, $`C_{h\mathrm{},60}^{\mathrm{QRA}}`$, nearly independent on $`\mathrm{}`$ and dependent on $`N_\mathrm{s}`$ (we use here $`N_\mathrm{s}=60`$). The value $`C_{h\mathrm{},60}^{\mathrm{QRA}}`$ can be used as a convenient estimator of the quantization effect on the power spectrum at large multipoles and can be quite well estimated by the average of the power excess for $`\mathrm{}>\mathrm{}_{\mathrm{crit}}`$:

$$C_{h\mathrm{},60}^{\mathrm{QRA}}=\frac{_{\mathrm{}=\mathrm{}_{\mathrm{crit}}}^{\mathrm{max}(\mathrm{})}C_{\mathrm{},60}^{\mathrm{QRA}}}{\mathrm{max}(\mathrm{})\mathrm{}_{\mathrm{crit}}+1}$$
(15)
where $`\mathrm{max}(\mathrm{})750`$ and 1500 respectively for the 30 and 100 GHz frequency channels.

The power spectra of absolute and relative quantization error maps, in terms of both $`C_{\mathrm{}}`$ and $`\delta T_{\mathrm{}}`$ are shown in Fig. 5 for different values of $`\sigma /q`$. The power excess varies roughly as a power law of $`\sigma /q`$, in good agreement with eq. (6).

Table 3 gives $`C_{h\mathrm{},60}^{\mathrm{QRA}}`$ for three quantization steps for a single radiometer at the 30 and 100 GHz frequency channels. The simulations from which $`C_{h\mathrm{},60}^{\mathrm{QRA}}`$ are derived assuming the nominal PLANCK same scanning strategy and $`1/f`$ noise knee frequency. In addition we report two cases obtained changing the knee frequency and the scanning strategy, which will be discussed at the end of this section.

A fit of the data obtained by simulations with eq. (6) yielded a residual log-log correlation coefficient $`|1r_{\mathrm{log}\mathrm{log}}|<10^5`$, which indicates an excellent agreement with the expected $`q^2`$ dependence. This allows to parametrize 888In this work we have not studied numerically (because of computation time limits) the dependence of $`C_{\mathrm{}}^{\mathrm{QRA}}`$ on $`N_\mathrm{s}`$. So the $`1/N_\mathrm{s}`$ dependence of eq. (5) $`C_{\mathrm{}}^{\mathrm{QRA}}`$ as a function of $`\sigma `$, $`\sigma /q`$, $`N_\mathrm{s}`$ has been assumed.:

$$C_{\mathrm{}}^{\mathrm{QRA}}=C_{\mathrm{},1}^{\mathrm{QRA}}\frac{\sigma ^2}{N_\mathrm{s}}\left(\frac{\sigma }{q}\right)^2.$$
(16)
The normalization factor $`C_{\mathrm{},1}^{\mathrm{QRA}}`$ can be easily determined from the above-mentioned best-fit; a good approximation can be obtained by replacing $`C_{\mathrm{},1}^{\mathrm{QRA}}`$ with $`K_f/N_{\mathrm{pix}}`$, where $`N_{\mathrm{pix}}`$ represents the number of pixels on the map and $`K_f`$ is a normalization factor $`0.72<K_f<1`$ that depends on details such as: the scanning strategy, the pixelization scheme and the geometry of the instrument. A further approximation $`K_f=1`$ yields:

$$C_{\mathrm{}}^{\mathrm{QRA}}\frac{\sigma ^2}{12N_\mathrm{s}N_{\mathrm{pix}}}\left(\frac{\sigma }{q}\right)^2,$$
(17)
which is a useful relationship to estimate upper limits to the quantization effect on the recovered power spectrum.

Finally we can determine the ratio $`K_{wn}^2/K_{\mathrm{qra}}^2`$ defined in eq. (5) by fitting $`\mathrm{log}C_{\mathrm{},60}^{\mathrm{QRA}}`$ versus $`\mathrm{log}(\sqrt{12}\sigma /q)^2`$ and taking the ratio with the power spectra on the unquantized noise. This yields $`K_{wn}^2/K_{\mathrm{qra}}^2=1.024`$, that represents a good approximation of the ideal value $`K_{wn}^2/K_{\mathrm{qra}}^2=1`$. This is a further confirmation that eq. (5) and eq. (6) well describe the quantization effect on the power spectrum. From these simulations the total power excess due to white noise plus the QRA processing is:

$$C_{\mathrm{}}^{\mathrm{exc}}=C_{\mathrm{},\mathrm{WN}}\left[1+\left(3.55\frac{\sigma }{q}\right)^2\right].$$
(18)
On the other hand, as evident in Fig. 4 referring to realistic simulations including the $`1/f`$ noise and the application of the destriping algorithm, for $`\mathrm{}<\mathrm{}_{\mathrm{crit}}`$, the power excess is no longer constant, but shows an increase for decreasing $`\mathrm{}`$. However, for $`\sigma /q2`$ the deviation from a constant power excess at low $`\mathrm{}`$ is not large. For example, for the 30 GHz channel it is equivalent to $`\stackrel{_<}{_{}}16\%`$ increase of the power excess (50% for the 100 GHz). For $`\sigma /q1`$ the deviation is significant up to $`7.5C_{\mathrm{}}^{\mathrm{QRA}}`$ at 30 GHz and $`16.8C_{\mathrm{}}^{\mathrm{QRA}}`$ at 100 GHz. Although It is not easy to parametrize this effect, $`C_{\mathrm{}}^{\mathrm{QRA}}`$ may be well described by a second or third order polynomial of $`\mathrm{log}\mathrm{}`$ in the log-log space (Maris 2002a ) for $`\mathrm{}>1`$, which for $`\mathrm{}>70`$ well fits the constant power excess previously assumed. However the results of this fit can not be easily generalized as in the previous case, since details such as the polynomial degree are related not only to the mission parameters but also to the assumptions implicit into the destriping procedure and the parameters of the $`1/f`$ noise. This is supported by the comparison of simulations obtained for different values of $`f_{\mathrm{knee}}`$ (see Tab. 4 for the 30 GHz frequency channel). No significant differences appear for $`\mathrm{}>30`$ by changing $`f_{\mathrm{knee}}`$from 0.1 Hz to 0.5 Hz, while at $`\mathrm{}<30`$ the power excess changes up to a factor of two.

A change in the scanning strategy instead, affects the power excess at both low and large $`\mathrm{}`$. In Tab. 4 is shown what happens allowing for a precession of the spin axis during the mission. $`C_{\mathrm{}}^{\mathrm{QRA}}`$ increases by about 7.6% at high $`\mathrm{}`$, and as much as 82% at low $`\mathrm{}`$. Since the same sequence of random noise samples has been used for both the precessed and un-precessed scanning strategy, these changes are due to the different weighting of the samples in the sky induced by the precession of the spin axis.

## 4 Removing the effect of the QRA error in the power spectrum estimation

We have seen that quantization effects on the power spectrum are significantly smaller than the noise and that they are well represented by few parameters. Here we discuss the possibility to remove quantisation excess noise in the data analysis from the final estimated power spectrum. Of course, this is of particular interest at middle and large multipoles while at low multipoles the noise power spectrum is much smaller. The estimation of the sky power spectrum can be derived from the power spectrum obtained from the observed map by subtracting the sum of the expected noise power spectrum and of the expected QRA power spectrum:

$$C_l^{\mathrm{sky}}=C_l^{\mathrm{Obs}}\mathrm{E}[C_{\mathrm{}}^{noise}]\mathrm{E}[C_{\mathrm{}}^{\mathrm{QRA}}].$$
(19)
The accuracy of the subtraction of the QRA effect is determined by the accuracy of the model for $`C_{\mathrm{}}^{\mathrm{QRA}}`$and by the dispersion of the true values of $`C_{\mathrm{}}^{\mathrm{QRA}}`$about their expectation. The typical accuracy of this subtraction at large $`\mathrm{}`$, as evaluated from the last column of Tab. 3, is of about 5% (6%) for the 100 (30) GHz channel, leaving a residual (i.e. after subtraction) QRA error equivalent to an unsubtracted QRA error with a quantization step at least four times smaller than that really used.

The RMS of the expected QRA power spectrum decreases with $`\mathrm{}`$, analogously to the case of the RMS of white noise that exibits a typical $`1/\sqrt{\mathrm{}}`$ dependence. We expect then that the accuracy of the subtraction of the QRA effect improves with increasing $`\mathrm{}`$. For example, tests with the simulated power spectra reported in section 3 show that the RMS of the QRA power spectrum ranges between $`10^{11}`$ mK2 and $`\mathrm{few}\times 10^{10}`$ mK2 at $`\mathrm{}1500`$ (750) for the 100 (30) GHz channel when $`\sigma /q`$ ranges between 2 and 1/2. As a consequence, the power excess is reduced up to a factor 45 (34) for the 100 (30) GHz channel at very large $`\mathrm{}`$, equivalent to an unsubtracted QRA error with a quantization step reduced by a factor $`67`$. On the contrary, at $`\mathrm{}`$ of few tens the RMS of the expected QRA power spectrum is of the same order of magnitude of the power spectrum itself, $`C_{l,60}^{\mathrm{exc}}`$. For example at $`\mathrm{}15`$ it ranges between $`\mathrm{some}\times 10^{10}`$ mK2 and $`\mathrm{some}\times 10^8`$ mK2 for the 100 GHz channel, and between $`\mathrm{few}\times 10^{10}`$ mK2 and $`\mathrm{few}\times 10^8`$ mK2 for the 30 GHz channel. Compared with the averaged power excess for $`\mathrm{}\stackrel{_<}{_{}}30`$ this is equivalent to a reduction of a power excess between $`50`$% and a factor of two.

We searched for possible improvements of the removal procedure by with a polynomial fit of $`\mathrm{log}C_{l,60}^{\mathrm{exc}}`$ versus $`\mathrm{log}\mathrm{}`$. We find that the best fitting relation depends on the frequency channel, the value $`\sigma /q`$, the scanning strategy and the beam position so that a proper fit for each receiver and scanning strategy would be required. However, our simulations show that for $`\sigma /q2`$ the values of $`C_{l,60}^{\mathrm{exc}}`$ are very dispersed around the mean for $`\mathrm{}<30`$, while for $`\mathrm{}>70`$ the goodness of fit does not depend significantly on the degree of the polynomial. Therefore a simple linear fit is adequate.

In order to evaluate the final impact of the accuracy of the subtraction of the QRA quantization effect in the final estimation of the sky power spectrum it is crucial to compare the RMS of the power spectrum of QRA quantization error with the RMS of the noise power spectrum that sets the fundamental uncertainty in the sky power spectrum recovery.

We require that

$$\mathrm{RMS}[C_{l,60}^{\mathrm{exc}}]<\mathrm{RMS}[C_l^{\mathrm{noise}}]/\lambda ,$$
(20)
where the constant $`\lambda `$ has to be set large enough to reduce the residual QRA quantization error to an acceptable level (for example, $`\lambda >5`$ implies a residual QRA quantization error less than 20% of the unsubtractable instrumental noise).

After having applied the destriping algorithm, the instrumental noise is essentially white noise dominated and we can then assume

$$\mathrm{RMS}[C_l^{\mathrm{noise}}]\mathrm{RMS}[C_l^{\mathrm{WN}}]\mathrm{E}[C_l^{\mathrm{WN}}]/\sqrt{\mathrm{}}.$$
(21)
As discussed before, $`C_{l,60}^{\mathrm{exc}}`$ depends on the considered receiver and scanning strategy; of course, the same holds for $`\mathrm{RMS}[C_{l,60}^{\mathrm{exc}}]`$. It is then impossible to rewrite in general the condition (20) in terms of $`\sigma /q`$.

To circumvent this problem we have searched for a law, $`U_{\mathrm{}}`$, which gives at each $`\mathrm{}`$ the maximum among the ratios $`\mathrm{RMS}[C_{l,60}^{\mathrm{exc}}]/C_{l,60}^{\mathrm{exc}}`$ found for the set of simulations considered here. We find that the approximate law

$$U_{\mathrm{}}0.93\mathrm{}^{0.57}$$
(22)
is a quite good approximation for the whole range of multipoles relevant here.

The condition expressed by eq. (20) is then satisfied provided that

$$0.93\mathrm{}^{0.57}\mathrm{E}[C_{l,60}^{\mathrm{exc}}]\stackrel{_<}{_{}}\mathrm{E}[C_l^{\mathrm{WN}}]/\sqrt{\mathrm{}}/\lambda .$$
(23)
At middle and large multipoles the ratio between the variances of the white noise and of the quantization error equals that between the two corresponding power spectra. By using eq. (7), the condition (23) can be then easily rewritten in terms of the $`\sigma /q`$ ratio as:

$$\left(\frac{\sigma }{q}\right)^2\stackrel{_>}{_{}}\frac{0.93}{12}\mathrm{}^{0.07}\lambda .$$
(24)
At $`\mathrm{}\stackrel{_>}{_{}}30`$ (500) this condition is clearly satisfied even for $`\lambda \stackrel{_>}{_{}}4`$, 16 or 65 (5, 20, 80) respectively for $`\sigma /q1/2`$, 1 or 2.

In conclusion, this analysis demonstrates that it is possible to subtract the quantization error impact in the sky power spectrum estimation from Planck data at a level of accuracy which makes its residual effect very small and significantly smaller that the intrinsic unsubtractable uncertainty due to the noise.

However it should be noted that these results are obtained with the assumption of stationary noise. In particular the application of this correction requires to consider every changes in the detector calibration which affects $`q`$ along the mission. Moreover, the previous result about the sensitivity of $`C_{\mathrm{}}^{\mathrm{QRA}}`$ on the scanning strategy imposes to consider that the final $`C_{\mathrm{}}^{\mathrm{QRA}}`$ to be subtracted from the power spectrum shall be evaluated by numerical simulations to take into account the real scanning strategy. This is particularly relevant for a multi-horn instrument like Planck, where a full evaluation of $`C_{\mathrm{}}^{\mathrm{QRA}}`$ shall be obtained only with a full simulation of the mission with each horn observing somewhat differently the sky.

## 5 Conclusions

An unprecedented amount of high quality data are expected from the new generation of CMB experiments designed for extended, high resolution imaging of the microwave sky. These data sets will be analyzed to measure the angular power spectrum to high precision, and also to search for low amplitude statistical properties of the CMB distribution, such as non-gaussianity or high-order autocorrelation functions. It is therefore of interest to accurately consider the effect of instrumental systematics that could significantly perturb the statistical properties of the noise and signal distribution both in the map and in the angular power spectrum recovery.

For a space mission, the need of minimizing the telemetry rate requires that the signal digitization in terms of the $`\sigma /q`$ ratio be as low as possible. On the other hand, a poor signal sampling, combined with the adopted scanning strategy, may produce subtle effects impacting the statistics. In this paper we have analyzed the impact of signal digitization, reconstruction and averaging of a full-sky CMB survey. The effect of the QRA process has been analyzed at the level of data streams, maps and power spectra, both by analytical approximations and numerical simulations. We have simulated the effect for the case of the Planck/LFI experiment, but the methods presented here and the main conclusions can be applied to any other experiment. The simulated signal includes white and $`1/f`$-type noise and a realization of CMB anisotropy, while the effect of QRA process on high order statistics has been compared with that from main Galactic and extragalactic foregrounds.

Because the signal is dominated by white noise, the $`\sigma /q`$ ratio is the main parameter describing the QRA scheme. At the level of data streams we studied the distribution of QRA errors as a function of $`\sigma /q`$ and of the number $`N_\mathrm{s}=3`$, 4, $`\mathrm{}`$, 60 of repeated pixel observations obtained by the PLANCK scanning strategy. As expected, when $`N_\mathrm{s}`$ is large the QRA error is well approximated by a normally distributed random variable, and standard error propagation may be applied.

For large $`N_\mathrm{s}`$ most of the numerical results at the data-stream level may be well approximated by analytical formulae. At map and $`C_{\mathrm{}}`$ level this is not always true, since the effect of the scanning strategy combined with destriping and map-making algorithms can not be modelled in a simple way. For these cases, we have derived semi-analytical relations whose parameters are calibrated through simulations.

The QRA process increases the RMS noise per pixel at the level $`1\%`$ for $`\sigma /q2`$. The distribution of the RMS induced by the QRA process depends on the scanning law and resembles the distribution of white noise in the map. At low and middle ecliptical latitudes, for Planck/LFI the RMS decreases toward the ecliptical poles scaling approximately as $`\sqrt{\mathrm{sin}\theta }`$ where $`\theta `$ is the ecliptical colatitude. The details of the distribution of the QRA additional noise on the map depend on the frequency channel, the location of the horn inside the focal plane and the scanning law. We have derived the parameters that characterize this dependence for Planck/LFI. We find that the skewness and the kurtosis measured all over the map are modified up to a few percent for $`\sigma /q2`$ when the $`1/f^\alpha `$ noise is subtracted with destriping algorithms. Although not large, this effect is not negligible compared to the sensitivity of the kurtosis index recovery for a pure CMB sky. On the other hand, it is much less relevant than residual foreground contamination which accurate modelling turns to set the fundamental uncertainty on the recovery CMB anisotropy high order moments.

Finally, we find that the quantization introduces a power excess, $`C_{\mathrm{}}^{ex}`$, that, although related to the instrument and mission parameters, is weakly dependent on the multipole $`\mathrm{}`$ at middle and large $`\mathrm{}`$ and can be quite accurately subtracted. For $`\sigma /q2`$, the residual uncertainty, $`\mathrm{\Delta }C_{\mathrm{}}^{ex}`$, implied by this subtraction is of only $``$1–2$`\%`$ of the RMS uncertainty, $`\mathrm{\Delta }C_{\mathrm{}}^{noise}`$, on $`C_{\mathrm{}}^{sky}`$ reconstruction due to the noise power, $`C_{\mathrm{}}^{noise}`$.

Only for $`\mathrm{}\stackrel{_<}{_{}}30`$ the quantization removal is less accurate; in fact, the $`1/f`$ noise features, although efficiently removed, increase $`C_{\mathrm{}}^{noise}`$, $`\mathrm{\Delta }C_{\mathrm{}}^{noise}`$, $`C_{\mathrm{}}^{ex}`$ and then $`\mathrm{\Delta }C_{\mathrm{}}^{ex}`$. Anyway, at low multipoles $`C_{\mathrm{}}^{\mathrm{sky}}\mathrm{\Delta }C_{\mathrm{}}^{\mathrm{noise}}>\mathrm{\Delta }C_{\mathrm{}}^{ex}`$ and the uncertainty introduced by the QRA effect is therefore in any case much less than the unavoidable uncertainty due the cosmic variance.

###### Acknowledgements.

We gratefully acknowledge K.M. Górski and all the people involved in the realization of the tools of HEALPix pixelisation. We warmly thank ??? for having provided us with the anisotropy map from thermal SZ effects in clusters of galaxies and G. De Zotti and L. Toffolatti for numberless discussions on foreground contamination.

## Appendix A calculation of $`\alpha _{N_\mathrm{s}}`$

The $`\alpha _{N_\mathrm{s}}`$ coefficient in eq. (13) is defined as the ratio between the standard error for the QRA process after $`N_\mathrm{s}`$ averages $`\stackrel{~}{\delta }_{N_\mathrm{s}}`$ and the RMS for the same process:

$$\alpha _{N_\mathrm{s}}=\frac{\stackrel{~}{\delta }_{N_\mathrm{s}}}{\sigma _{N_\mathrm{s}}},$$
(25)
where $`\sigma _{N_\mathrm{s}}`$ is the RMS of the QRA error distribution after $`N_\mathrm{s}`$ averages: $`u_{N_\mathrm{s}}(X)`$, $`\stackrel{~}{\delta }_{N_\mathrm{s}}`$ is obtained from the definition of standard error:

$$\begin{array}{ccc}I_{N_\mathrm{s}}(\stackrel{~}{\delta }_{N_\mathrm{s}})\hfill & =\hfill & _{\stackrel{~}{\delta }_{N_\mathrm{s}}}^{\stackrel{~}{\delta }_{N_\mathrm{s}}}𝑑tu_{N_\mathrm{s}}(t)\hfill \\ & =\hfill & \frac{1}{\sqrt{2\pi }}_1^{+1}𝑑te^{t^2/2}0.6827.\hfill \end{array}$$
(26)
We are interested in determine the central moments and the integral of $`u_{N_\mathrm{s}}(\overline{\delta }_{N_\mathrm{s}})`$ i.e. to determine the characteristic function of $`u_{N_\mathrm{s}}(\overline{\delta }_{N_\mathrm{s}})`$, from the characteristic function for $`S_{N_\mathrm{s}}=N_\mathrm{s}\overline{\delta }_{N_\mathrm{s}}`$. Where, $`S_{N_\mathrm{s}}`$ represents the distribution of the sum instead of the average of $`N_\mathrm{s}`$ random variables. Taking in account that both distributions are even and for $`N_\mathrm{s}=1`$ they coincide with the top-hat distribution: $`u(x)=1`$, for $`|x|1/2`$, or 0 for $`|x|>1/2`$, it is possible to demonstrate that the characteristic function of $`u_{N_\mathrm{s}}(\overline{\delta }_{N_\mathrm{s}})`$ is:

$$\mathrm{\Phi }_{N_\mathrm{s}}(\omega )=\left[\mathrm{sync}\left(\frac{\omega }{2N_\mathrm{s}}\right)\right]^{N_\mathrm{s}}.$$
(27)
where $`\mathrm{sync}(x)=\mathrm{sin}(x)/x`$. From eq. (27) it is immediate to recover the central moments of $`u_{N_\mathrm{s}}(\overline{\delta }_{N_\mathrm{s}})`$. In particular it is relevant to note that all the odd moments are null, while: $`m_2=\sigma _{N_\mathrm{s}}^2/q^2=1/(12N_\mathrm{s})`$ and $`m_4=(5N_\mathrm{s}2)/(240N_\mathrm{s}^3)=\mathrm{kurtosis}/q^4`$.

The evaluation of $`\stackrel{~}{\delta }_{N_\mathrm{s}}`$ requires to solve eq. (26). This is obtained by the following formula:

$$\begin{array}{cc}I_{N_\mathrm{s}}(\stackrel{~}{\delta }_{N_\mathrm{s}})=\hfill & R(\omega _{\mathrm{max}})+\hfill \\ & \\ & \frac{4\stackrel{~}{\delta }_{N_\mathrm{s}}}{\pi }_0^{\omega _{\mathrm{max}}}\mathrm{sync}\left(\omega \stackrel{~}{\delta }_{N_\mathrm{s}}\right)\left[\mathrm{sync}\left(\frac{\omega }{2N_\mathrm{s}}\right)\right]^{N_\mathrm{s}}𝑑\omega ,\hfill \end{array}$$
(28)
where $`\omega _{\mathrm{max}}`$ represents a cut-off in an otherwise indefinitely extended integration range. The residual $`R(\omega _{\mathrm{max}})`$ is upper bounded by: $`R(\omega _{\mathrm{max}})4(2N_\mathrm{s})^{N_\mathrm{s}}/\left(\pi N_\mathrm{s}\omega _{\mathrm{max}}^{N_\mathrm{s}}\right)`$ which allows to determine $`\omega _{\mathrm{max}}`$ once the absolute accuracy of integration $`ϵ`$ is given:

$$\omega _{\mathrm{max}}=2N_\mathrm{s}\left[\frac{4}{\pi ϵN_\mathrm{s}}\right]^{1/N_\mathrm{s}}.$$
(29)
Note that it would be also possible to solve eq. (26) through direct integration of $`u_{N_\mathrm{s}}(\overline{\delta }_{N_\mathrm{s}})`$, but this is not a satisfactory method since the integral converges slowly for $`N_\mathrm{s}>10`$.

Values of $`\alpha _{N_\mathrm{s}}`$ as a function of $`N_\mathrm{s}`$ from equations (25) and (26) are tabulated, with four digits of accuracy and for $`N_\mathrm{s}=1`$, 2, 3, 4, 5, 6, 10, 15, 30, 60, in the last column of Tab. 5. While the values of $`\stackrel{~}{\delta }_{N_\mathrm{s}}`$, and $`\sigma _{N_\mathrm{s}}`$ are tabulated in the second and the third columns of the table respectively. For $`N_\mathrm{s}30`$ the following relation holds:

$$\mathrm{log}_{10}(\alpha _{N_\mathrm{s}}1)0.923359\mathrm{log}_{10}N_\mathrm{s}1.129003.$$
(30)
with an accuracy compatible with the accuracy of the tabulated values.

## Appendix B the ADC quantization problem

In this paper we have extensively discussed the effect of the software quantization of the data before lossless compression. This quantization is only the last step for the on-board data-processing which may affect the data quality, being the data acquired by the on-board computer by an analog to digital converter (ADC) which introduces at least three sources of noise: $`i)`$ the intrinsic ADC quantization noise proportional to the ADC quantization step $`q_{\mathrm{ADC}}`$; $`ii)`$ the read-out noise $`\sigma _{\mathrm{RON}}`$; $`iii)`$ non linearities due to residual differences in $`q_{\mathrm{ADC}}`$ as a function of the input value seen by the ADC. Detailed analyses of these problems are outside the scope of this paper and will be addressed in forthcoming works. However, in the context of this study on the signal digitisation effect it is relevant to discuss the conditions under which the ADC quantization effect is not critical compared to the software quantization effect.

The impact of the ADC quantization depends on the steps of the on-board data processing. Assuming $`(a)`$ to acquire separately the sky and the reference signal with an ideal ADC (i.e. considering negliglible the sources of noise $`i)`$ and $`ii)`$), $`(b)`$ to have a perfectly balanced radiometer (i.e. both the sky and load detectors have the same weight), $`(c)`$ the software quantization, and then the total QRA error, to be independent from the ADC quantization, and $`(d)`$ to have $`N_\mathrm{h}`$ ADC samples that are averaged to obtain a given instrumental sample (i.e. the hardware sampling rate is $`N_\mathrm{h}`$ times higher than the software sampling rate), then the global RMS quantization error, $`\sigma _g`$, including both hardware and software quantization, is given by:

$$\sigma _g=\sigma _q\sqrt{\frac{2}{N_\mathrm{h}}\left(\frac{q_{\mathrm{ADC}}}{q}\right)^2+1},$$
(31)
where $`\sigma _q`$ is the RMS of the software quantization error. This relation is substantially correct at the level of TODs as well as at the level of the maps, since the repeated observation of the sky pixels and the averaging determined by the scanning strategy operate in the same way in the case of data affected by software quantitazion alone or by hardware plus software quantization, analogously to the considerations of Sect. 2.1. Clearly, $`\sigma _g`$ is larger than the expected QRA error reported by eq. (1) which represents a very good approximation in the limit $`qq_{\mathrm{ADC}}`$ and/or for large values of $`N_\mathrm{h}`$.

In order to understand the conditions under which the ADC quantization effect is not critical compared to the software quantization effect we focus on the impact on the angular power spectrum estimation. At middle and large multipoles the power spectrum of the global and the software quantization errors are proportional to the square of the corresponding RMS quantization error, i.e.:

$$\frac{C_l^{\mathrm{QRA},\mathrm{g}}C_l^{\mathrm{QRA}}}{C_l^{\mathrm{QRA}}}\frac{\sigma _g^2\sigma _q^2}{\sigma _q^2}.$$
(32)
By following an approach analogous the that described in the second part of Sect. 4, we require that

$$\frac{C_l^{\mathrm{QRA},\mathrm{g}}C_l^{\mathrm{QRA}}}{C_l^{\mathrm{QRA}}}\stackrel{_<}{_{}}\frac{1}{\lambda }\frac{\mathrm{RMS}[C_l^{\mathrm{QRA}}]}{C_l^{\mathrm{QRA}}}.$$
(33)
Differently from the case discussed in Sect. 4, we have now to search for a law, $`L_{\mathrm{}}`$, which gives at each $`\mathrm{}`$ the minimum value among the ratios $`\mathrm{RMS}[C_{l,60}^{\mathrm{exc}}]/C_{l,60}^{\mathrm{exc}}`$ found for the set of simulations considered here. We find that the approximate law

$$L_{\mathrm{}}0.46\mathrm{}^{0.5}$$
(34)
is a quite good approximation for the whole range of multipoles relevant here.

By using the above relation between $`\sigma _g`$ and $`\sigma _q`$ and this expression for $`L_{\mathrm{}}`$, the above condition is satisfied provided that

$$\frac{1}{N_\mathrm{h}}\left(\frac{q_{\mathrm{ADC}}}{q}\right)^2\stackrel{_<}{_{}}\frac{0.23\mathrm{}^{0.5}}{\lambda }.$$
(35)
In the case of Planck/LFI, $`N_\mathrm{h}50`$ and even in the worst case of $`q0.5`$ mK and $`q_{\mathrm{ADC}}1.2`$ mK, this condition is satisfied at $`\mathrm{}\stackrel{_>}{_{}}30`$ with a value of $`\lambda 20`$, which corresponds to a deviation from the value of expectation of the angular power spectrum of pure software quantization error of only few per cent of its RMS. Of course, given an accurate description of the ADC quantization we can easily include it in the data analysis. Anyway, this estimate shows that its effect is typically very small.

Finally, it is worth to note that just replacing $`2/N_\mathrm{h}`$ with $`1/N_\mathrm{h}`$ in eq. (31) accounts for a different conceptual acquisition scheme where the sky and reference signals are differenced before to perform the ADC quantization. This scheme has a better propagation of the quantization error than the scheme assumed in eq. (31), implying that the above condition can be satisfied with a value of $`\lambda `$ two times larger.

## Appendix C Statistical moments of simulated sky maps

We report here a comprehensive tabulation of the statistical moments of the CMB simulated map adopted in this work compared with the statistical moments of simulated maps of the most relevant Galactic and extragalactic foregrounds. Our maps of Galactic foregrounds are simulated according to Maino et al. (2002), that of extragalactic source fluctuations according to the model by Toffolatti et al. (1998), while for thermal Sunyaev-Zeldovich effects (SZ) from clusters of galaxies we have used the template available at the MPI web site ???. It is infact interesting to compare the modification on the statistical moments of a pure CMB anisotropy sky due to the digisation effect with that due to the foreground contamination.

The HEALPix scheme (Gòrski et al. 1998) is here adopted and different resolutions, identified by the parameter $`N_{\mathrm{side}}`$ ($`12\times N_{\mathrm{side}}^2`$ is the number of pixels in the map), are considered.

We consider here maps with and without beam convolution (FWHM of $`33^{}`$ as appropriate to the LFI 30 GHz channels). The input maps have been first simulated at $`N_{\mathrm{side}}=1024`$ without beam convolution. Convolution and, possibly, degradation are then applied.

We report the statistical moments referring separately to each component and to combination of CMB and foreground maps, by including or not masks to exclude regions at low Galactic latitudes and/or high signal pixels, as indicated in the tables. In fact, the regions at low Galactic latitudes have to be avoided in CMB anisotropy analysis as well as pixels significantly contaminated by SZ effects or extragalactic sources have to be previously detected and removed. More precisely, for uniformity reasons, we apply the same 5$`\sigma `$-clipping threshold, $`\sigma `$ being the rms of CMB fluctuations at the considered resolution for the corresponding case (including or not beam convolution), independently of the considered kind of map composition.

Note that, while in the cases of unconvolved maps with large skewness and kurtosis indices the convolution decreases both these estimators, in the cases of unconvolved maps with moderate or small skewness and kurtosis indices the convolution may produce a weak increase of these estimators because the higher order moments are relatively less reduced than the variance.


format: text

text:
# On deciding whether a Boolean function is constant or not 11footnote 1Work supported in part by Istituto Nazionale di Fisica Nucleare, Sezione di Trieste, Italy

## 1 Introduction

Query complexity theory is mainly concerned with the computational cost required to determine some specific property of functions. The cost is measured by the number of queries which can be addressed to a “black box” device (the oracle) that outputs instantaneously an answer to a query. Various results concerning query complexity of Boolean functions, i.e. functions with a finite-size domain which have two possible output values only, exist in the literature . The first of this sort of problems which has been addressed by the quantum computational theory, is represented by the Deutsch-Josza algorithm . It involves the decision whether a given Boolean function of $`n`$ binary digits is constant (that is whether it outputs the same binary digit for every input value) or balanced (that is whether it outputs $`1`$ on exactly half inputs and $`0`$ on the other half) and it can be solved with certainty by querying a quantum oracle just once.

This task is successfully achieved by exploiting two essential features of quantum mechanics: its linearity, which allows the simultaneous evaluation of a function on linear superpositions of its input values, and the interference between amplitudes, which raises the probability of obtaining the desired results.

The a-priori knowledge that the function computed by the quantum oracle is either constant or balanced, is crucial in solving the problem exactly, i.e. with zero probability of error. In this paper, we address the question of deciding the constancy of a generic function chosen at random among the set of $`2^{2^n}`$ $`n`$-bit Boolean functions. In such a case, the iteration of the same algorithmic procedure permits one to decide about the constancy of such generic function with a probability of error which depends on $`n`$ and on the number $`k`$ of quantum oracle-queries.

Unfortunately, this probability tends to $`1`$ in the worst case scenario, namely when the given function outputs the same value for every input but one. Nevertheless, it is our purpose to determine the efficiency of the quantum algorithm and the range of values of $`n`$ and $`k`$ for which it is more efficient compared to the unique classical procedure which involves querying the oracle successively until two different outputs are encountered.

Not too surprisingly, due to classical conditioning of each output on the previous ones, in some specific situations (the worst cases), for large number of inputs and nearly as many oracle-queries, the efficiency of the classical algorithm turns out to overcome that of the quantum one. Moreover, in the average case, i.e. when the $`n`$-bit Boolean function is sorted at random, the quantum algorithm is always preferable to the classical one.

## 2 The worst case analysis

The general problem we are going to tackle is as follows. Suppose that we are given an oracle which instantaneously computes an arbitrary n-bit Boolean function $`f:\{0,1\}^n\{0,1\}`$ on inputs that, as usual, will be enumerated by the integers $`\{0,1,\mathrm{},2^n1\}`$. Our goal is to devise a quantum algorithm, involving oracle queries, in order to determine, with some probability of error, whether such a function is constant or not.

Since, by hypothesis, we have no a-priori information or knowledge about $`f`$, we can only suppose that it has been sorted completely at random within the set $`𝒮_n`$ of $`2^{2^n}`$ n-bit Boolean functions.

If we have at our disposal a classical oracle which outputs the value of $`f`$, and which can be queried at most $`k`$ times, the only conceivable classical algorithm able to solve the problem is : “query the oracle $`k`$ times or until two output values are different: in the latter case stop querying and declare that the function is not constant, otherwise say that $`f`$ is constant”.

Of course, this simple algorithm fails to give the right answer only when the function is actually non-constant and when we have obtained $`k`$ consecutive equal digits in the querying procedure. In fact, the function $`f`$ might output a different digit on any of the untried inputs, so that one has to evaluate the probability of such an event.

We start by considering the worst possible scenario, that is the unknown function $`f`$, drawn from the set $`𝒮_n`$, actually outputs the same bit, say $`0`$, for all of its inputs but one (for simplicity of notation we will say that such a function is of $`f_1`$-type). The probability $`p_1(k,n)`$ of declaring erroneously constant such a function, coincides with the probability of obtaining $`k`$ consecutive digits $`0`$ after querying the oracle $`k`$ times:

$$p_1(k,n)=\text{Prob}[\underset{k}{\underset{}{0\mathrm{}0}}]=\text{Prob}[0]\left(\underset{\mathrm{}=1}{\overset{k1}{}}\text{Prob}[\mathrm{\hspace{0.17em}0}|\underset{\mathrm{}}{\underset{}{0\mathrm{}0}}]\right),$$
(2.1)
where we have indicated with $`\text{Prob}[\mathrm{\hspace{0.17em}0}|\underset{\mathrm{}}{\underset{}{0\mathrm{}0}}]`$, $`1\mathrm{}k1`$, the conditional probability of getting $`0`$ after the $`(\mathrm{}+1)`$-th oracle query if a sequence of $`\mathrm{}`$ consecutive zeroes has already been obtained. Since the function $`f`$ is of $`f_1`$-type, one easily calculates $`\text{Prob}[0]=(2^n1)/2^n`$ and $`\text{Prob}[0|\underset{\mathrm{}}{\underset{}{0\mathrm{}0}}]=(2^n\mathrm{}1)/(2^n\mathrm{})`$, whence

$$p_1(k,n)=\underset{j=1}{\overset{k}{}}\frac{2^nj}{2^n(j1)}=1\frac{k}{2^n}.$$
(2.2)
This probability distribution is equal to zero for $`k=2^n`$; indeed, due to statistical correlations between subsequent outputs, in such a case we can determine with certainty if the function is constant or not. As we will see, this is not the case with the Deutsch-Josza quantum algorithm which, unless the function is either constant or balanced, is not able to give the correct answer with certainty. First, let us briefly recall the basics of that algorithm whose associated quantum circuit is the following:

We have indicated with $`H^n`$ the $`n`$-fold Hadamard transformation on a $`n`$-qubit state and with $`U_f`$ the quantum oracle associated to $`f`$. It is defined by $`U_f|x|y=|x|yf(x)`$, where $``$ denotes the sum modulus $`2`$, $`|x`$ is an $`n`$-qubit state such that $`x\{0,1,\mathrm{},2^n1\}`$, and $`|y`$ is a one-qubit state. If the initial input states are those indicated in the figure, it is possible to prove  that the $`n`$-qubit output state $`|\psi `$ in the figure above, turns out to be :

$$|\psi =\underset{z=0}{\overset{2^n1}{}}\left(\underset{x=0}{\overset{2^n1}{}}\frac{(1)^{zx+f(x)}}{2^n}\right)|z,$$
(2.3)
where $`zx`$ is the bitwise inner product of $`z`$ and $`x`$, i.e. $`zx=_{j=1}^nz_jx_j`$. Therefore, according to Eq. (2.3), a measurement process performed onto the $`n`$-qubit computational basis $`\left\{|z\right\}`$ yields the result $`z=0`$ with the following probability :

$$\text{Prob}[z=0,|\psi ]=\left|\frac{1}{2^n}\underset{x=0}{\overset{2^n1}{}}(1)^{f(x)}\right|^2.$$
(2.4)
If the function $`f`$ is constant, the result $`z=0`$ is obtained with certainty; this in turn implies that, if any other value $`z0`$ is obtained, the function cannot be constant. Yet, if the function is of $`f_1`$-type, the result $`z=0`$ is obtained with a probability strictly smaller than $`1`$, given by:

$$\text{Prob}[z=0,|\psi |f=f_1]=\left|\frac{1}{2^n}\underset{x=0}{\overset{2^n1}{}}(1)^{f(x)}\right|^2=\left|\frac{1}{2^n}(2+2^n)\right|^2=\left(1\frac{1}{2^{n1}}\right)^2.$$
(2.5)
The simplest quantum procedure to determine whether the function $`f`$ is constant or not, consists of a simple iteration of the quantum circuit described above and it amounts to: “iterate the Deutsch-Jozsa algorithm $`k`$ times or until a measurement result $`z0`$ appears: in the latter case stop querying and declare that the function is not constant, otherwise say that $`f`$ is constant”.

We denote by $`q_1(k,n)`$ the probability that, after obtaining $`z=0`$ in $`k`$ consecutive measurements, we are wrong in asserting the constancy of the $`f_1`$-type function. According to Eq. (2.5), one gets:

$$q_1(k,n)=\left(1\frac{1}{2^{n1}}\right)^{2k}.$$
(2.6)
It is worth noting that, contrary to the classical case where the probability of obtaining the same output value changes, by conditioning, after each oracle query, the consecutive (quantum) oracle-queries are independent from each other. Therefore, the probability of obtaining $`z=0`$, $`k`$ times consecutively, is the product of the probabilities of the uncorrelated events; quantum conditioning on past results would certainly improve the efficiency of the algorithm, however we do not know of any such technique.

In order to compare the quantum querying procedure with respect to the classical one in this worst case, we define the following efficiency function $`\mathrm{\Delta }_1(k,n)`$:

$$\mathrm{\Delta }_1(k,n)p_1(k,n)q_1(k,n)=1\frac{k}{2^n}\left(1\frac{1}{2^{n1}}\right)^{2k}.$$
(2.7)
Given $`n`$ (the number of bits) and $`k`$ (the number of oracle-queries), a positive value of $`\mathrm{\Delta }_1(k,n)`$ indicates that the quantum algorithm is more efficient than the classical one. For $`n2`$ and $`1k2^n1`$, $`\mathrm{\Delta }_1(k,n)`$ attains its absolute maximum at $`n=2`$ and $`k=1`$:

$$\underset{k,n}{\mathrm{max}}\mathrm{\Delta }_1(k,n)=\mathrm{\Delta }_1(k=1,n=2)=0.5.$$
(2.8)
For fixed $`n`$, the function $`\mathrm{\Delta }_1(k,n)`$ has relative maxima at

$$k^{}0.5\frac{\mathrm{log}\left(\frac{2^{1n}}{\mathrm{log}(12^{1n})}\right)}{\mathrm{log}(12^{1n})}\mathrm{\hspace{0.25em}\hspace{0.25em}0.35}2^n\text{for}\mathrm{\hspace{0.25em}\hspace{0.25em}2}^n1,$$
(2.9)
with $`\mathrm{\Delta }_1(k^{},n)0.40`$, for $`2^n1`$.

In Fig. 1, $`\mathrm{\Delta }_1(k,n)`$ is plotted as a function of $`k`$ for different values of $`n`$; the curves displayed from left to right are associated with increasing values of $`n`$, $`5n12`$:

As it can be argued from the picture plotted above, $`\mathrm{\Delta }_1(k,n)`$ is a concave function in the variable $`k`$ for every fixed value of $`n3`$. Therefore, the worst case efficiency of the quantum algorithm increases, with respect to the classical one, for $`kk^{}`$ (where $`k^{}0.352^n`$ for $`2^n1`$), until it reaches its (relative) maximum value. Such a value exhibits a weak dependence on the variable $`n`$ and it is approximately equal to $`0.40`$ for $`2^n1`$. On the contrary, for $`k>k^{}`$, $`\mathrm{\Delta }_1(k,n)`$ decreases to zero, thus indicating a substantial equivalence of the two methods for large $`k`$’s.

The above results show that, in the worst case, it is almost always better to resort to the quantum rather than to the classical algorithm. However, the improvement of quantum queries over the classical ones decreases with the number of queries. Furthermore, $`\mathrm{\Delta }_1(k,n)`$ becomes negative for $`k2^n`$, namely there is a (small) range of values of $`k`$ for which the classical querying procedure gives less probability of error than the quantum one.

Not surprisingly, this is due to the fact that the $`2^n`$-th classical query is able to determine with certainty if the function is constant or not, since $`p_1(k=2^n,n)=0`$. Instead, the quantum algorithm, as it has been conceived, can never answer with certainty since $`q_1(k,n)`$ never vanishes.

## 3 The average case analysis

Until now we have restricted our analysis to the worst possible scenario, i.e. when the function $`f`$ outputs the same value for all possible inputs but one. In the general case, the function may be any of the $`2^{2^n}`$ possible $`n`$-bit Boolean functions belonging to the set $`𝒮_n`$. This set contains $`2`$ constant functions (one always outputs $`0`$, while the other outputs $`1`$), $`\left(\genfrac{}{}{0pt}{}{2^n}{2^{n1}}\right)`$ balanced functions (they output the same value for exactly one-half of the possible inputs) and $`2\left(\genfrac{}{}{0pt}{}{2^n}{2^nm}\right)`$ $`f_m`$-type functions (they output $`0`$ ($`1`$) on $`2^nm`$ inputs and $`1`$ ($`0`$) on the remaining $`m`$ ones).

We repeat the analysis of the previous section and suppose we are given one particular $`f_m`$-type function. The probability of obtaining always the same output bit by querying the classical oracle $`k`$ times, and thus of concluding, erroneously, that the function is constant, is equal to $`\text{Prob}[\underset{k}{\underset{}{0\mathrm{}0}}]+\text{Prob}[\underset{k}{\underset{}{1\mathrm{}1}}]`$. The latter are the probabilities of obtaining $`k`$ consecutive digits $`0`$ and $`1`$ respectively, and they explicitly read

$$\text{Prob}[\underset{k}{\underset{}{0\mathrm{}0}}]=\underset{j=0}{\overset{k1}{}}\frac{2^nmj}{2^nj}=\frac{(m2^n)_k}{(2^n)_k},$$
(3.1)
$$\text{Prob}[\underset{k}{\underset{}{1\mathrm{}1}}]=\underset{j=0}{\overset{k1}{}}\frac{mj}{2^nj}=\frac{(m)_k}{(2^n)_k}.$$
(3.2)
In the above, for simplicity of notation, the Pochhammer symbol $`(a)_na(a+1)\mathrm{}(a+n1)`$ has been used. Therefore, the classical error probability $`p_m(k,n)`$ in declaring constant an $`f_m`$-type function is:

$$p_m(k,n)=\frac{(m2^n)_k}{(2^n)_k}+\frac{(m)_k}{(2^n)_k}.$$
(3.3)
The corresponding probability $`q_m(k,n)`$ for the quantum case is much easier to calculate. In fact, given an $`f_m`$-type function, the conditional probability of obtaining $`k`$ consecutive measurement results $`z=0`$ via the iterated Deutsch-Jozsa procedure is, according to Eq. (2.4):

$$\text{Prob}[k\text{times}z=0,|f=f_m]q_m(k,n)=|\frac{1}{2^n}\underset{x=0}{\overset{2^n1}{}}(1)^{f(x)}|^{2k}=(1\frac{m}{2^{n1}})^{2k}.$$
(3.4)
As in the previous section, we can now define a function $`\mathrm{\Delta }_m(k,n)`$ which quantifies the relative efficiency of the quantum algorithm with respect to the classical one:

$$\mathrm{\Delta }_m(k,n)p_m(k,n)q_m(k,n)=\frac{(m2^n)_k}{(2^n)_k}+\frac{(m)_k}{(2^n)_k}\left(1\frac{m}{2^{n1}}\right)^{2k}.$$
(3.5)
$`\mathrm{\Delta }_m(k,n)`$ does not dramatically change with $`n`$; therefore, we fix $`n=7`$ and present typical behaviours for different values of $`m`$.

We notice that $`\mathrm{\Delta }_m(k,n)`$ is symmetric around $`m=2^{n1}`$, and that, when $`m`$ gets close to $`2^{n1}`$, the function $`f_m`$ tends to be balanced and $`\mathrm{\Delta }_m(k,n)`$ rapidly decreases. This behaviour indicates that the quantum algorithm can perform much better than the classical one, but only for small numbers of oracle-queries. On the contrary, the two procedures become equally efficient for larger values of $`k`$. For instance, in Fig. 2, when $`m=60`$ and the function $`f_m`$ is almost balanced, the two procedures are nearly indistinguishable for $`k>5`$. Instead, in the general case of $`m`$ significantly different from $`2^{n1}`$, the dependence of $`\mathrm{\Delta }_m(k,n)`$ on $`k`$ is much smoother and presents a maximum and a change of concavity. Further, for $`k2^n`$, $`\mathrm{\Delta }_m(k,n)`$ assumes negative values, invisible in the above pictures. This indicates a better performance of the classical algorithm over the quantum one, as already noticed when dealing with the $`f_1`$-type function in the previous section. However this phenomenon is much less pronounced than in the worst case situation.

### 3.1 Randomly chosen Boolean functions

We now tackle the general case and suppose we are given a randomly chosen Boolean function belonging to $`𝒮_n`$. Then, the average classical error probability $`\overline{p}(k,n)`$ in declaring such a function constant, after $`k`$ consecutive queries with identical outputs, is:

$$\overline{p}(k,n)=2\left[\underset{m=1}{\overset{2^{n1}1}{}}\frac{\left(\genfrac{}{}{0pt}{}{2^n}{m}\right)}{2^{2^n}}\left(\frac{(m2^n)_k}{(2^n)_k}+\frac{(m)_k}{(2^n)_k}\right)\right]+\frac{\left(\genfrac{}{}{0pt}{}{2^n}{2^{n1}}\right)}{2^{2^n}}\frac{2(2^{n1})_k}{(2^n)_k}.$$
(3.6)
In the previous expression, the factor $`\left(\genfrac{}{}{0pt}{}{2^n}{m}\right)/2^{2^n}`$, with $`m2^{n1}`$, gives the probability of sorting, within the set $`𝒮_n`$, an unbalanced $`f_m`$-type function which outputs a number of zeroes less than $`2^{n1}`$ (the extra factor $`2`$ in front of the sum takes into account the contribution coming from functions which output a number of zeroes which is greater than $`2^{n1}`$). Analogously, the factor $`\left(\genfrac{}{}{0pt}{}{2^n}{2^{n1}}\right)/2^{2^n}`$ gives the probability of sorting a balanced function.

On the contrary, the corresponding quantum error probability $`\overline{q}(k,n)`$ of choosing at random a Boolean function within $`𝒮_n`$ and declaring it constant, due to $`k`$ consecutive measurement results $`z=0`$, is:

$$\overline{q}(k,n)=\underset{m=1}{\overset{2^n1}{}}\frac{\left(\genfrac{}{}{0pt}{}{2^n}{m}\right)}{2^{2^n}}\left(1\frac{m}{2^{n1}}\right)^{2k}.$$
(3.7)
As before, we introduce the average efficiency function $`\overline{\mathrm{\Delta }}(k,n)`$:

$$\overline{\mathrm{\Delta }}(k,n)\overline{p}(k,n)\overline{q}(k,n)=2\left[\underset{m=1}{\overset{2^{n1}1}{}}\frac{\left(\genfrac{}{}{0pt}{}{2^n}{m}\right)}{2^{2^n}}\left(\frac{(m2^n)_k}{(2^n)_k}+\frac{(m)_k}{(2^n)_k}\right)\right]+$$

$$\frac{\left(\genfrac{}{}{0pt}{}{2^n}{2^{n1}}\right)}{2^{2^n}}\frac{2(2^{n1})_k}{(2^n)_k}\underset{m=1}{\overset{2^n1}{}}\frac{\left(\genfrac{}{}{0pt}{}{2^n}{m}\right)}{2^{2^n}}\left(1\frac{m}{2^{n1}}\right)^{2k}.$$
(3.8)
This function quantifies the (average) relative efficiency of the quantum algorithm with respect to the classical one: the larger $`\overline{\mathrm{\Delta }}(k,n)`$, the more preferable the quantum procedure with respect to the classical one.

The dependence of $`\overline{\mathrm{\Delta }}(k,n)`$ on $`k`$ and $`n`$ is shown in Fig. 3 below, where its different behaviour with respect to $`\mathrm{\Delta }_m(k,n)`$ in Fig. 2 emerges.

It is worth commenting on the behaviour of $`\overline{\mathrm{\Delta }}(k,n)`$. First, we notice that this function is approximately zero for most of the values of $`n`$ and $`k`$, indicating a substantial equivalence of the two algorithms; indeed, in all but a few situations, they are able to solve the problem with negligible probability of error.

Second, for fixed $`n`$, $`\overline{\mathrm{\Delta }}(k,n)`$ is monotonically decreasing in $`k`$; for $`k=1`$, the quantum error probability $`\overline{q}(k=1,n)`$ falls rapidly to zero, while the classical one, $`\overline{p}(k=1,n)`$, approaches unity.

The decrease to zero of $`\overline{q}(k=1,n)`$ is due to the fact that the number of balanced (or almost balanced) functions increases considerably for $`n`$ large and the quantum algorithm is able, as we already know, to recognize them with certainty (or almost certainty). On the contrary, the number of worst case functions, i.e. those which output almost always the same value, gets rapidly to zero and the probability of sorting one of them becomes negligible. This argument does not apply to the classical probability $`\overline{p}(k=1,n)`$, which instead tends to $`1`$; indeed, the classical algorithm fails to give the right answer for balanced (or almost balanced) functions which form the overwhelming majority in $`𝒮_n`$.

Finally, unlike the previous cases, the numerical analysis indicates no negative values of $`\overline{\mathrm{\Delta }}(k,n)`$; thus, in the average case, the quantum algorithm is expected to be always better than the classical one.

## 4 Conclusions

In this paper we have shown how a trivial iteration of the original Deutsch-Jozsa algorithm, can be applied to solve the general problem of deciding if a function $`f`$, sorted at random among the set of $`n`$-bit Boolean functions, is constant or not. We have compared the error probability of the quantum procedure with that of the classical one, the latter consisting of consecutive oracle-queries performed until two different outputs are encountered. We have analyzed the range of values of $`n`$ and $`k`$ for which the quantum procedure is more efficient than the classical one, both on average and in the worst possible situation. From a numerical analysis, one concludes that, by iterating the quantum algorithm, one can always solve the problem, on average, better than with the classical method, and much better in those situations where the number of queries $`k`$ is small and $`2^n`$ large. Moreover, we have seen that, when $`f`$ is a $`f_m`$-type function with $`m`$ not much greater than $`1`$ (the worst cases), then the classical method is asymptotically (that is, when $`2^n`$ is large and $`k2^n`$) preferable to the quantum one. This peculiar effect is due to the fact that the Deutsch-Jozsa algorithm, unlike the classical querying, can never solve the general problem with certainty since it does not allow for conditioning on previous results.