THE ATACAMA COSMOLOGY TELESCOPE: DATA CHARACTERIZATION AND MAPMAKING

For our purposes, we will define a point source as any object with an angular size comparable to or smaller than the beam size of the telescope. The beam full width at half-maximum (FWHM), θ_1/2, is 137 for 148 GHz. Planets are the most important point sources as their high surface brightness makes them useful for both calibration and beam measurements. Distant galaxies can also be approximated as point sources and are helpful for the pointing solution.

As the telescope scans and the sky rotates, a single detector traces a zigzag on the sky, sampling the point source every time it intersects with the beam. The number of times a point source appears depends on the scan period, the beam size, and the central scan azimuth. A point source appears as a succession of "blips" in the time stream. The shape of these blips is a slice through the telescope point-spread function and depends upon the angular separation between the center of the beam and the location of the source. The angular speed of the scan on the sky is given by $\dot{\theta } = v_{\mathrm{scan}}\cos {(50{\mbox{$.\!\!^\circ $}}5)}$. The 2008 scan speed of v_scan = 15 s⁻¹ implies an angular sky speed of 572 s⁻¹. Assuming a Gaussian-like beam of equivalent width and neglecting the transit speed of the source, the 3 dB cutoff frequency is

$\begin{equation} f_{-3\,\mathrm{d}\mathrm{B}} = \frac{2\ln {(2)}\:\dot{\theta }}{\pi \theta _{1/2}}. \end{equation} \tag{ 1 }$

This means that the contribution from point sources to the TOD is limited to frequencies below 18 Hz, as can be seen in Figure 3.

Zoom In Zoom Out Reset image size

Given our scan speed, nominal elevation, sky rotation, and sampling rate, we expect nearly 10 samples per beam on the sky at 148 GHz. The time constant of the detectors is fast enough that its effect in the beam shape is small, as will be discussed in Section 9.

Using CMB simulations, we can estimate the expected TOD response to extended sky features. Averaging the power spectra from many such synthetic TODs, we can estimate the CMB power spectrum in TOD space. Figure 3 shows the average power spectrum of a simulated 148 GHz observation of the CMB sky and the average data power spectrum from one 15 minute file. Because of Silk damping, the CMB has a relatively sharp cutoff in its characteristic size: unlike the point-source observations shown in the figure, the CMB has little signal at frequencies higher than 10 Hz.

Given ACT's angular scan speed, the TOD frequency associated with a sky feature of angular size θ_f is

$\begin{equation} f_{\rm TOD} = \frac{\dot{\theta }}{2\theta _f}, \end{equation} \tag{ 2 }$

where the factor of two in the denominator comes from considering that the angular scale is the size of a positive or negative temperature bump, which is half of a wavelength on the sky. Given that multipole moment ℓ relates to angular scale as ℓ ≃ π/θ_f for scales small compared to the full sky, an approximate conversion between multipole and TOD frequency is

$\begin{equation} \ell \simeq \frac{2\pi f_{\rm TOD}}{\dot{\theta }} = \frac{2\pi f_{\rm TOD}}{v_{\mathrm{scan}}\cos {(50{\mbox{$.\!\!^\circ $}}5)}}. \end{equation} \tag{ 3 }$

Note that this relation is only a rule of thumb; the actual mapping of TOD frequency into multipole space depends on the specific scan strategy. This relation shows that the CMB power spectrum in the TOD can be shifted in frequency by changing the scan rate. As a reference, the TOD frequency corresponding to ℓ = 3000 was 7.9 Hz in the 2008 season.

The current signals from the bolometers are amplified in 100-SQUID series arrays (SAs) operated at 3 K (Swetz et al. 2011). Slow temperature changes in the SAs produce slow drifts in the TODs. During the first 10 hr of the night, the SA temperature drifts down by 250 mK, and rises up by 150 mK in the last 2 hr. In terms of equivalent sky temperature at 148 GHz, this corresponds to a drift of nearly 4 K in signal. In frequency space, the drift imprints a 1/f signature on the data, which meets the detector noise level at a knee of nearly 1 Hz. Despite small differences in the responses of different SA amplifiers, most of this signal appears as a common mode to all the detectors. As the coupling occurs at the readout circuit, this signal is also present in the dark detectors, simplifying its identification and eventual removal.

In most cases this signal is subdominant to the atmosphere signal, but when the PWV is low enough the thermal drift becomes comparably significant. This is not evident in Figure 4(a) since the examples of low-atmospheric power get averaged down when grouping many TODs.

In contrast, drifts in detector temperature cause no measurable effect because it is servo-controlled to within less than 1 mK error.

Mechanical accelerations, occurring mainly at the scan turnaround, can couple optically and thermally to the detector and result in an undesired instrumental response. The optical coupling can be produced by a mechanical motion of the detector coupling layer, which is a 40 μm thick silicon layer placed 100 μm away from the detectors for optical impedance matching. The coupling efficiency is sensitive to the distance between the coupling layer and the detectors, so small vibrations can cause significant effects, especially at higher frequency bands where the coupling layer is thinner and the coupling efficiency is more sensitive to changes in the distance. The effect is also larger near the center of the detector array, presumably because the coupling layer vibrates in its fundamental mode. In the TOD, vibrations manifest themselves as a series of spikes visible in the common mode at the scan turnarounds, with opposite signs for opposite directions, leading to lines in the power spectra at several harmonics of the scan frequency. In the 148 GHz "waterfall plot" of Figure 6, the vibrations contribute to the scan harmonics and some resonant lines at higher frequencies. The latter are also seen for stare observations, and their central frequencies are shared among all three arrays, suggesting that they correspond to natural frequencies of the whole system. To mitigate the mechanical effect, the turnaround acceleration was reduced from its value in the 2007 season to the value shown in Table 1 after 2008 October 8. The coupling layers for 218 GHz and 277 GHz were removed after the 2008 season.

Zoom In Zoom Out Reset image size

Thermal perturbations of the cryostat, revealed in spectral analysis of the bath temperature, also add to the low frequency scan harmonics seen in Figure 6. This effect is column dependent because the coupling occurs at the readout circuit level and differs from the coupling layer effect in that the signal is not stronger at the center of the array.

In general, high frequency spectral lines in the TOD from different detectors destructively interfere when projected into map space, mostly canceling out when many observations are combined. On the other hand, low-frequency harmonics of the scan may produce non-negligible bar-like features in the maps perpendicular to the scan direction, if not treated properly.

Electromagnetic pickup couples to the readout circuit, producing various signatures in the data. These signals appear correlated among subgroups of detector TODs, particularly among detectors from the same column or row. We call these "column" or "row" correlations, respectively, each corresponding to a different source of electromagnetic pickup. Narrowband signals appear correlated among detectors within the same column. They couple in somewhere after the first SQUID stage of the readout circuit, as the subsequent circuitry is shared by the detectors in a column. On the other hand, broadband signals appear correlated among detectors within the same row or even a few rows apart. We believe these are caused by rapidly varying signals, such as spikes from strong current transients: given our time-multiplexed readout scheme, in which the same row is read from all columns simultaneously, these signals become correlated among detectors in the same row.

Figure 7 shows the detector correlation matrix for a given 15 minute TOD file, where every element in the matrix corresponds to the correlation between the TODs of detectors i and j, which index the elements of each axis. For instance, the diagonal elements correspond to i = j, so they are all equal to one. In Figure 7(a), adjacent detectors on each axis belong to the same column in the array, with a black line separating detectors from different columns, while Figure 7(b) is organized by rows. We can clearly see the kinds of correlations described above affecting either columns or rows of detectors. As would be expected from a stochastic distribution of transients, the patterns seen in the row-correlated matrix change in time. To quantify the correlations, a quality factor is defined as the mean of the squared off-diagonal elements of the correlation matrix:

$\begin{equation} Q = \frac{2}{N(N-1)}\sum _{i> j}^N{c^2_{i,j}}, \end{equation} \tag{ 6 }$

where c_{i, j} is the correlation between the TODs of detectors i and j, and N is the number of detector TODs. Lower Q factors indicate less correlated noise.

Zoom In Zoom Out Reset image size

The column-correlated signals were significantly reduced by carefully shielding the system from electromagnetic noise in the environment. The row-correlated signals, instead, were found to be related to the switching power supplies feeding the readout electronics. They were replaced by linear power supplies at the end of the 2008 season.

If not treated properly, broadband row-correlated signals can produce features in the maps, mostly seen as features perpendicular to the scan direction for that particular observation.

The readout circuit of the detectors uses SQUIDs, which are very sensitive to magnetic fields. Despite significant magnetic shielding, some residual pickup remains as the telescope sweeps through Earth's magnetic field. Although the scan is almost triangular, the observed signal is chiefly sinusoidal, which can be explained by eddy current losses and hysteresis in the magnetic shielding. The magnetic signal amplitude and phase are fairly stable across detectors in the same column, but differ significantly between columns. We believe that these phase shifts are related to the complexity of SQUIDs and electrical loops in the system. A typical amplitude of this signal is equivalent to nearly T ≃ 7 mK in CMB units.

As this is a readout-related signal, it can be measured using the dark detectors and removed by fitting and subtracting a sinusoid. Moreover, the scan frequency is significantly below the science band and the magnetic pickup is heavily down weighted in the mapping.

The nightly variations in atmospheric loading and daily cryogenic cycling necessitate a rebiasing of the detectors at the beginning of each night's observations, which affects the detector response. Additionally, the opacity of the atmosphere varies from night to night. Changes in loading and opacity throughout the night produce smaller variations in the system response.

In order to account for the variations in system response between nights, a responsivity calibration is derived from the analysis of I − V curves taken during the nightly rebiasing. An I − V curve is the relation between the output current response (I) of the detector's SQUID readout circuitry to a slowly ramping detector bias voltage (V). For the 148 GHz band, the median deviation in responsivity between nights as measured by the I − V curves is 3.0%.

A second test, known as a bias step (Niemack 2008), is performed three times per night to detect changes in detector responsivity through the night. This is achieved by recording the detector response to a series of small, square-wave pulses applied to the detector bias voltage. The responsivities can be obtained by analyzing the detector responses to this signal and are found to be in agreement with responsivities obtained from the I − V curve analysis. For the 148 GHz band, the median deviation in responsivity over a night, as measured by the bias step, is 1.0%. For more information about the ACT detectors, the biasing routine, and responsivity, see Fisher (2009), Battistelli et al. (2008), Swetz et al. (2011), and Zhao et al. (2008). The time dependence on the atmospheric conditions will be discussed below together with the overall system calibration.

To determine the detectors' relative gain coefficients, we use the large common-mode signal provided by the variations in atmospheric brightness. This is analogous to a flat-field calibration done with an optical CCD. For each 15 minute data file, we compute the gain factor that best fits the detector drift to the common mode drift. The dispersion of this factor per detector over the season, averaged over all detectors, is better than 2% rms for the 148 GHz array. This includes the variability of the time-dependent calibration step. We correct for this by multiplying each detector time stream by its corresponding gain factor.

The conversion to sky temperature also depends on the atmospheric transmission. This is estimated with the ATM model (Pardo et al. 2001), which uses PWV measurements from APEX corrected for the ACT site and other Atacama-specific parameters. When PWV measurements are not available, the season-average transmission ($\mathcal {T} = 0.976$) is used. For the 148 GHz band, the rms of the transmission during the 2008 season was 2.3%. Given that large atmosphere loading can excite non-linearities in the detector responsivities, an extra degree of freedom is given to the planet fit as a function of PWV, producing a final correction to the calibration factor. The fit is done comparing Uranus flux measurements obtained under different atmospheric conditions. The resulting transmission yields

$\begin{equation} \mathcal {T}(\mathrm{w}) = \exp \left[ -\left(\tau _d + \tau _w\,w + \tau _x\left(w - \overline{w}\right)\right)/\sin \theta _{\mathrm{alt}}\right], \end{equation} \tag{ 7 }$

where θ_alt is the observation altitude, τ_d = 0.0093 is the "dry opacity," τ_w = 0.0190 mm⁻¹ is the "wet opacity," τ_x = 0.0138 mm⁻¹ is the fit parameter, w is the PWV measurement from APEX, corrected for the ACT site in millimeters, and $\bar{w} = 0.44\,\mathrm{m}\mathrm{m}$ is a pivot PWV used in the fit. Here, both τ_d and τ_w are fixed parameters provided by the ATM model.

The overall map calibration is compared to the WMAP seven-year map, by correlating multipoles in the range 400 < ℓ < 1000, providing an uncertainty of 2% in temperature (Hajian et al. 2011). Putting all pieces together, the overall system calibration for 148 GHz in the 2008 season is

$\begin{equation} \mathcal {C}(\mathrm{w}) = \frac{C}{\mathcal {T}(\mathrm{w})}, \end{equation} \tag{ 8 }$

where C = 19.41 K/pW.

The relative pointing between detectors was determined by modeling the beam in the time streams of planet observations. This analysis made use of approximately 30 observations of Saturn at the normal CMB observing altitude of 505, in two azimuth ranges corresponding to the rising and setting of the planet.

The telescope scans are slow enough that each detector samples the vicinity of the peak response to a planet several times in a single observation. It is thus possible to model the beam position and shape in two spatial dimensions, project this to a detector signal using the telescope encoder information, and fit the data in the time domain. For simplicity, we use a two-dimensional Gaussian as the beam model. The fit produces azimuth and altitude offsets for each detector relative to the telescope pointing encoders, as well as measures of the peak response, beam FWHM, and optical time constants.

When combining different planet observations, we first aligned them by applying an offset correction to each one. The position of each detector is taken as the mean of the positions from all observations, after rejecting outliers. This produces relative detector positions for each array that are in good agreement with design expectations (Fowler et al. 2007). A comparison between the offsets for rising and setting observations (which differ in azimuth by approximately 95°) shows no significant rotation or shear of the array relative to the local altitude and azimuth axes. This constrains the tilt in the telescope azimuth axis to be smaller than 1' and indicates that rotation of the array with respect to azimuth is negligible.

The uncertainty in the relative pointing is no greater than 15 for all three arrays. Given the large number of detectors in each array, this small uncertainty in the relative detector positions amounts to a negligible contribution to the pointing error in the final maps.

The azimuth and altitude of the telescope encoder readings must be corrected to account for their offset from the true boresight for each frequency band. The correction is different for each of the four central CMB observation azimuths, namely, the equatorial and southern stripes at rising and setting orientations.

The correction is obtained by projecting the data from a particular band and orientation into a map (after having applied the relative pointing solution described above), and comparing the positions of the bright point sources to catalog positions. These were obtained from the AT20G survey,³⁶ from ATCA. Some of the sources used were extended, but their effect in the result was negligible. The resulting offsets in equatorial coordinates are converted to offsets in the boresight position, and this procedure is iterated to ensure convergence. The boresight offset varies by about 1' over the season, primarily in elevation.

The dominant source of error in the offset correction comes from estimating the centroid positions of the point sources in the maps, before matching them to the catalog positions. This error scales inversely with the square root of the number of point sources available, which was 20 for the 148 GHz band, and with their signal-to-noise ratio (S/N). The uncertainty was found to be 26 for the southern stripe and 53 for the equatorial stripe at 148 GHz. These values were obtained by adding in quadrature the error in the fit from the rising and setting maps.

The telescope pointing is expected to vary slightly due to thermal deformation of the mirror structure. This variation is estimated from observations of Saturn taken at nearly the same azimuth and altitude on different nights. The rms pointing variation of such observations is 43. Since any pixel in the final season maps contains contributions from many different nights, this random variation does not contribute significantly to the pointing uncertainty.

Pointing deviations are significantly higher at dawn when the telescope temperature begins to change more rapidly, showing a trend that repeats every morning. The drift begins nearly 1 hr after sunrise, reaching nearly 50'' an hour later. Data taken more than 1 hr after sunrise are not included in the maps.

On top of this, the altitude was observed to drift by about 20'' over the course of the 120 day season, producing pointing errors that were correlated with right ascension. To remove this trend, linear corrections of 02 day⁻¹ and 015 day⁻¹ were applied to the rising and setting fields, respectively. These corrections have little effect on the maps other than to remove the small, R.A.-dependent pointing offset.

The net pointing uncertainty is thus dominated by the systematic uncertainty in the alignment of the rising and setting maps, along with some residual variation due to temperature-dependent mirror deformation. Comparing the positions of bright point sources in the maps to their catalog positions, we estimate the pointing error in the final maps to be 48, with no preferred orientation (Marriage et al. 2011b). Note that this uncertainty is much smaller than the beam size and would cause errors of less than 1% in the measured flux for point sources.

The detectors are classified in three groups.

• 1.  
  Live detector candidates. Those that have the potential to be used for mapmaking.
• 2.  
  Dark detector candidates. Those that do not couple to the sky (for example a broken pixel) but can be used to diagnose systematics. This includes the subset of the 32 dark detectors which work properly and the subset of defective detectors with working readout circuit.
• 3.  
  Broken detectors. Those defective detectors that cannot be used to probe systematic effects. They include those former live and dark detectors with defective readout circuits, and slow detectors (τ > 10 ms). Many of them must be turned off while observing to prevent them from interfering with other detectors.

Several methods were used to classify detectors into these categories. These include assessing correlation with the atmospheric signal (by far the largest signal), searching for consistently oscillating detectors,³⁷ and finding biasing problems. Table 5 gives a summary of the number of detectors in each of the three groups for each array.

Once the live and dark detector candidates have been identified, a number of possible pathological behaviors may still justify removing part of the data from the reduction pipeline.

Out of the total number of data files acquired, we rejected files for the following reasons and in the following order.

• 1.  
  Files that correspond to planet observations and other calibration or engineering tests.
• 2.  
  Data taken more than 1 hr after sunrise to avoid pointing and beam errors caused by telescope deformations as it was thermally settling.
• 3.  
  Files with fewer than 400 effective detectors, as they were considered likely to be pathological.
• 4.  
  Bad weather: PWV greater than 3.0 mm (transmission below 90%).
• 5.  
  Poor cryogenic performance: detector base temperature more than 7 mK above the nominal temperature or when it changed more than 1 mK within the 15 minute file.
• 6.  
  Poor calibration: if the relative gain dispersion of the detectors was more than 10%.
• 7.  
  Data for which the analysis software failed.

The final amount of data available for analysis is summarized in Table 6.

Detectors are affected by sporadic pathologies. Depending on the kind of pathology, it may be necessary to reject a section or the full length of a detector TOD from a given data file.

The main causes for these pathologies are quantum jumps of the magnetic flux of a SQUID in the readout circuit (V−ϕ jumps), excessive detector noise, conducted noise from oscillating detectors, excessive electromagnetic pickup, and mechanical contamination, which can be optically or thermally coupled into the signal.

The following tests were performed over 15 minute data files to detect these pathologies.

• 1.  
  Drift test. This probes low-frequency deviations of a detector TOD from the atmospheric signal. The data are first low-pass filtered above 50 mHz and calibrated into units of power, as described in Section 7. Then the thermal drift is removed by de-projecting both the dark detector common mode and the housekeeping temperature of the detectors. Finally, the atmosphere signal is removed by de-projecting eight multi-common modes, as explained in the Appendix. The drift error is defined as the standard deviation of the residual data. Detectors in 148 GHz are cut as outliers if their drift error is higher than 0.35 fW.
• 2.  
  Correlation test. The drift test is complemented by finding the correlation between the detector drift and array common mode drift for TOD frequencies below 50 mHz. Detectors that correlate less than 98% are excluded.
• 3.  
  Gain test. The drift test assesses the shapes of the TODs, but not their amplitudes. The relative detector gain can be quantified by the factor that best fits the atmospheric drift to every single TOD at frequencies below 50 mHz. For this the atmospheric drift is estimated as the array common mode after removing the thermal contamination. Detectors are cut whenever their gains differ from the mean by more than 15%.
• 4.  
  Mid-frequency noise test. Some pathologies, mainly associated with mechanical contamination, become prominent at frequencies between 0.3 and 1.0 Hz, below the "science band." To isolate pathological detectors, we band-pass filter the data below and above those frequencies, de-project the array common mode (also filtered), and obtain the standard deviation of the residuals, which we call mid-frequency error. Detectors are cut whenever their mid-frequency error is more than 2.5 times the median value for the data file.
• 5.  
  High-frequency noise tests. At frequencies above 3 Hz the data start being dominated by detector noise, which is chiefly Gaussian. Non-Gaussianities above this frequency are mostly caused by electromagnetic pickup, conducted electrical noise or opto-mechanical perturbations. To isolate them, we high-pass filter the data below 5 Hz and then test for non-Gaussianities by computing the standard deviation, skewness, and kurtosis. This is done within sections of the TOD of one scan period, with a characteristic length of 10.2 s. Using the transformation proposed by D'Agostino & Belanger (1990), the last two statistics yield a normal distribution in the case of Gaussian noise, which is verified for most of our data. Outliers are identified by how much they deviate from the mean. Only sections of the TOD with noise rms lower than 0.9 fW are kept.
• 6.  
  Scan test. Motion defects and encoder errors are also detected and cut, affecting sub-sections of the TOD. This is done by examining the azimuth time stream and searching for interruptions in the scan pattern.
• 7.  
  Glitch test. The data are also affected by spike-like glitches, for example from cosmic rays. Spikes larger than 10 times the noise rms are cut, including a 0.5 s buffer on either side. Also, whenever two such cuts are separated by less than 5 s they are stitched together into a single cut. If more that five glitches are found in a single detector TOD, then the whole detector TOD is cut.
• 8.  
  Calibration bolometer test. During 2008 observations a calibration bolometer was used to load the detectors with radiation of roughly 400 mK every 24 minutes, each event lasting for 1.3 s. For these, a window of nearly 3 s is excised around the event.³⁸ We included this within the "Glitch" cuts in Table 5.

As a general rule, if more than 20% of the detector TOD would be cut, then the full detector TOD is cut instead.

Dark detectors were selected in an analogous but simplified way. In this case, only full detector TOD cuts were performed. The cut criteria were the drift error of the dark detectors, the gain with respect to the dark common mode,³⁹ and the noise rms, all given in raw data units.

After applying all the selection criteria given above, and considering only the data files available for analysis, the average number of effective detectors was 593 ± 95 in the 2008 season at 148 GHz. The error here is the standard deviation over different data files.

Table 5 shows a summary with the number of detectors in each of the three detector groups for the 148 GHz array, as well as the cut contribution from the eight main selection criteria: drift, gain, mid-frequency noise, noise rms, noise skewness, noise kurtosis, and partial cuts with more than 20% of the TOD cut. The bottom line is the average number of effective detectors in the season. Figure 8 is a diagram of the 148 GHz array showing the fraction of the time that each detector is uncut. Figure 9 shows the daily number of the effective detectors during the 2008 season, along with the weather conditions indicated by the PWV.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

On top of these cuts, a small fraction of the remaining data is cut during the mapmaking process, as described in the following section.

To make maps from ACT data, we solve for the best-fit sky given the noise in the data. In particular, we find the sky map that minimizes χ² given a model for the noise, and a model for what the data should look like:

$\begin{eqnarray} \mathbf {d} &=& \mathbf {M}\mathbf {x}+\mathbf {n}. \end{eqnarray} \tag{ 10 }$

Here $\mathbf {x}$ is the model for which we wish to solve, $\mathbf {M}$ describes how the data depend on the model, and $\mathbf {n}$ is the particular realization of the noise in the ACT data. Traditionally, $\mathbf {x}$ is a vector whose components are the sky map pixels and $\mathbf {M}$ is the pointing matrix. In its simplest conceivable form, each data point sees a single pixel in the map, so each row of $\mathbf {M}$ (corresponding to a single data point) has a single 1 in the column corresponding to the map pixel at which it was pointed. However, our model for the data, $\mathbf {x}$, can contain many contributions in addition to the sky map: components include atmospheric noise, correlated signals in the data, and missing (cut) data samples as discussed below. The mapping formalism easily generalizes to cover multiple components as long as the data depend on them linearly:

$\begin{eqnarray} \mathbf {M} &=& \left[\begin{array}{@{}cccc@{}}\mathbf {M}_1 & \mathbf {M}_2 & \cdots & \mathbf {M}_c\end{array}\right], \end{eqnarray} \tag{ 11 }$

$\begin{eqnarray} \mathbf {x} &=& \left[\begin{array}{@{}c@{}} \mathbf {x}_1 \\ \mathbf {x}_2 \\ \vdots \\ \mathbf {x}_c\end{array}\right]. \hspace{60pt} \end{eqnarray} \tag{ 12 }$

With an estimate of the noise covariance $\mathbf {N} \equiv \left<\mathbf {n} \mathbf {n}^T \right>$, which is an n_sample by n_sample matrix (n_sample ≃ 10⁹), we wish to find the model that maximizes the likelihood function

$\begin{equation} \mathcal {L} = \exp \big[{-}\case{1}{2}(\mathbf {d} - \mathbf {M}\mathbf {x})^T \mathbf {N}^{-1}(\mathbf {d} - \mathbf {M}\mathbf {x})\big], \end{equation} \tag{ 13 }$

where $\mathbf {d}$ is a vector containing all the data samples.

The standard linear least-squares solution is

$\begin{equation} \mathbf {M}^T\mathbf {N}^{-1}\mathbf {M}\mathbf {x} = \mathbf {M}^T\mathbf {N}^{-1}\,\mathbf {d}. \end{equation} \tag{ 14 }$

The matrix $\mathbf {M}^T\mathbf {N}^{-1}\mathbf {M}$ is too big to be practicably inverted directly (we typically have 10⁷ map pixels), so we instead iteratively solve the least-squares equation for $\mathbf {x}$ using a Preconditioned Conjugate Gradient (PCG) scheme (Wright et al. 1996; Hinshaw et al. 2003; Press et al. 2007). Preconditioning involves introducing $\mathbf {P}$, an approximate inverse of $\mathbf {M}^T \mathbf {N}^{-1} \mathbf {M}$, in order to speed up convergence of classic Conjugate Gradient, and solving the better-conditioned system:

$\begin{equation} \mathbf {P}\,\mathbf {M}^T\mathbf {N}^{-1}\mathbf {M}x = \mathbf {P}\,\mathbf {M}^T\mathbf {N}^{-1}\,\mathbf {d}. \end{equation} \tag{ 15 }$

To map ACT data, we use the mapmaking code Ninkasi and run on the Scinet General Purpose Cluster (GPC; Loken et al. 2010). The mapmaking algorithm for the maps in this release is improved from that used in previous papers (Fowler et al. 2010; Marriage et al. 2011b, 2011a; Dunkley et al. 2011; Das et al. 2011a; Hajian et al. 2011; Sehgal et al. 2011) in four ways: (1) rather than solving for detector-correlated noise (including atmospheric noise) explicitly, we put the correlations in the noise matrix, producing a substantial noise improvement in the damping tail; (2) we explicitly solve for cut data (time-stream gaps), which made the problem symmetric; (3) we subtract models for the point sources in the ACT maps directly from detector time streams, recovering part of the power that otherwise was biased low by the mapper; and (4) we re-estimate the noise after subtracting an initial map to remove signal-induced bias in the noise estimation.

The mapping is done in a cylindrical equal-area projection with a standard latitude of δ = −535 and pixels of 30'' × 30'', roughly one-third of the beam FWHM.

In addition to data selection, calibration, and pointing, there are a few other pre-processing steps that are done to the data before solving for the maps.

We remove the median of each detector time stream for each 15 minute period, and subtract a single array-wide slope across the period so that the ends of the time streams roughly line up. We do this to reduce ringing in Fourier transforms and to facilitate searching for correlations among the detectors. We linearly interpolate across gaps in the data, such as those arising from cosmic-ray hits (again, to reduce Fourier artifacts: the data in the cuts are otherwise not used). We deconvolve the time streams by the anti-alias filter (described in Section 5) and the detector time constants (discussed in Section 9).

Next, the non-optical contamination signals are reduced using the dark detectors. This is done as follows: we take the time streams from each dark detector and subtract a mean and a slope. We then high-pass filter the dark time streams, filtering out any signal below 5 Hz. Most of the remaining signal is in only a few independent modes. We take those corresponding to the seven largest eigenvalues of the resulting covariance matrix and do a linear least-squares fit to the live detector data, which had been similarly processed (having removed the mean and slope, and high-pass filtered). Finally, we subtract the reconstructed fit from the data. We can do this because the modes subtracted are not correlated to the sky signal.

In Fowler et al. (2010), we downsampled the time streams from 400 Hz to 200 Hz, using a time-domain triangular kernel of the form [0.25 0.50 0.25]. The ACT signal band (nearly 30 Hz according to Section 3) is well away from the downsampled frequency limit. We find that downsampling does affect the raw power spectrum at up to several percent for ℓ  >   ∼ 5000. This should be mostly corrected for when using a beam measured using maps derived from downsampled data. However, we also find that the S/N on point sources is 2% higher for non-downsampled maps, triggered by the anti-aliasing filter used for downsampling. Consequently, we do not downsample the time streams in the results presented here.

Also in contrast to Fowler et al. (2010), we carry out one further step in the time domain. With the 30'' pixels somewhat undersampling the ACT 137 beam, we find that to get percent-level accuracy in source fluxes, we must deal with the bulk of the source flux directly in the time streams. To do this, we make first-pass maps in which sources are found and their fluxes estimated. Simulations show that these fluxes are typically accurate to a few percent, with source flux systematically underestimated by approximately 4%. We then take these source fluxes and model them directly in the time streams using the full ACT beam, and subtract the model from the time streams. We do a simple source-only projection into a map, and add this to the (mostly) source-free maximum likelihood map. We find that, in simulations, this recovers mean source fluxes to 1% accuracy, at which point it is subdominant to the calibration uncertainty. We do no additional filtering in the time domain.

The noise structure is modeled in both frequency and time domains. We include, as a term in the noise, a time domain windowing of the first and last 20 s of each time stream of the form (1 − cos x)/2, again to reduce Fourier ringing. We then search for correlations across the data by examining their covariance matrices, split at 4 Hz (ℓ ≃ 1500). We find all eigenvectors in the low-frequency covariance matrix with corresponding eigenvalues greater than 3.5² times the median (the eigenvalues of the data covariance matrix Δ^TΔ are the squares of the corresponding data singular values, so this corresponds to an amplitude of 3.5 times the median in the time streams). We then project those eigenvectors out of the high-frequency covariance matrix, and again find all eigenvectors with eigenvalues larger than 3.5² times the median in the remaining high-frequency matrix. We find that the maps are not particularly sensitive to the exact threshold values chosen. This procedure typically finds 15–20 low-frequency eigenvectors and one or two high-frequency ones. The eigenvectors are response patterns of correlations across the array, and usually correspond to things like a common mode, gradients across the array, and the row-correlated noise. They provide the linear combination of detector time streams needed to produce a correlated mode, represented by the vector $\mathbf {\hat{v}}$ in Equation (A1). This projection corresponds to removing around 20 modes out of the approximately 600 available.

With the shapes of the array correlations in hand, we can use them to complete the description of the noise. We solve for the correlated modes corresponding to the eigenvectors and subtract them from the time streams. We then Fourier transform the (de-correlated) detector time streams and the correlated modes, and model them using frequency bins. In each frequency bin, we find the average power in each detector time stream and in each correlated mode. If the detector noises are denoted by the diagonal matrix $\mathbf {N}_{D}$, the detector covariance eigenvectors by the n_detector by n_mode matrix $\mathbf {V}$ and their corresponding bin-wise noise powers by $\mathbf {N}_{V}$, then the bin-wise Fourier-space noise is simply $\mathbf {N}_{f} \equiv \mathbf {N}_{D}+\mathbf {V}\,\mathbf {N}_{V}\,\mathbf {V}^T$. Here $\mathbf {N}_f$ is an n_detector by n_detector matrix acting on a single frequency bin. This matrix can be quickly inverted using the Sherman–Woodbury formula (Duncan 1944; see Hager 1989 for a review). If we denote the time-domain edge tapering operator by $\mathbf {W}$ and the Fourier transform operator by $\mathbf {F}$, then our entire noise inverse becomes

$\begin{equation} \mathbf {N}^{-1} \equiv \mathbf {W}^{T}\,\mathbf {F}^{T}\,\mathbf {N}_{F}^{-1}\,\mathbf {F}\,\mathbf {W}, \end{equation} \tag{ 16 }$

where the operator $\mathbf {N}_{F}^{-1}$ is conformed by all the $\mathbf {N}_f^{-1}$ such that it acts on every frequency bin in the Fourier transform. In the previous equation, every operator can be represented by an n_sample by n_sample matrix, in which detector operations like $\mathbf {N}_f$ are expanded such that all samples within a single detector are treated in the same way. To minimize our sensitivity to any low-frequency non-Gaussian component of the atmospheric noise, we taper the Fourier weights at frequencies below 0.5 Hz, with the weight explicitly set to zero below 0.25 Hz. We note that the noise in the mapmaking equation (Equation (15)) can be interpreted as a set of weights. The map is optimal if the weights are perfect, but the map remains unbiased as long as identical weights are used for the left and right sides of the equation, and the weights are uncorrelated with the sky signal. This is the essential part of the method. Time-stream filters are by design biased and must be accounted for with simulations. On the other hand, in our treatment, modes in the map that might cause problems are de-weighted as opposed to being filtered, producing an unbiased solution through careful consideration of the noise structure of the data.

One further difference remains between these maps and those used in Fowler et al. (2010). Here, for every single time-stream sample cut, we fit for its (unknown) amplitude as part of the map solution, as suggested by, for example, Patanchon et al. (2008). In addition to cuts from cosmic rays and the like, we explicitly cut the first and last second of each TOD to give the mapping algorithm the freedom to match the TOD beginning and end as smoothly as possible. Since these data are already highly downweighted by the time-domain taper, the additional amount of data loss is negligible.

With these pre-processing steps, form of the noise, and form of the solution specified, we must then actually solve for the map. We use the classic PCG algorithm, using the hit-count map as a Jacobi (diagonal) preconditioner for the sky map part of the solution. To recover the scales around ℓ ≃ 200 and below, we need a few hundred PCG iterations. The mapmaking is quite CPU-intensive, with each iteration taking about a minute on 1600 2.53 GHz Nehalem cores on the Scinet GPC, or a full run to 1000 iterations taking about a day of wall clock time, or 5 yr of CPU time.

We note that care must be taken when estimating the noise to make sure it does not bias the map. Consider the following thought experiment. Two detectors are scanned across a source, and their weights are estimated from the internal scatters of their time streams. In general, noise in the detectors will lead to one of them observing a lower flux from the source. If the detector weights are then measured and the source is not removed from the time streams, the detector which happened to measure the lower flux will usually have a lower variance in its data, and will therefore on average receive more weight than the other detector. When the detectors are combined, the more heavily weighted low-flux measurement will lead to the source flux being systematically biased low. The magnitude of the bias is set by the S/N in the individual time streams, not the final S/N. The bias can be mitigated if an estimate of the signal is removed from the time streams before noise estimation. Since the CMB is highly subdominant to the noise in ACT time streams, this effect is small, but simulations show that it is not negligible. Therefore, first we make an initial map using the full data set and the noise measured directly from the time streams, then we subtract that map from the time streams to reduce the sky signal in them before estimating the noise for a second time, and finally we use this improved estimate of the noise to make the final maps.

ACT must take a bit more care than all-sky experiments because of the large change in sensitivity across the map. In particular, near the map edges, only a small amount of data contributes to the map, and so removing the raw map would lead to an artificial reduction in the noise of the data near the edges, which would lead to an artificial increase in their weight. To prevent this from happening, the starter maps used are first filtered and apodized. The apodization is generated by first rescaling the hit count map to the 0–1.0 range, then setting all values higher than a threshold to unity, and finally smoothing the resulting weight map with a Gaussian window. Furthermore, from the starter data map we filter out scales larger than ℓ = 300 where atmospheric noise is very large, and scales smaller than ℓ = 3000 where the map is highly noise dominated and simulations show that the bias effect is truly negligible. Finally, we generate the starter map by multiplying the filtered data map with the apodization window. To test this procedure, we carry out the two-step procedure on both the real data and the real data with a simulated map injected into it. We compare the difference of these two sets of maps to the original simulation and show the resulting transfer function in Figure 10. At the low-pass scale of ℓ = 300, the maps are unbiased to 0.5%, and rapidly become unbiased to better than a part in 10³ on smaller scales. Moreover, we find consistency of the ACT×WMAP map cross spectra with the ACT auto-spectrum and the WMAP binned power spectrum for scales in the range 300 < ℓ < 1000, as shown in Hajian et al. (2011).

Zoom In Zoom Out Reset image size

The resulting maps cover an area of 845.6 deg² on the sky, ranging between 20^h43^m and 7^h53^m in right ascension and −481 and −572 in declination. A total of five maps were made: one using the full data set, plus four partial versions using independent subsets of the data for noise estimation purposes. Each map is accompanied with a hit map with the number of data samples that fell on every pixel of the corresponding map.

The noise in the maps is strongly scale dependent and fairly distributed according to the hit counts. In general, noise at smaller scales is mostly uncorrelated between pixels and proportional to the square root of the number of hits, while at larger scales correlations produced by the atmosphere and systematics become more evident. Correlated noise is better behaved in Fourier domain where its covariance is approximately diagonal, so it can be represented as a power spectrum. Examples of noise analysis done to these maps can be found in Fowler et al. (2010), Marriage et al. (2011b), and Das et al. (2011a).

These maps are an overall improvement over those used in Fowler et al. (2010), Marriage et al. (2011b, 2011a), Dunkley et al. (2011), Das et al. (2011a), Hajian et al. (2011), and Sehgal et al. (2011), marked particularly by lower noise on angular scales between 500 < ℓ < 2000, reducing the variance by up to a factor of two or three. However, the original maps were already nearly sample limited on these scales, so errors on the power spectra are only modestly reduced. These changes are due to slightly different data selection cuts, improvement in the noise treatment during the mapping process, and improved per-detector calibration.

Figure 11 shows the final 2008 southern map, overlaid by contours of iso-sensitivity (35, 50, and 65 μK arcminute), as different regions of the map have different integration times. These sensitivities were computed using the number of samples per unit area and the NET given in Table 3. We run for 2000 PCG iterations on this map to ensure that it is truly converged, though as mentioned we find that all science scales are converged in far fewer PCG iterations.

Zoom In Zoom Out Reset image size

Figure 12 shows the power spectrum computed for this map, compared to previous versions already published in Fowler et al. (2010) and Das et al. (2011a). The uncertainty in the spectrum reveals an improvement in the quality of the map.

Zoom In Zoom Out Reset image size

The parameters obtained using a fit to the spectrum released with this paper are slightly shifted with respect to the values presented in Dunkley et al. (2011). The largest shifts in the primary parameters were for n_s and Ω_b h² at Δσ = 0.3 and Θ_A with Δσ = 0.7, although the change in the inferred dark energy density changed by a negligible 0.03σ. The other parameters changed by less than 0.07σ. The shift in Θ_A is dominated by changes around the third peak of the power spectrum, caused by the new calibration, cuts, and pointing used in these maps. The secondary parameters were affected by the new mapping procedure, decreasing the Poisson point source level from 13.7 μK² to 12.5 μK². The upper 95% confidence limits of the correlated point source amplitude and SZ amplitude are instead unchanged.

Also, a cross-correlation analysis to the BLAST maps (http://blastexperiment.info/) shows correlations detected at a 25σ level, implying a detection of emission by radio and dusty star-forming galaxies at high ℓ, as described in Hajian et al. (2012).

Finally, the flux of sources as estimated in the data release is approximately 5% higher than those reported in Marriage et al. (2011b), due to the new source treatment in the mapmaking and a change in the beam solid angle. The best representation of the beam shape, as well as the window function for power spectrum analysis, is released together with the maps.

A.0.1. Dark Mode Removal