Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Scuba half-degree extragalactic survey : data reduction by Iterative deconvolution Lepage, Kyle Quentin 2004

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2004-0253.pdf [ 13.43MB ]
Metadata
JSON: 831-1.0091468.json
JSON-LD: 831-1.0091468-ld.json
RDF/XML (Pretty): 831-1.0091468-rdf.xml
RDF/JSON: 831-1.0091468-rdf.json
Turtle: 831-1.0091468-turtle.txt
N-Triples: 831-1.0091468-rdf-ntriples.txt
Original Record: 831-1.0091468-source.json
Full Text
831-1.0091468-fulltext.txt
Citation
831-1.0091468.ris

Full Text

Scuba Half-Degree Extragalactic Survey: Data Reduction by Iterative Deconvolution by Kyle Quentin Lepage BASc, The University of British Columbia, 2000 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in The Faculty of Graduate Studies (Department of Physics and Astronomy) We accept this thesis a$ conforming to the required /tjandard THE UNIVERSITY OF BRITISH COLUMBIA April 21, 2004 © Kyle Quentin Lepage, 2004 Library Authorization In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Name of Author (please print) Date (dd/mm/yyyy) Title of Thesis: , . , L „ C y U l f - D ^ r ^ E * 4ra ^.!« r-l,V ^>ceo KS-O \ r^ > Degree: H ^ - U r <v£ S p ' i ^ Y e a r : ^&c>h Department of \>uy*<_s. & ^ V ^ c ^ - v The University of British Columbia Vancouver, BC Canada Abstract n Abstract A generalization of an iterative algorithm used to reduce large data sets consisting of observations of the cosmic microwave background is presented. In particular, the al-gorithm is modified to reduce data sets obtained by arbitrary observations. The gen-eralized algorithm is applied, in simulation, to a data set consisting of complicated observations taken during a Scuba HAlf Degree Extragalactic Survey (SHADES). The algorithm is found to be unstable. The generalized time-order algorithm is reformulated as a difference equation whose solution gives the Nth iteration reduced data map in analytical form. Within the analytical framework, the convergence properties of the algorithm are determined, a stable, known work-around explained, and the converged map shown to be optimal in the least-squares sense. The known work-around, a modification of the generalized iterative algorithm, trades stability for ideal behaviour at the edge of the reduced data map. Using simulations of a sky populated by a realistic source count model, it is shown that the deleterious effect of the modified algorithm towards the edge of the map is negligible as compared to the statistical noise for most of the half-degree square map. Simulations of the modified algorithm applied to data taken with the SHADES measurement strategy with Gaussian noise added, reveal the presence of streaking and large scale fluctuations in the converged map. The streaking which is demon-strated to be affected only by the measurement strategy, is indicative of a band-diagonal pixel-pixel covariance matrix, and is absent in noiseless simulations. The large scale fluctuations are indicative of non-negligible, far off-diagonal terms in the pixel-pixel covariance matrix, and are also absent in noiseless simulations. The ef-fect of different observation strategies on the streaking is investigated by simulation. The presence of the amplified large scale noise is explained by linear systems theory. Abstract iii The determination of the optimal observation strategy is formulated as the task of determining the linear combination of absolute power measurements which will yield a desired pixel-pixel covariance matrix. This problem is discussed. An un-deconvolved map of real SHADES data is produced. Sources are detected in the un-deconvolved map using Bayesian techniques. Contents iv Contents Abstract ii Contents iv List of Tables vii List of Figures ix Acknowledgements xviii I Thesis 1 1 The Scuba Half-Degree Extragalactic Survey 2 1.1 Introduction To The Submillimeter Galaxy Population 3 1.2 What Is The Cosmic History Of Massive Dust-Enshrouded Star For-mation Activity? 4 1.3 Are SCUBA Sources The Progenitors Of Present-Day Massively El-liptical Galaxies? 5 1.4 What Fraction Of SCUBA Sources Harbour A Dust-Obscured Active Galactic Nucleus? 5 1.5 Legacy Value 6 1.5.1 Data Compression 6 1.5.2 Data Correlation: Integrated Sachs-Wolfe Effect 6 2 Observing The Scuba HAlf-Degree Extragalactic Survey 8 2.1 Submillimeter Common User Bolometer Array . 8 2.2 Data Acquisition 9 Contents v 2.3 SCUBA/JCMT Measurement Noise 11 2.3.1 Brief Overview Of Noise 11 2.3.2 SCUBA Noise 13 2.4 Data Preprocessing 14 2.4.1 Calibration 14 2.4.2 Nuisance Signals And Outlier Rejection 18 2.5 Data Reduction 21 2.5.1 Interpolative Reduction 23 2.5.2 Emerson Reconstruction 27 2.5.3 Direct Inversion 28 2.5.4 Iterative Inversion 29 3 Iterative Reduction of S H A D E S Data 31 3.1 Wright's Algorithm Generalized 31 3.1.1 Algorithm Checks 32 3.2 Generalized Algorithm Simulations 33 3.2.1 SHADES Simulation 34 3.2.2 Small Field Simulations 37 3.2.3 Minimum Comparable Size 43 3.3 Iterative Inversion As Newton-Raphson Root Finding 45 3.4 Formulation of Iterative Inversion as a Difference Equation 46 3.4.1 Formulation 46 3.4.2 Solution 48 3.4.3 Solution Check #1 49 3.4.4 Solution Check #2 50 3.4.5 Solution Check #3 50 3.4.6 Difference Equation Solution Simulation 52 3.4.7 Stability 54 3.4.8 Noise Properties 57 3.5 Algorithm Stabilization 59 3.5.1 Unstable Eigenvector Removal 59 Contents vi 3.5.2 Boundary Pixel Modification 61 3.5.3 Measurement And Pixel Update Decoupling 62 4 Central Update Algori thm 71 4.1 Edge Effects 71 4.2 Gaussian Noise Simulation 81 4.3 More Realistic Noise Simulation 87 4.4 Chopping Effects On Pixel-Pixel Correlation 92 4.5 Chapter Results 92 5 S H A D E S Subaru X M M - D e e p Field 95 5.1 Data Reduction 95 5.2 Source Detection 98 5.3 Deconvolved Map Source Detection 106 6 Discussion 113 6.1 Main Results 113 6.2 Follow-Up Observation 114 6.3 Measurement Strategy 116 6.4 Optimal Measurement Strategy 118 6.5 Noise &; Correlation Removal 119 A Jordan Form 121 Bibliography 123 List of Tables vii List of Tables 3.1 Correlation Mat r ix : For a seven pixel, seven measurement, single-difference simulation 40 3.2 Correlation Mat r ix : Seven pixel, thirteen measurement, double-difference simulation 42 3.3 Correlation Mat r ix : Seven pixel, fifteen measurement, double-difference simulation 42 3.4 Eigenvalues of the characteristic matrix for the seven pixel, five double-difference measurement simulation 52 3.5 Eigenvectors of the characteristic matrix in the seven pixel, five double-difference measurement simulation. Each row represents a separate eigenvector. The first row is the eigenvector corresponding to the left-most eigenvalue in Table 3.4. For clarity, the eigenvec-tor components have been rounded to three significant figures. The eigenvectors are not linearly independent 52 3.6 Eigenvalues of the characteristic matrix for the seven pixel, 15 double-difference measurement simulation 53 3.7 Eigenvectors of the characteristic matrix for the seven pixel, 15 double-difference measurement simulation. Each row is a separate eigenvector; the eigenvector in the top row corresponds to the left most eigenvalue in Table 3.6 53 3.8 Eigenvalues of the characteristic matrix for the six pixel, four double-difference measurement simulation with two fixed pixels. There is an absolute eigenvalue greater than one 56 List of Tables viii 3.9 Eigenvalues of the characteristic matrix in the simulation with six pixels, two double-difference measurements and two absolute mea-surements. There is an absolute eigenvalue greater than one 56 3.10 Eigenvalues of the characteristic matrix for the six pixel, four double-difference measurements and two absolute measurements simulation. There is an absolute eigenvalue greater than one 57 3.11 Eigenvalues of the characteristic matrix for the six pixel, four double-difference measurement and two absolute measurement. There are no absolute eigenvalues greater than one 58 3.12 Eigenvalues of the characteristic matrix for the seven pixel, five double-difference measurement simulation 61 3.13 Characteristic eigenvalues for the edge-pixel enlarged pointing matrix 62 3.14 Characteristic matrix eigenvalues for the central update algorithm in the seven pixel, five double-difference measurement simulation. All of the absolute eigenvalues are less than or equal to one 63 5.1 Subaru X M M Deep Field Source Detection List - Part 1 110 5.2 Subaru X M M Deep Field Source Detection List - Part II .... I l l 5.3 Subaru X M M Deep Field Source Detection List - Part III ... . 112 List of Figures ix List of Figures 2.1 SCUBA 850 micrometer bolometer array. The plot is obtained from the positions of actual data on the sky. Note the spacing between equal distances along the two axes are different, the height and width of the array are equal 9 2.2 Tripos Coverage 11 2.3 Eight representative scuba bolometer timestreams and spectra for the Subaru XMM Deep Field after flatfielding, and extinction correction. The timestreams are non-stationary and strongly correlated by at-mospheric emission. The timestream for bolometer 21 extends past the range of the y-axis. The standard deviation of the bolometer timestreams varies substantially from bolometer to bolometer. The magnitude of the Fourier transform of the bolometer timestreams are shown in the lower half of the plot. Note the flat spectra and the spikes on bolometers 17, 19, 20, 21, and possibly 23 and 24 16 2.4 Flux conversion factor as a function of observation file number. The observation files are not in time order. Note the considerable fluctu-ation of the FCF with observation number. Only calibrations with CRL618, Uranus and Mars are used. The chop throw of the calibra-tions is ignored 17 List of Figures x 2.5 Correspondence between sources detected at greater than 3a signif-icance from calibrated data and from uncalibrated data. The cali-brated data uses Method 3 in Section 2.4.1 to compute flux conver-sion factors. The source detection is described in Section 5.2, but differs in that the detected sources have not been iteratively cleaned from the map 19 2.6 Eight bolometer timestreams and their corresponding Fourier trans-form magnitudes for the Subaru XMM Deep Field after preprocessing. Note the lack of outliers in the preprocessed timestreams as well as the reduced rms when compared to the timestreams in Figure 2.3. The Fourier transform magnitudes are flatter than their non-preprocessed counter-parts. (Figure 2.3) 22 2.7 A measurement map of simulated measurements taken with the single pointing SHADES observation strategy of a fake sky with a single, centred source. Note the presence of negative sidelobes centred upon the six chop throws surrounding the source. Each square pixel in this map is three arcseconds square 25 3.1 Residual maps of a simulation of the generalized iterative algorithm. The measurements are performed on a statistically plausible fake sky using the full SHADES measurement strategy. The measurements are noiseless. High spatial frequency structure is removed from the residuals in the early iterations. In later iterations the unstable spa-tial frequencies begin to grow. The fake sky map saturates at pixel values of zero and ten. Pixels are three arcseconds square 35 3.2 Residual maps for iterations 1 through 6 and 25 in a simulation of the generalized iterative algorithm applied to measurements taken of a single source. The measurements comprise a single pointing in the SHADES measurement strategy. 36 List of Figures xi 3.3 Top: Fake sky map. Middle: Fake sky map estimates for iteration 1, 11 and 21. Bottom row: Residual maps for iteration 1, 11 and 21. The colour bar applies to the residual maps only 64 3.4 Residual maps for iteration 1, 15 and 25. Note the reduction of high-frequency spatial structure between iteration 1 and 15 and the re-appearance of high-frequency structure in iteration 25. The algorithm slowly diverges 65 3.5 Single-Difference: 7 Pixel Simulation, 6 Measurements 66 3.6 Single-Difference: 7 Pixel Simulation, 7 Measurements 67 3.7 Double-Difference: 7 Pixel Simulation, 5 Measurements 68 3.8 Double-Difference: 7 Pixel Simulation, 13 Measurements 69 3.9 Double-Difference: 7 Pixel Simulation, 15 Measurements 70 4.1 Standard deviation of the residual map pixels as a function of itera-tion for a central update algorithm data reduction simulation. The simulated data consist of measurements of a single, lOmJy, centred source, using the full SHADES measurement strategy. The compu-tation neglects pixels within 150" from the map edge. The standard deviation of the residual map pixels is 0.0008 mJy at iteration 101. At iteration 1001 the standard deviation of the residual map pixels is less than 0.0001 mJy 72 4.2 Residual maps estimated on various iterations of a simulation of the central update algorithm using the full SHADES measurement strat-egy. The fake sky map is a single source, centred in the map, with a peak flux of ten. The pixel values are in mJy and the contrast has been stretched on the lower three panels. The first iteration map is a measurement map. Note the convergence and the absence of a visible effect caused by the edges of the map. The maximum absolute pixel value of the residual map on the 75th iteration is less than a tenth of a percent of the source amplitude 73 List of Figures xii 4.3 Standard deviation of the residual map pixels as a function of iter-ation during a central update algorithm simulation. The simulated measurements are of an edge only, statistically plausible fake sky map and are noiseless. Pixels within 150" of the map edge are not used in the computation. At iteration 401, the residual map pixel standard deviation is within 0.02 mjy of the converged value of 0.07 mjy. The steps in the curve between iterations 401 and 801 are due to limited precision when reading and writing data to file. Note the maximum residual map pixel standard deviation of 0.185 mjy near iteration 60. At iteration one, the initial map estimate is the null map, and all pixels used in the computation of the standard deviation are equal to the fake sky map. Thus the standard deviation of the residual map pixels is zero, and rises with iteration as pixels inside the periphery are modified due to flux in the chop positions of the measurements centred upon pixels at the edge of the periphery. 75 4.4 An edge-only fake sky map and residual maps for a noiseless, full SHADES simulation of the central update algorithm. The fake sky map is populated with sources drawn from a realistic source count model. The edges are 150" wide. The pixel values are in mjy. Note the growth of the edge error. The maximum residual pixel standard deviation occurs around iteration 60. Note how the mean of the residual map pixels grows increasingly negative, reducing the stan-dard deviation. At iteration 101, the absolute pixel values are less than two and a half percent of a 10 mjy source, 360" from the map edge. The white lines are spaced 200 pixels apart, each pixel is 3" wide. 76 4.5 The 801s* residual map for the statistically plausible, edge only, cen-tral update algorithm simulation. The residual due to the edge effect is less than 0.15 mjy, 150" from the map edge 77 List of Figures xiii 4.6 Residual map standard deviation as a function of iteration during a simulation of data reduction using the central update algorithm. The data consist of measurements of a statistically plausible fake sky ob-tained with the full SHADES measurement strategy. Convergence to within 0.005 mJy of the converged standard deviation occurs by itera-tion 50. The converged standard deviation of 0.068 mJy is 0.002 mJy from the the converged standard deviation of the residual pixels in the edge-only simulation (Figure 4.3). The minimum of 0.0655 mJy in the standard deviation of the residual pixels occurs near iteration 40 79 4.7 Residual maps for a full SHADES simulation of the central update algorithm. The fake sky map is populated with sources drawn from a realistic source count model. Between iterations 50 and 90 the variance of the residual pixels reduces in the centre of the map and slightly increases toward the map edges. The pixel values are in mJy. Notice the high spatial frequency residuals, visible at iteration one, are largely gone by iteration 40 80 4.8 Residual map average and standard deviation as a function of itera-tion during a simulation of data reduction using the central update algorithm. Gaussian random noise is added to simulated measure-ments of a fake sky containing no flux. At iteration 100, the resid-ual pixel standard deviation is 1.95 mJy higher than the expected ^measurement/ViV The central-update algorithm converges when re-ducing measurements with Gaussian additive noise 82 4.9 Measurements with Gaussian noise are iteratively reduced using the central update algorithm. The maps are, clockwise from upper left, the 1st iteration map, the 101st iteration map, the power spectrum of the 1st iteration map and the power spectrum of the 101st and 1st iteration map. Note the structure that appears in the 101st iteration maps 83 List of Figures xiv 4.10 No source streaking, other than the existing cross-hatch pattern, is apparent, in a simulated reduction of measurements of four large sources containing Gaussian additive noise using the central update algorithm. The pixel values are in mjy 85 4.11 Measurements of a statistically plausible fake sky with additive Gaus-sian noise are iteratively reduced using the central update algorithm. The maps are, clockwise from upper left, the 1st iteration map, the 101st iteration map, the power spectrum of the 1st iteration map and the power spectrum of the 101st iteration map. The presence of sources appears as low frequency power in the 1st iteration power spectrum map 86 4.12 First four plots are of the randomly extended noise sequences for bolometers 1, 2, 9, and 35. The bolometer selection is arbitrary. Last four plots are the corresponding Fourier transform magnitudes. Note the lack of j noise in the spectrums. There is a bump at approx. 0.007 Hz in the bolometer 35 spectrum, possibly the Fourier transform spike described in chapter 2 88 4.13 Residual map standard deviation as a function of iteration in a central update algorithm simulation with a statistically plausible fake sky map and the full SHADES measurement strategy 89 4.14 Clockwise from upper left: Sky estimate maps for iteration one, 101, power spectrum map for iteration 101, power spectrum map for it-eration one. Note the similarity between these maps and the maps produced in simulations with Gaussian additive noise. (Figure 4.12) The off-diagonal terms of the timestream noise-covariance matrix play a minor role in the pixel-pixel correlations present in the 101st itera-tion map. The pixel-pixel correlations are predominantly due to the off-diagonal terms in the PTP matrix for the SHADES measurement strategy. The PTP matrix was introduced in chapter two 91 List of Figures xv 4.15 Simulation of data reduced using the central update algorithm. The data are simulated measurements of a statistically plausible fake sky with blank field extended noise samples. One half of the chops are rotated 45 degrees from the east-west and north-south chop direc-tions. Much of the cross-hatching visibile in the 101st iteration map in Figure 4.14 has disappeared. The strong cross-shape in the 101st power spectrum map in Figure 4.14 has been reduced. Large scale structure with half-wavelengths around 600" seem unaffected 94 5.1 The Subaru XMM-Deep Field measurement map. Units are in Jy, pixels are 3" square. The horizontal and vertical axes are specified in units of pixels 96 5.2 Subaru XMM-Deep Field noise map and hit map, top and bottom, respectively. The noise map units are in Jy 97 5.3 Mask used in source detection to fit the amplitude of a point centred on each pixel. Note the negative sidelobes. The central Gaussian has a full-width half-maximum of 14.7". Each pixel is 3" square. The mask is 150" square. The mask mean pixel intensity is zero 100 5.4 Subaru XMM Deep Field: Clockwise from upper-left, most-probable amplitude map, amplitude uncertainty map, and ratio map. The source detection acts as a statistically weighted low-pass filter. The funny pattern along the edges of both the most probable amplitude map and the ratio map is due to the averaging of only a few noisy values in the assignment of the pixel values near the edges. The saturated, mask sized regions in the upper right and lower left of the most probable amplitude map and the ratio map are due to very large valued pixels along the edge of the Subaru measurement map. . . . 101 List of Figures xvi 5.5 Subaru XMM Deep Field most probable amplitude pixel histograms with noise model. Top, histograms produced from pixels with noise less than.4mJy. Bottom, histograms made from pixels with noise less than 3mJy. There is an excess of positive, and negative flux in the uncleaned most probable amplitude map, due to the presence of sources. The apparent lack of false detections in the top panel is due to the greater pixel noise populating the top histogram's positive tail. 103 5.6 Subaru XMM Deep Field source detections. The black crosses mark the location of detected sources. The closest detections are around 70" from the map edge 104 5.7 Subaru XMM-Deep Field source list correspondence. There are 26 unpaired UBC detections, 25 unpaired detections from the SHADES community, and 35 matches 105 5.8 Sky estimate map pixel standard deviation as a function of iteration. The standard deviation is within 0.01 mjy of its converged value at iteration 100 106 5.9 401st sky estimate map. The pixel-pixel correlations are apparent. Measurements centred upon pixels involved in the centre of measure-ments less than thirty times are discarded. Source detection on this map will have to deal with the large scale structure, and if it is to be optimal, the cross-hatching 107 5.10 Deconvoluted source detection correspondence with the sources de-tected on the measurement map (1st iteration map) and with the source list provided by the rest of the SHADES community using the data acquired up to February 1st, 2004. There is no indication which source detection method is more reliable 109 List of Figures xvii 6.1 Simulation of data reduced using the central update algorithm. The data are simulated measurements of a statistically plausible fake sky with blank field extended noise samples. An extra chop's worth of data has been added to measurements made with the full SHADES measurement strategy. The chop has a size of 44" and is rotated 45 degrees with respect to the original chop directions. Clockwise from upper left, first iteration sky estimate map, 101st iteration sky estimate map, 101st iteration power spectrum map, first iteration power specrum map. The power along the dominant cross in the 101st power spectrum map in Figure 4.14 is somewhat reduced, though not as much as in Figure 6.1, and the cross-hatch pattern modified. . . . 115 6.2 Simulation of data reduced using the central update algorithm. The data are simulated measurements of a statistically plausible fake sky with blank field extended noise samples. Two extra chop's worth of data has been added to measurements made with the full SHADES measurement strategy. The chops have a size of 44" are rotated 45 degrees with respect to the original chop directions and are mutu-ally orthogonal. Clockwise from upper left, first iteration sky esti-mate map, 101st iteration sky estimate map, 101st iteration power spectrum map, first iteration power specrum map. Power along the prominant cross in Figures 6.1 and Figure 4.14 is further reduced. At low frequencies the power has been concentrated to narrow lines along the direction of the prominant cross 117 Acknowledgements xviii Acknowledgements I would like to acknowledge the people whose suggestions and insight were invalu-able: Mark Halpern, Colin Borys, Douglas Scott, and Guillaume Patanchon. Part I Thesis Chapter 1. The Scuba Half-Degree Extragalactic Survey 2 Chapter 1 The Scuba Half-Degree Extragalactic Survey The Scuba HAlf Degree Extragalactic Survey (SHADES) is a major extragalac-tic survey of 850/xm radiation obtained using the Submillimetre Common- User Bolometer Array (SCUBA) at the James Clerk Maxwell Telescope (JCMT). The survey is designed to cover half of a square degree to a 4cr detection limit of 8 mJy1 sources. By constraining the evolution, and the energy source of the population of galaxies detected at 850 /jm, the survey is designed to maximize the impact of JCMT on cosmology over the next three years. To do this, three specific questions are addressed. [7] 1. What is the cosmic history of massive dust-enshrouded star-formation activ-ity? 2. Are SCUBA sources the progenitors of present-day massive ellipticals? 3. What fraction of SCUBA sources harbour a dust-obscured active galactic nu-cleus (AGN)? In addition to these three main goals it is desirable to maximize the legacy value of the survey. ^ m J y = 1 0 - 2 6 W m - 2 H z _ 1 Chapter 1. The Scuba Half-Degree Extragalactic Survey 3 1.1 Introduction To The Submillimeter Galaxy Population Submillimetre radiation ranges from 200/xm to about one millimetre. It is advanta-geous to observe the early universe at these wavelengths because of an effect called positive k-correction. Positive k-correction is the compensation for inverse square distance attenuation of submillimetre wavelength detected radiation by the increase in rest frame luminosity of the emitted radiation detected in the submillimetre as a function of redshift. For observations at 850 /im, the k-correction largely com-pensates for the inverse square attenuation of electromagnetic radiation and the detected flux from an eight hundred fifty micron source is approximately constant between redshifts two and eight. Similarly, because of the negative k-correction at wavelengths other than the submillimetre band, high-redshift, submillimetre de-tected galaxies often lack counterpart detections. [4] High-redshift submillimetre bright sources are believed to be very luminous, dusty galaxies which account for a significant percentage of all energy produced by all galaxies during the age of the universe. [4] Understanding the clustering and evolution of the high-redshift population of submillimetre detected galaxies will constrain early universe structure formation models. To constrain the clustering of the submillimetre population requires a large survey of submillimetre flux and modest, unbiased redshift estimates of detected sources. [7] Obtaining redshifts for the high redshift submillimetre detected galaxy population is difficult, as the inferred submillimetre energy distribution of the popu-lation is a modified blackbody2 spectrum, lacking features and having few detectable spectral lines. The detected spectral lines consist mostly of CO emission, detections emanating predominantly from strongly gravitationally lensed galaxies. [4] If a galaxy detected at submillimetre wavelengths can be detected at other wavelengths, and the spectral emission distribution is known, it is possible to determine 2 "grey-body" Chapter 1. The Scuba Half-Degree Extragalactic Survey 4 where T is the temperature of the radiating grey-body, and z is its redshift.[4] Redshift estimates require temperature constraints. Methods of constraining the spectral energy distributions and the temperature of the submillimetre galactic population involve studying: 1. Theoretical models of galactic spectral energy distributions. 2. Galacticly lensed high redshift galaxies. 3. Low redshift submillimetre galaxies. 4. Low-redshift analogous galaxies, such as ULIRGs. Large, local galaxies with high star formation rates are bright in far infrared radia-tion, though stars themselves radiate negligibly in the far infrared. The far infrared radiation arises through absorption and re-radiation of the stellar radiation by heavy elements in the interstellar medium. The heavy elements are called "dust". Studies of the above constrain high redshift submillimetre galaxy spectral energy distribu-tions. The distribution of the temperatures of ULIRGs might constrain the expected temperature for high-redshift submillimetre galaxies, allowing for redshift estima-tion. Regardless, [10] [3] perform extensive Monte-Carlo simulations and claim the feasibility of determining the redshift of sources detected in the SHADES survey, to Sz = ±0.5, given follow-up observations from the Balloon Borne Large-Aperture Submillimeter Telescope (BLAST). BLAST will observe SHADES fields at three additional wavelengths, 250, 350, adn 500/mum. A submillimetre survey spanning a large fractin of a square degree is required to constrain the clustering of the submillimetre galaxy population. [7] 1.2 What Is The Cosmic History Of Massive Dust-Enshrouded Star Formation Activity? Star formation activity is determined by extrapolating observations of local star forming galaxies, detected in the infrared[15], to submillimetre observations. The Chapter 1. The Scuba Half-Degree Extragalactic Survey 5 suggested star-formation rate is SFR = 2.1xlO~iU x ) (1.2) LeOfim is determined by assuming a spectral energy distribution for the submillimetre detected sources and extrapolating to sixty micrometer luminosity. Note that small errors in the submillimetre source temperature can lead to large errors in Leofim because of the T 4 dependence of blackbody luminosity on temperature. Nonetheless, when combined with weak redshift constraints (6z = ±0.5), an estimate of the star formation rate as a function of redshift should be possible.[16] 1.3 Are SCUBA Sources The Progenitors Of Present-Day Massively Elliptical Galaxies? Currently, observations of high-redshift submillimetre galaxies are consistent with luminous, massive, dusty galaxies, similar to the local universe ULIRG, and also consistent with observations of moderate sized galaxies with optical and ultra-violet counterparts experiencing a period of rapid star formation. [22] One tenth of a per-cent of all local galaxies are ULIRGs, and if the former type of galaxy is the pre-dominant type of SCUBA source, then it is likely that they are the progenitors of today's massive elliptical galaxies. Massive elliptical galaxies trace the regions of over-density in the universe and experience strong clustering. Thus, if SCUBA sources are the progenitors of the massive ellipticals, they also, are expected to be highly clustered. 1.4 What Fraction Of SCUBA Sources Harbour A Dust-Obscured Active Galactic Nucleus? The spectral energy distribution (SED) of sources detected in the submillimetre is largely thermal and lacks characteristic features. The absence of sharp features in the SED makes it difficult to determine the source of radiation which is absorbed and thermally re-radiated by the dust. The possibilities are dust-obscured stars and Chapter 1. The Scuba Half-Degree Extragalactic Survey 6 the accretion disks surrounding dust-obscured black holes with masses greater than one billion solar masses. To detect active galatic nuclei (AGN), rest-frame soft x-ray-observations of submillimetre sources are required. 1.5 Legacy Value The legacy value of the SHADES survey is maximized by placing the data in a form maximally conducive to useful analyses. It is difficult to predict future use; the data need to be reduced in a careful, yet practical, fashion. Reduced data are compressed from the original data to an extent which strikes an acceptable compromise between information content and data quantity. The reduction process should introduce minimal, well understood spatial correlations to the data.3 1.5.1 D a t a Compress ion The raw data consist of many measurements of the SHADES region, and are too numerous to be easily used. Reduction is careful removal of the time component, and careful interpolation of the non-uniform data samples of the flux emanating from the SHADES field onto a uniformly spaced grid of points. Each point is associated with a small surrounding region of space and is assigned a value. The region of space, point, and value combination is a pixel. A pixelized image is a histogrammed approximation to a two dimensional function which is the flux radiated by the region observed as a function of spatial coordinates. 1.5.2 D a t a Cor re la t ion : Integrated Sachs-Wolfe Effect The accuracy of the reduced data image, or map, should be adequate and its pre-cision quantified. The degree to which spatial correlations present in the data are due to the data reduction algorithm should also be understood. If this is the case, then it is possible to detect physical effects which result in spatial correlations. An example of such a detection in maps of the Cosmic Microwave Background (CMB), is the detection of the integrated Sachs-Wolfe effect, an effect detectable in 3 This is a statement regarding the quantitative description of the data certainty. Chapter 1. The Scuba Half-Degree Extragalactic Survey 7 CMB images spanning much larger angular scales than one-half degree. The Sachs-Wolfe effect is the net gravitational doppler shift experienced by a photon as it traverses a region of high density. If the gravitational well expands at a constant rate during the photon's traverse then the doppler shift of the photon is the redshift the photon would experience during motion through flat, expanding space. If, however, the growth of the gravitational well changes during the traverse of the photon, then there is an anomalous doppler shift in addition to the redshift experienced by the photon due to constant universal expansion. By correlating the redshift's of photons radiating from the sphere of last scattering which pass through overdense regions of space, it is possible to significantly detect an acceleration to the universe.4 This detection would be difficult, if not impossible, if the images of the Cosmic Microwave Background (CMB) were significantly spatially correlated by the algorithm used to reduce the CMB data. [12] 4 Q A > 0 , 9 5 % C L , s t a t i s t i ca l er rors o n l y Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 8 Chapter 2 Observing The Scuba HAlf-Degree Extragalactic Survey The SHADES observation strategy, executed on the JCMT using SCUBA, is de-signed to provide sufficient and uniform coverage of a half-degree square patch of sky. The use of acquired data requires calibration prior to reduction. The subse-quent sections describe SCUBA, and the acquisition, pre-processing, and reduction of the data acquired during SHADES.1 2.1 Submillimeter Common User Bolometer Array SCUBA is a sparsely-sampled array of hexagonaly spaced, cryogenic bolometers, situated behind filters whose pass-band is centred on 450 and 850 micrometer wave-length radiation. Radiation at these wavelengths are relatively unattenuated by the atmosphere. The array is located at the JCMT focal plane. Figure 1 shows the spacing of the 850/xm SCUBA bolometers. 1 T y p i c a l l y i n the a s t r o n o m i c a l c o m m u n i t y d a t a p re -p rocess ing is r ega rded as pa r t o f d a t a re-d u c t i o n . In th i s thes is the d i s t i n c t i o n be tween d a t a r e d u c t i o n , as the ac t of r e d u c i n g the n u m b e r of pa rame te r s w h i c h descr ibe the d a t a , a n d p re -p rocess ing w h i c h is c a l i b r a t i o n , is m a d e for ease o f d iscourse . Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 9 850 Micron Bolometer Array -100 8 - 1 5 0 o - 2 0 0 h - 2 5 0 -300 - 1 5 0 -100 - 5 0 0 Right Ascension (arcseconds) Figure 2.1: SCUBA 850 micrometer bolometer array. The plot is obtained from the positions of actual data on the sky. Note the spacing between equal distances along the two axes are different, the height and width of the array are equal. 2.2 Data Acquisition The SHADES survey consists of 180 nights of observation in grade two weather, where grade one is best and grade six is the worst. Data are collected at 1 Hz by averaging the readings from two locations on the sky, which are observed alternately at a frequency of 7.825 Hz. The sky locations are varied by chopping, a process which refers to rapidly tilting the secondary mir-ror, or chopper. Every 16 seconds 2, the primary mirror is moved such that the position "chopped onto" is the position originally pointed to by the telescope. The data acquired in this new position are subtracted from the "single-differenced" data previously acquired, prior to the primary mirror movement. The result is "double-differenced" data, which removes any offsets, and signals of constant slope along the 2 T h e r e is a th ree second lag wh i l e the p r i m a r y m i r r o r is m o v e d . Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 10 direction of the chops, aiding in the removal of signal from atmospheric emission. Each sample from the double-differenced data is the normalized3 subtraction of two single-difference samples acquired 19 seconds apart. The focal plane of the JCMT is undersampled by SCUBA both at the 450 pm and 850 fxm wavelengths. Thus, a small, 1 Hz motion of the secondary mirror, called "jiggling" is used to acquire data located between the bolometer detectors in the focal plane, and "fully sample the array"4. In particular, the jiggle pattern is a 16 point pattern repeated every nineteen seconds designed to Nyquist sample the 850 detectors. A 64 point pattern, which repeats every four primary mirror movements, is used to nyquist sample the 450 //m detectors. The difference in the number of points in the jiggle patterns for the two wavelengths is due to the nar-rowing of the telescope point spread function from the 850 /zm wavelength to the 450 fim wavelength, ie. The spatial Nyquist frequency is increased for detection of 450 /jm radiation. The jiggle pattern is specified in elevation and azimuth, and also rotates with respect to the celestial sphere. Observation of each SHADES pointing involves a 64 point jiggle pattern with three chop "throws" specified in two directions on the celestial sphere.5 The chop throws are 30", 44", and 68" in the east-west direction and 30", 44", and 68"6 in the north-south direction. The chop throws are specified in declination and right ascension. These chop throws are chosen such that each chop throw measures a different spatial frequency harmonic and links many different pixels in different measurements. More on this topic will be discussed in section 6.3 on measurement strategy. The SHADES field is covered by many pointings of the type described above. A relatively uniform observation coverage of the field is obtained through "tripos" pointings, a scheme whereby three pointings, whose pointing centres are offset to form three overlapping circles of sky coverage, or a tripos, are panned over the SHADES region. Figure 2.2 is a cartoon sketch of the coverage of a single tripos 3The double-difference is divided by two. 4Nyquist sample. 5Right Ascension/Declination (J2000 coordinates). 61" is one arcsecond. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 11 pointing. [17] Each circle in the tripos pointing is actually hexagonal and spans the 2 Figure 2.2: Tripos Coverage array size, roughly 150". As previously described, each pointing in the tripos consists of the six chop throws in the two directions, centred in the centre of one of the circles in Figure 2.2. 2.3 S C U B A / J C M T Measurement Noise The measurements acquired by the SCUBA detector are noisy. The precision of the reduced data depends upon the noise present in the unreduced measurements as well as how the reduction process affects the measurement noise. Much of the design of JCMT and SCUBA is motivated by a desire to mitigate the effects of noise. 2.3.1 Brief Overview Of Noise The measurement of a physical quantity is composed of the quantity of interest plus one or more contaminating or nuisance signals associated with the measurement process. In principle, the nuisance signals can be computed based upon classical physics and removed.7 In practice, a physical description of the nuisance signals 7 T h i s d i scuss ion assumes a l l measu remen ts c a n be desc r i bed b y c lass ica l phys i cs . Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 12 is very difficult and probablilistic descriptions of these signals are employed. A nuisance signal described by probabilistic properties is called noise. The spectrum of typical measurement noise is a combination of two processes8 which are described by a particular function of frequency. The two processes are " j " noise and "white" noise, j exhibits a monotonic decrease with increasing frequency and white noise is constant over all frequencies, j noise dominates at low frequencies and white noise dominates at frequencies greater than the " j knee", or the frequency at which the white and j noise components to the noise spectrum are equal, j noise indicates power at low frequencies and corresponds to correlations in the noise. A typical assumption made of the two processes comprising the time stream spectrum are that they are Gaussian distributed. That is, the random variables comprising the process are described by Gaussian probability density functions. This approximation, when describing noise, is supported by the Central Limit Theorem. The Central Limit Theorem states that probability density functions of the sum of many random variables, regardless of their density functions, tends to a Gaussian probability density function as the number of random variables in the sum tends to infinity. [14] This claim is consistent with physical effects causing the noise; that is, the random nature of the noise is due to complexity. Alternative justification for the use of Gaussian probability density functions comes from a maximum entropy arguement, whose result is that a Gaussian probability density function assumes the least amount of knowledge about a process under certain, commonly encountered, conditions. [20] Samples from a random process9 can be dependent. Probability density func-tions for the random variables comprising a process with dependent samples are multivariate10 and non-separable. Data sampled from such a random process results in correlated samples. When detecting signals in a noisy data sequence containing correlations, it is necessary to account for the correlations. An optimal method of fitting to data with noise sampled from a Gaussian process is given in section 2.5 on 8 A random process is a collection of random variables[14]. When discussing a sequence of data points, the random process is the entire sequence and each data point is a random variable. 9Collection of random variables. 1 0 A multivariate probability density function is called a joint-probability density function. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 13 data reduction. Gaussian probability density functions of a single variable are com-pletely characterized by their mean value and their variance. With a multivariate Gaussian random variable, the process is completely described by its mean and its covariance matrix. The diagonal entries of the covariance matrix give the variances of the variables, and the off-diagonal terms give a measure of the amount with which each variable varies with the others. For a sequence of data values sampled from a random process, given by the vector n, the covariance matrix is N = (nnT) (2.1) and the correlation matrix is N~l. Note that N is symmetric and has dimension equal to the length of the vecotr n. Here, each data sample is considered a random variable in the Gaussian random process which is sampled to produce the vector of data samples, ft. When attempting to detect weak signals or large signals to high precision, one performs a low signal to noise ratio experiment. In such experiments, effort is made to reduce the noise correlations such that averaging TV measurements results in a reduction of the standard deviation of the averaged measurements by a factor of VN. 1 1 If all noise correlations have been removed from the signal, an arbitrar-ily large signal to noise ratio can be obtained given a sufficiently large number of measurements. A noise signal possessing correlations, such as a slowly varying drift, cannot be averaged to achieve an arbitrarily high signal to noise ratio because the average of the slow drift will eventually dominate the error. The drift appears as j noise in the power spectrum of the signal. 2.3.2 SCUBA Noise The dominant noise present at the SCUBA detector comes from atmospheric emis-sion. The remaining sources of noise are comprised of confusion noise, photon detec-tion statistical noise, Johnson noise, noise due to non-ideal electronic components n T h e variances of Gaussian random variables add as their inverse square. Thus, the variance of the sum of N measurements is a factor of y/N larger if all of the measurements are sampled from the same probability density function. Division by N to obtain the average value reduces the standard deviation of the average of the measurements by a factor i V ~ 0 ' 5 . Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 14 (ie. physical components), quantization noise, etc. Combined, the effect of the noise is to give j noise and Gaussian white noise components to the signal. Atmospheric emission, to good approximation, is constant over the patch of sky subtended by the SCUBA bolometer array.[1] Astronomical observation of the SHADES region is a low signal to noise ratio experiment. The instantaneous signal to noise ratio is approximately 3% for a 10 mjy source, and considerable averaging is required to detect SHADES sources at the desired 4cr level, where u is the standard deviation of the beam convolved pixels, and is desired to be 2.5mJy.[17] 2.4 D a t a P r e p r o c e s s i n g Prior to data reduction the data must be calibrated, outliers rejected, and nuisance signals understood, estimated, and if possible, removed. 2.4.1 Calibration SCUBA measurements are calibrated to account for unequal bolometer gains, at-mospheric absorption and to acquire the correct conversion from detector signal in Volts to detected flux in milliJanskies. 1 2 The SHADES survey is granted time in grade two and upper grade three weather. Weather, in this context, refers to the current opacity of the atmosphere to 850 /um radiation. The weather is quoted in terms of a fitted parameter in a model of atmospheric absorption. The model is / = I0e-TA (2.2) where IQ is the intensity of 850/im radiation incident on the atmosphere, / the intensity of 850 //m radiation incident on the focal plane, r the 850 /an zenith opacity, and A the airmass.13 Grade two to upper grade three weather corresponds to a zenith tau measured by the CSO14 at 225 GHz between 0.05 and 0.10.[17] The 12lJansky = 1 0 - 2 6 W m ^ H z - 1 1 3 Airmass is inversely proportional to elevation and incorporates the increased amount of atmo-sphere along the line of site into equation 2.2. 1 4Caltech Submillimeter Observatory Chapter 2. Observing The. Scuba HAlf-Degree Extragalactic Survey 15 225 GHz zenith tau is related to the 850 /im tau through relations available on the JCMT webpage. The relations are determined by correlating the CSO and have changed with time. The data are adjusted by the JCMT standard SURF procedures, flatfield and extinction, correcting for unequal bolometer gains and atmospheric absorption, re-spectively. The flux conversion factor, or FCF, in units of Janskies per Volt, is computed by reducing calibrators, fitting a two dimensional Gaussian to the result-ing map of the calibrator, and comparing the integrated flux within the full-width half maximum of the fitted Gaussian to a known value. For a given night of obser-vation there are several calibration observations, and the FCF computed from these observations is not constant. Figure 2.4 is a plot of the FCF multiplying the data as a function of data observation file. These FCFs are all computed using CRL618 using Method 3, described below. There is no discrimination between FCFs com-puted using calibration observations taken with different sized chop throws. The observation files contain measurements of the Subaru XMM-Deep Field from De-cember 18, 2002 to June 22, 2003. There are 292 observation files each containing 320 double-difference measurements for each of the 37 bolometers in the SCUBA array. The observation files are not in time order. Deciding which FCF to apply to a given observation file of data is complicated. The following methods were tried: 1. Set the FCF to a nominal value of 225 Jy V - 1. 2. Average the first and last computed FCF for a night of observation. Use a running average of the FCFs if the FCF is above 400 Jy V - 1 or below lOOJyV-1.16 3. Average the FCFs computed from calibration observations performed most recently with respect to the time of the observation to be calibrated, and 1 5 A running average is an average of every FCF that is actually used on the data during the data reduction up to the current point in the program execution. The running average is re-computed with every computed FCF that will be used in the conversion of the data from V to Jy. 1 6 FCFs above 400 Jy V - 1 and below 100 Jy V - 1 are deemed bad and neglected from the com-putation. 400 and 100 are chosen as many standard deviations from the average. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 16 Dote? ^eSi~^!Ee' Qnl£>rn$-s' 1 ' 3.3*11? j-o-ia5, LSx^a* Bctio Sequence; Bolometer 1ft J a.'B«iD* i.a*riffi^' i.sma* Dose Sequence: Bolometer 20 3.a.fio* i.Qwin* i.a«io> s.c^ cfl Doioi Sequence; Bnlorn&tef 21 a.c*io* i.c^ic^ i.a-oa* U H Qai« Sequence; Snipmeler 2.S 3.3*!ry ijQ*'H»* i.iiaa* i'.e^ff Dalo Sequence: Solomcicr 24 I . O - I B * !.i»'ia* 1 ,0**0* a.o*!D* I J O - W I S * s.ax'.a1 Sg^clryrn Mgqnij^ rJB^  g^ornftlf;r 1 9 o.-oooi C - J M I Q a e n c o o.toan Spectrum Mcqniiudc: BoJorntjier 19 O.CDtn DiODIO 0.01430 Spectrum Moqniiiide: Bdorr'Cis r 20 1 j | p vmmt a.floai arena a.oiflo c.iooc U D D Q Sfj*u;tfum iySaqnity«»i-; EM^ metsr 2'f D-OS10 0.0 5 CO 0.! COEJ Sp^ryrr? MpqnlUi^p: gi^omsfor ?.Z aomo a.rjico o.irxiC' Sp^-run-i v.rjcniiuiic: 3<srigrrBrtcr £3 C.DI'BO 0.5030 i»-*E , . . ' to-4. . . . . D.ocia o.o'.co Q . I D O S &OOIQ 0.0130 G.IGOD ^3CO0 Figure 2.3: Eight representative scuba bolometer timestreams and spectra for the Subaru XMM Deep Field after flatfielding, and extinction cor-rection. The timestreams are non-stationary and strongly correlated by atmospheric emission. The timestream for bolometer 21 extends past the range of the y-axis. The standard deviation of the bolome-ter timestreams varies substantially from bolometer to bolometer. The magnitude of the Fourier transform of the bolometer timestreams are shown in the lower half of the plot. Note the flat spectra and the spikes on bolometers 17, 19, 20, 21, and possibly 23 and 24. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 17 Flux Conversion Factor vs Observation Number 300 250 200 150 100 200 Observation Number 300 Figure 2.4: Flux conversion factor as a function of observation file number. The observation files are not in time order. Note the considerable fluctuation of the FCF with observation number. Only calibrations with CRL618, Uranus and Mars are used. The chop throw of the calibrations is ignored. the soonest to be performed after the observation to be calibrated has been completed. If there are only FCFs before or after the observation, use the FCF computed from a calibration observation that occurred closest in time to the observation to be calibrated. Use only FCFs which were computed after the most recent focus. Do not use any FCFs above 400 J y V - 1 or below 100 JyV - 1. If under these conditions there are no FCFs, use the running average. The locations and number of 3cr sources detected in a measurement17 map of the Subaru XMM -Deep Field demonstrated a slight dependence upon the method of FCF computation. To quantify the effect, the sources detected in two maps of the Subaru XMM Deep Field are plotted in Figure 2.5. One of the maps is produced 1 7 A measurement map is a non-deconvolved map obtained by grouping all of the measurements centred upon a pixel, averaging them, and assigning the average to the pixel value. This is done for all pixels of the map. Measurement maps are described more fully in section 2.5.1 on Interpolative Reduction. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 18 using measurements calibrated with FCFs computed using the third method above, the other with measurements calibrated using a nominal FCF value of 225 JyV - 1. The source detection algorithm is described in section 5.2. The correspondence be-tween sources detected from data reduced with Method 1 calibration and sources detected from data reduced with Method 3 calibration is 90%. Here correspon-dence referes only to location and not to the significance of the detection. For the remainder of this thesis, data reduction is performed with an FCF set to 225 Jy V - 1. Despite the correspondence demonstrated in 2.5, the usefulness of the calibra-tions depends on the degree of fluctuation of the FCF values, as well as on the intended analyses. For instance, the effect of calibration error on clustering compu-tations is unknown and not constrained by this test. If it were desirable to estimate the flux to greater precision, and the statistical error were not dominant, then ac-counting for the fluctuations in the FCF would be necessary. To estimate the effect of FCF fluctuation on the analysis intended, realistic simulations must be performed. Because of the quantity of calibration data taken during the SHADES observations, the statistical properties of the FCF can be determined and used in subsequent simulation.18 2.4.2 Nuisance Signals And Outlier Rejection For SCUBA data taken by double differencing, the dominant nuisance signal is atmo-spheric emission. Atmospheric emission is approximately constant for all bolometers on any bolometer array sample. [1] Atmospheric emission and outliers present in the data are removed iteratively by removing 5cr outliers from the timestreams, then, for every array sample, subtracting the average of the bolometer samples, and re-peating until convergence. [19] This process is further complicated by the presence of spikes at 1/16 Hz in the magnitude of some of the bolometer timestream Fourier transforms. 1 9 When subtracting the average of the array, if there is a coherent component to 1 8 This idea was communicated by Tim Jenness to the SHADES collaboration through email. 1 9 In this context, units of Hertz are not completely accurate, due to the delays in data acquisition between telescope nods. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 19 ~i 1 1 r - i — i — i — i i r - | i i i r -+ Calibrated O Uncalibrated 2h o c 73 Q O - 2 - 1 0 _i i i_ -5 0 5 Right Ascension (arcmin) Figure 2.5: Correspondence between sources detected at greater than 3a signifi-cance from calibrated data and from uncalibrated data. The calibrated data uses Method 3 in Section 2.4.1 to compute flux conversion factors. The source detection is described in Section 5.2, but differs in that the detected sources have not been iteratively cleaned from the map. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 20 the spike, then the value subtracted from the bolometer timestreams will contain the coherent average of the spiky bolometer values. To avoid this, when subtracting the average of the array at every time step from the individual bolometer readings, the array is divided into two parts. One part of the array is composed of bolometers which do not exhibit the spectrum spike, and the other part of the array is composed of bolometers which do. The spike is detected by banding the magnitude of the Fourier transform of each bolometer timestream into three equally sized bands in frequency. There are no gaps between the bands of frequencies, and the central band is centred upon ^  Hz. Within each band, the sum of the Fourier transform magnitude is performed over the band frequencies. The sums for the first frequency band and the last frequency band are averaged and compared to the sum for the central frequency band. If the central band sum is greater by a threshold, then the spike is considered detected. The average of the bolometers in the different axray parts is computed and subtracted from the bolometers which belong to the part of the array from which the subtracted average is computed. When computing the array averages, bolometers exhibiting a standard deviation three standard deviations away from the mean of the bolometer standard deviations for the observation are not included in the computation. I will refer to this computation as array average subtraction for the remainder of this thesis. The frequencies at which astrophysical sources appear in the data are at the harmonics of the 64 point jiggle pattern. That is, at l/64Hz, 1/32 Hz, 1/16 Hz, etc. The harmonics are formed from the repetition of the jiggle pattern every 64 seconds. Note that in the absence of the jiggle pattern, when pointed at a source, the chopping and nodding double-difference measurement strategy will produce a constant signal, which appears at zero frequency in the double-differenced data. Additional Fourier space structure will be caused as the bolometers move on and off source due to the jiggle pattern. This is expected to broaden the spectrum centred on the 1/64 Hz harmonics. Thus, bandstop filtering, with the stopband centred upon the 1/16 Hz spike is not an option. Lastly, the array average subtraction is performed iteratively along with 5a out-lier rejection. That is, the array average subtraction is performed, then measure-Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 21 merits five standard deviations from the mean of the measurements for any given ob-servation file (320 double-difference measurements per observation file) are removed. These two computations are repeated four times. Figure 2.5 is a plot of the clean timestreams and the magnitude of the Fourier transform for a representative sample of the bolometers. The preprocessed timestreams have a smaller standard deviation than their non-preprocessed counterparts and have no outliers. The magnitude of the Fourier transform of the preprocessed bolometer timestreams is flatter than the non-preprocessed bolometer timestream Fourier transform magnitudes. (Figure 2.3) 2.5 Data Reduction The final step in obtaining useable data is data reduction. The desired properties of the reduced data are stated in section 1.4. This section describes common methods of data reduction, their strengths and weaknesses, and relevant issues. It is beneficial when discussing the methods of data reduction to describe the process using linear algebra. To this end, beginning with the single-difference mea-surement, the telescope measurements are described as equations relating measure-ments to linear combinations of the observed sky pixels. When doing this, a useful approximation is to approximate the telescope point-spread function as a Dirac-delta function.20 A double-difference measurement becomes a relationship between three pixels, and a single-difference measurement is given by: Mlc = -Po+pi (2.3) where po is the pixel in the chop location of the single difference, pi is the pixel at which the telescope points, and M\c is the left-chop single-difference measurement. A second single-difference measurement is formed after the primary mirror is moved: MTC = -Pl+p2 (2.4) Now the chop is looking at the position first pointed to by the telescope, and the telescope points at a new pixel, p2- The double-difference measurement is obtained 2 0 This approximation can be relaxed by maintaining the Dirac-delta telescope point-spread func-tion, and convolving the true sky map with the actual telescope point-spread function. This is done in simulations later on. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 22 Dfli« Sequence: BnlnmMer M i j 0 a.a-fcin' J , Q * I B * i.a*io* Sflio Haoycnce; EJoipmeter *9 1 j ;.L _ _ J Quia Sei^ yssnce: Botorne!;*?1 a.fl.«iD* u * ! c s i . s * i a » z.a*'*(P Daio Sequence: Boiofnclef ?1 S .a-lD^ '.0**2* ^ft-ilO* fi-DS'ioa Dai-o Sentience: ti^ 'or'"*?^^ 23 :iL i.a».<c3 L S M O * ?.v*iiP to-* ! -Spcctfurrt Msijoriuar:; B r i o ^ ' " i P L J 1 !}-*-£. o.scei B .OO ! c a.sicra &. was S^ctrM^^ ? l*PgrtJfff? c : Hctlofneicf 10 O.CK!3 CUB 1C 0.0 '< OB 3 a CO Stwc^rn Mcgniiijdi?; Bot&mdsr #1 O.OCQ! D.OJM 0 fl.'B f CM C, 10GO i . 3 * i D J ! .Q " IC5 I.&K'.CJ* Z.&xifX-[Join Secuflncc: Bulnrncic _ , ,—. , ,— , S.C-'D' ; jO»lC f t i.5*?C* 2.rj«tC Dale Sflquflnc*;; Bolometer 22 <##» 3.3*10* LO*tC* L5«1t f i Da In 5eny«nco; Bolfifncis?1 24 1 5 .0* 1 C {-C« 1 C * J . 3 * 5 0 * '2 Sir-1 CP _ SpcciHjm IVo j^iiJijde: B^yrrisigi JB fT^II!li!ll |p | ( | | | j 10" OJKJt C.OaiC Q.CIOD 5.WCO o.aooi oxflic o.aiao ciaco 1B~ arafo cunon M O C C Spflctru^t MrjgfiHugrr: Bpiorrysiftr 7J-f>. 00.10 0.2103 Figure 2.6: Eight bolometer timestreams and their corresponding Fourier transform magnitudes for the Subaru XMM Deep Field after preprocessing. Note the lack of outliers in the preprocessed timestreams as well as the reduced rms when compared to the timestreams in Figure 2.3. The Fourier transform magnitudes are flatter than their non-preprocessed counter-parts. (Figure 2.3) Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 23 by subtracting the two single- difference measurements and dividing by 2, obtaining: M i c - M r c -po + 2pi - p2 . . Mdd = g = 2 ^ ' Note that in actuality, po is the telescope point-spread function weighted sum of sky flux centred upon the centre of pixel po. The above approximation, using the Dirac-delta telescope point-spread function, is one which introduces a small amount of error into reduced data maps, and is used in all data reduction documented in this thesis. [19] Lastly, realistic measurements contain noise: where n is the noise. Some methods of reducing observed data are: 1. Interpolative Reduction 2. Emerson Reconstruction^] 3. Direct Inversion 4. Iterative Inversion Methods one through four will be discussed in the following sections. 2.5.1 Interpolat ive Reduc t ion Interpolative reduction is performed by fitting a function to measurements of the sky within some region of the pixel whose value is to be determined. The value of the fitted function at the pixel centre is assigned to the pixel. Fitting and pixel assignment are repeated for every pixel in the map. Commonly used fitting functions are two-dimensional Gaussians, planes, and step functions. Note that the act of fitting to a step function which is pixel shaped and non-zero only inside the pixel of interest is the same as averaging all measurements of the sky located within the pixel of interest. To see this consider x2: (2.6) N (2.7) i=i Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 24 where y is the pixel value, or the height of the step function, di is the ith data sample located within the pixel, and N is the number of data samples within the pixel. Taking the derivative of x2 and setting it to zero, obtain: N N = (2-8) i=l i=l and v = (2-9) which gives the pixel value as the average of the N measurements lying within the pixel. Maps created using the above technique will be called "measurement maps" for the duration of the thesis. Figure 2.6 is a measurement map created by reducing simulated measurements of a source. The simulated measurements are taken of a fake sky consisting of a single source, using the SHADES observation strategy for a single SHADES pointing. There is an equivalent and informative, linear algebraic formulation of interpola-tive data reduction. The linear system, Px = M (2.10) where P is an Nm measurements by Np pointings matrix, x is an iVp pixel vector of pixel values and M is an Nm vector of measurements. Equation 2.7 describes noise free measurement of the sky estimate, x, to get the data, M. In the case of data reduction, M is known and the task is to determine x. This is the inverse problem. Since the linear system is overdetermined and the measurements noisy, an optimal map, x is desired. The inverse solution, optimal in a generalized least-squares sense, is x = (PTN-1P)~1PTN-1M (2.11) Where N is the noise covariance matrix defined in equation 2.1. The element of the Nm by Nm noise covariance matrix for the ith row and jth column is given by the covariance of nj with rij where nj is the noise added to the ith measurement.[9] With the correct choice of P, the least-squares solution, x, given by equation 2.6, is equivalent to determining the interpolated pixel values. Consider the case Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 25 Figure 2.7: A measurement map of simulated measurements taken with the single pointing SHADES observation strategy of a fake sky with a single, cen-tred source. Note the presence of negative sidelobes centred upon the six chop throws surrounding the source. Each square pixel in this map is three arcseconds square. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 26 where the interpolating function is a two-dimensional step function, zero everywhere, except for the region of the pixel, where its value is one. Let the noise be uncorrelated and have unit variance. Then N is equal to the identity matrix. Further, let P be composed of rows that are all zero except for a single column in the each row, which is one. Now, consider the jth row of Px. The location of the one in the jth r o w of P is such that upon performing the multiplication, Px, the one in the row of columns multiplies the pixel upon which the measurement comprising the jth element of M is centred. More concisely, (PS)j mod Np = X3 mod Np = Mj (2.12) Thus, PTx forms a vector of length Np whose ith entry is the sum of all measurements centred upon the ith pixel, and (PTP)~x is an Np by Np diagonal matrix whose entries are -»r—, where -/V0bs is the number of times each pixel is centred upon in a measurement. Thus, the assigned pixel values are averages of the measurements centred upon each pixel. If we generalize to uncorrelated noise with some variance, the noise covariance matrix becomes diagonal but not the identity, and the resulting pixels are the inverse-variance weighted sum of the measurements which were taken when the telescope was pointing at each pixel. It is now desirable to determine the uncertainty of the pixel values when the data are noisy. Let M = Ps + n (2.13) where s is the flux emanating from the sky sampled at the pixel centres, M is the pointing matrix of the above section, and n is the vector of noise samples of length Nm. Every nf1 noise sample comes from a probability density, which is in general multivariate, and non-Gaussian . Assuming that the noise is sampled from a Gaussian process with zero mean, then the noise probability distribution is completely specified by the noise covariance matrix, as previously mentioned in section 2.3.1. If N is diagonal (uncorrelated noise), and the noise distributed with a variance, a2, then the pixel-pixel noise covariance matrix is given by E = a2(PTP)~1 (2.14) Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 27 In general, the pixel-pixel covariance matrix is given by (xxT) = ssT + £ (2.15) where the pixel-pixel noise covariance matrix, E, £ = (MTN~lM)~l (2.16) Here £ gives the uncertainty of the pixel parameters. 2 1 2.5.2 Emerson Reconst ruc t ion Emerson reconstruction^] is a method of removing the negative side-lobes around sources from a measurement map formed from differential measurements. The method uses Fourier techniques and requires the size and direction of the chop throws to be constant for all observations. The measurement map, x, can be con-sidered as a two-dimensional convolution of a measurement function with the sky flux: x = c®s + n (2.17) where c is the convolution function and n is the noise. By inverse Fourier trans-forming the quotient of x and c one obtains the deconvolved map. p = s + T~1(^j (2.18) where n is the temporal Fourier transform of the noise and c is the Fourier transform of the convolution function. In the pencil point spread function approximation for a one-dimensional map, c = -0.55 (x-x0- Ax) + 6(x-x0) - 0.56 (x - x0 + Ax) (2.19) where x0 is the position pointed at by the telescope and Ax is the chop size. Thus, c = exp ikx0 [exp ikx — 1 cos kAx\ (2.20) 2 1 See, for e x a m p l e , [20] for t he c o n n e c t i o n be tween l inear regress ion a n d the m a x i m a l l y l i ke ly so lu t i on . Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 28 and zeros in the magnitude of c cause spikes in j, introducing sinusoids into the deconvolved noise, correlating the pixels. This method of obtaining deconvolved reduced-data is not considered further in this thesis. For a full critique of this method see [11] [5] 2.5.3 D i rec t Inversion Direct inversion is solving for x given M in equation 2.10 by brute force, where P, x, and M are given in the interpolative reduction example, but with a pointing matrix P, given by the actual measurement strategy. That is, the ,7th row of P has a —0.5 in each column multiplying the chop pixels in the j t h measurement, and a 1 in each column multiplying the pixel centred upon in the jth measurement. There are two problems with this approach: 1. PTP is singular. 2. PTP is huge. The unmeasured modes of the sky manifest themselves as singularities in PTP. The inversion of PTP can be computed, up to a linear combination of the singular eigenvectors of PTP by computing the pseudo-inverse of PTP. [21] A pseudo-inverse is an inverse of the space spanned by the matrix, and yields the inverse of P up to a linear combination of the eigenvectors spanning the null space. The second problem is more serious. PTP is Np rows by Np columns and contains iV 2 elements. For three arcsecond pixels, SHADES consists of 360000 pixels.22 Thus, the computation involves the inversion of a 1.3 x 1011 element matrix. If each element is represented by eight bytes, the matrix is stored in 1000 Gigabytes of memory. This is 10 100 GByte hard drives. Standard matrix inversion techniques require order N3 operations to perform, in our case, this is roughly 1033 floating point operations. On a GigaHertz machine, assuming 109 floating point operations per second, one completes the inversion in approximately 1024 seconds, or in 5.5 • 1016 months.23 Exploitation of symmetries in the inversion may reduce 2 2 The survey is 1800 axcseconds by 1800 axcseconds in extent. 2 3Millions of Hubble times. Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 29 the number of operations required by a factor, but will not reduce the O (N3) dependence.[13] 2 4 Supercomputers are an option, but are not well suited to serial computations such as matrix inversion. They rely upon performing many parallel and independent computations simulateneously. The time required to perform the brute-force inversion is prohibitive, and due to the existence of a faster, iterative approach, it is not pursued. 2.5.4 Iterative Inversion In the mid-1990's Cosmic Microwave Background (CMB) experiments became large surveys of the sky. The number of pixel parameters to be determined in these surveys caused experimentalists to look for another way to solve equation 2.8. Ned Wright developed an iterative algorithm for determining the best fit pixel values which requires 0(Nm) operations, and the idea was adopted and extended by others.[23] The algorithm uses the data in timestream order, eliminating the need to store large matrices. For a single-difference experiment, the type of experiment in which the algorithm has been used, the iterative update scheme is given by where xf is the iih pixel of map x at the n t h iteration, Nm is the number of measure-ments, j is the measurement index, m,j is the jth measurement, <r2 is the variance of the j t h measurement, and A is the chop distance in pixels. [5] In words, when updating a pixel, assume the rest of the pixels comprising the map are correct. For each measurement involving the pixel to be updated, compute what the pixel must be if the other pixel in the measurement, and the measurement itself are correct. Thus, there will be a "correct" pixel update value for each measurement involving the pixel to be updated. Compute a weighted average of the "correct" pixel update values and assign the weighted average to the updated pixel. The weights are the variance of the jth measurement. This is Wright's algorithm generalized to noisy 2 4 Not completely true; circulant matrices are at least one exception. ,j=l (2.21) Chapter 2. Observing The Scuba HAlf-Degree Extragalactic Survey 30 measurements.25 A hallmark of this algorithm is that the next iteration pixel value does not depend on the current pixel value. Note that on any one iteration there are conflicting assumptions. The update of a pixel requires that all of the other pixels are correct. Thus, when updating another pixel, the pixel just updated is assumed correct, and one of the pixels previously assumed correct, is now assumed incorrect. It is thought that this is why iteration is necessary, and that iteration will converge to the correct answer, a map optimal in the generalized least-squares sense.26 Prior to this work, the algorithm described by equation 2.21 has been success-fully employed in reducing the data from many CMB experiments. Most notably, the algorithm has been successfully used in the reduction of Wilkinson Microwave Anisotropy Probe (WMAP) data. [9] The success of the algorithm has inspired its use in the reduction of data acquired from the SHADES field using SCUBA. The application of this algorithm to SHADES data comprises the main topic of this thesis. 2 5 The algorithm in this form does not take into account correlations between the data samples. 2 6Provided the timestream noise covariance matrix is diagonal. Chapter 3. Iterative Reduction of SHADES Data 31 Chapter 3 Iterative Reduction of SHADES Data In this chapter the iterative reduction algorithm, equation 2.21, is generalized for use with an arbitrary measurement strategy. The generalized algorithm is studied theoretically and by simulation. 3.1 Wright's Algorithm Generalized Equation 2.21 is valid for single-difference data and must be modified for double-difference data prior to its use in SHADES data reduction. The generalized al-gorithm retains the essence of the original algorithm. For each pixel involved in any one measurement, assume all of the remaining pixels in the measurement are correct. Compute the value of the pixel under consideration which will satisfy the equation between the remaining pixels , which are assumed correct, and the value of the measurement. Repeat this for all pixels involved in all measurements. For every pixel in the updated map, the above procedure will yield many estimates for each pixel. Averaging the estimates for each pixel gives the updated map. This is the same procedure used to determine equation 2.21, except now there are more pixels involved in a single measurement. The generalized version of equation 2.21 is: Equation 3.1 gives the n t h update of the i pixel in map x, where Nm is the number of measurements, j is the j t h element of the measurement vector M, and is the gain or factor multiplying the ith pixel in the computation (acquisition) of the (3.1) Chapter 3. Iterative Reduction of SHADES Data 32 (AP) t h measurement. The denominator is a sum of all gains, in all measurements, multiplying the i t h pixel. In this form, the iterative deconvolution update can be used to reduce data taken with an arbitrary point-spread function, and an arbitrary measurement strategy. For the reduction of SHADES data, assuming Dirac-delta telescope point-spread function, the gains are - 0 . 5 for gains multiplying the pixels at the chopped positions and 1 for gains multiplying the pixel upon which the measurement is centred. This is the double-difference measurement described in Section 2.6: M = -0.5pieftchop +Pcentre ~ 0.5p ri g nt chop (3-2) 3.1.1 A l g o r i t h m Checks There are two easy checks following from the considerations cited below. 1. The update of any pixel does not depend upon the previous value of the pixel being updated. This is a direct result of the assumptions made in the formation of equation 2.21 2. The limit of the pixel update should be zero, as the gain tends to zero for any single measurement. That is, when the pixel being updated is not involved in a measurement, the algorithm should be well-behaved. Check #1 Consider the update given by a single measurement, ie. Nm = 1, and consider the update for the ith pixel. The update is given by s? = ^ 2 - " ' +x?- X (3-3) \9i) The gains in the numerator and denominator cancel when multiplying the x™~1 term in the numerator of the update, with the result that the trailing x™~1 term is cancelled by the in the first term of the update. Since this occurrs for every measurement, the update of any pixel does not depend on the previous estimate of itself. Chapter 3. Iterative Reduction of SHADES Data 33 Check #2 Consider the contribution to the update of the x\h pixel from a single M-7 in a dataset with Nm measurements, where the jth measurement does not involve the pixel being updated. That is, consider the limit of equation 3.1 as gj tends to zero when g\ / 0, for all k ^ j: » m x n = & ^ A m * r + x r > ( 3 . 4 ) 9 i - ° (si) Breaking the first term into two fractions, obtain, B m ^ — ( 3 . 5 ) ^° Ef5 (si) Ef", (si)' When k = i in the second term, the second term will cancel the right-most term. Explicitly performing the cancellation yields, lim *? = - A T , T ^ ~ 2 (3-6) which equals zero. 3.2 Generalized Algorithm Simulations To characterize the generalized iterative algorithm, simulations of data reduction using equation 3.1 are performed on data acquired during mock experiments. Each mock experiment yields measurements of a fake sky taken according to a measure-ment strategy. 1 Fake skies measured vary from simple, seven bin, one dimensional discrete skies with a single Kronecker-delta source, to two dimensional discrete skies populated with Gaussian convolved point sources sampled from a realistic source count model.2 Like the fake skies, simulated measurement strategies range from simple, single-chop, single-difference measurement strategies, to the relatively com-plicated SHADES measurement strategy. 1 A fake sky is a map of the sky based upon some model of astrophysical flux. 2See [6] Chapter 3. Iterative Reduction of SHADES Data 34 In each simulation an initial map estimate is made, and equation 3.1 iteratively applied. Convergence is practically determined to be when the change in the sum of square pixel values is below a particular threshold. 3.2.1 SHADES Simulation The full SHADES measurement strategy is applied to a half-degree square fake sky populated with JCMT point-spread function convolved point sources sampled from a realistic submillimetre source count model.3 The original map estimate is taken to be the null map.4 Figure 3.1 shows the fake sky map, and the residual maps for various iterations. The colour bax at the bottom of the Figure applies to all of the maps except the fake sky map, which saturates at pixel values of zero and ten. A residual map is the difference between the estimate of the fake sky map and the actual fake sky map. There axe some noteworthy features in these maps. Small scale structure, present in the fake sky and in the residual maps, is gradually reduced in the residual maps during the first three iterations and is largely absent in the residual map at iteration 11. Between iterations 11 and 13, small-scale structure possessing a pattern different than the small-scale structure pattern at early iterations, begins to grow, and becomes the dominant residual pattern at iteration 25. Structure around the 150" scale remain constant during iteration. The edges of the maps also remain constant, because they are outside of the SHADES field and not measured. In this simulation the iterative algorithm diverges. The simulation described above is repeated. This time a single pointing of the SHADES measurement strategy is used to measure a fake sky populated by a single beam convolved point source with an amplitude of ten. Figure 3.2 depicts the residual maps for iterations 1 through 6 and iteration 25. As in the previous simulation, the residual flux oscillates and grows without bound. The initial map estimate might have an effect on algorithm convergence. To test this, the simulation used to produce Figure 3.1 is repeated, with an initial map 3The J C M T PSF is approximated as a Gaussian with a full-width half-maxmum equal to 14.5" and an amplitude of one. 4 A null map is a map whose pixel values are zero. Chapter 3. Iterative Reduction of SHADES Data 3 5 Figure 3.1: Residual maps of a simulation of the generalized iterative algorithm. The measurements are performed on a statistically plausible fake sky using the full SHADES measurement strategy. The measurements are noiseless. High spatial frequency structure is removed from the residuals in the early iterations. In later iterations the unstable spatial frequencies begin to grow. The fake sky map saturates at pixel values of zero and ten. Pixels are three arcseconds square. Chapter 3. Iterative Reduction of SHADES Data 36 20 40 60 80 100 20 40 60 80 100 2 0 4 0 6 0 8 0 1 0 0 Iteration 25 -40 -20 0 20 40 Figure 3.2: Residual maps for iterations 1 through 6 and 25 in a simulation of the generalized iterative algorithm applied to measurements taken of a single source. The measurements comprise a single pointing in the SHADES measurement strategy. Chapter 3. Iterative Reduction of SHADES Data 37 estimate equal to the fake sky map. Figure 3.3 shows the fake sky map, fake sky map estimates for iteration 1, iteration 11, iteration 21, and difference maps between the fake sky map estimate and the fake sky map for iteration 1, 11, and 21 for this simulation. The algorithm does not diverge from the perfect initial estimate. The initial map estimate is set to equal 1.1 times the fake sky map and the simulation is repeated. Figure 3.4 shows the residual maps for iteration 1, 15 and 25. By iteration 25 the unbounded mode in the final iterations of the simulation depicted in Figure 3.1 has re-appeared. This mode grows without bound and the algorithm diverges. It is interesting to note that the average pixel value in the residual maps of Figure 3.4 is around zero, whereas in Figure 3.1 the mean pixel value is between —0.9 and —0.4. The only change between the simulations producing Figure 3.1 and Figure 3.4 is the initial map estimate. Thus, to some extent, the initial map estimate affects the mean of the map pixels in this situation. At this stage, there is no noise, and all pixels involved in measurements are being updated. In this situation, assuming algorithm convergence, the edge pixels are well-defined, and no special treatment is required to determine the pixels at the map edges. In later chapters, not all pixels in each measurement will be updated and the measurements will be noisy. Under those circumstances, the edge behaviour is complicated and is addressed in section 4.1 on edge effects for a stable algorithm. The iterative algorithm instability might depend upon the size of the field mea-sured, the measurement strategy, the conditions near the edge of the measured field, the initial map estimate, and the sky being measured. The single pointing simula-tion has around 10000 pixels and the full SHADES simulation half a million. To investigate the effect of boundary conditions and to compute simulation pixel-pixel correlation matrices, the subsequent simulations are reduced to seven pixels. 3.2.2 Small Field Simulations The algorithm is known to converge when reducing single-difference measurements. A maximally measured seven pixel, single-difference simulation is performed. The fake sky map is sT = [ 0 0 0 1 0 0 0 Chapter 3. Iterative Reduction of SHADES Data 38 where the s*h element represents the flux from the ith pixel. The vector of measure-ments is obtained by, M = Ps where P is the pointing matrix and is given by: P = Thus, / - I 1 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 -1 1 0 V o 0 0 0 0 -1 1 / MT = [ ° 0 1 -1 0 0 (3.7) Note that there are six measurements and seven unknown pixel values, in the maxi-mum number of measurements on a one-dimensional, non-periodic, seven pixel patch of sky. This is a manifestation of the non-invertibility of PTP, and is due to the inability of the single-difference measurement scheme to measure a constant sky offset. Figure 3.5 shows the measured fake sky map, the fake sky estimate maps for iterations 45 to 50 for the iterative reduction of the measurements M, and the standard deviation and mean of the map pixel values as a function of iteration. The initial map estimate is the zero map. There are two estimates of the fake sky map in the middle plot of Figure 3.5. The algorithm oscillates between these two estimates with each iteration. The third plot in Figure 3.5 shows the mean and standard deviation of the map pixels constantly oscillating about some value. This is due to the oscillation in the fake sky estimate maps and continues with iteration indefinitely. The algorithm is marginally stable. The simulation is repeated with a single extra, "wrap-around", measurement. Chapter 3. Iterative Reduction of SHADES Data 39 The new pointing matrix is / - l l 0 0 0 0 - 1 1 0 0 0 0 - 1 1 0 0 0 0 -1 1 0 0 0 0 -1 0 0 0 0 0 \ 1 0 0 0 0 and the transpose of the measurement vector is 0 0 0 0 1 -1 0 0 \ 0 0 0 0 1 - 1 / (3.9) M1 0 0 1 - 1 0 0 0 (3.10) Figure 3.6 once again depicts the fake sky map, the fake sky estimate maps and a plot of the mean and standard deviation of the fake sky estimate maps with iteration. This time there is a single estimate of the fake sky map, correct up to a a constant offset. After a small number of osciallations the mean and standard deviation of the estimated map settle to a single value and the algorithm converges. The addition of the extra measurement changes both the simulation boundary conditions and the simulation measurement strategy, and increases algorithm stability. From P in equation 3.9 we can see that the extra measurement constrains the first and last pixel. It is difficult to determine if the effect is due to boundary conditions, or due to a change in measurement strategy since the boundary conditions refer to the way the edges of the map are measured. Besides the pointing matrix, P, the pixel-pixel correlation matrix, PTP, is very relevant in the solution to the linear system and is listed in Table 3.1.5 It is given here, for the pointing matrix, P, in equation 3.9, and referred to in later sections when comparing correlation matrices for different measurement strategies. The next step is to measure the same fake sky map with double-difference mea-surements. The first double-difference simulation uses the maximum number of measurements without wrap-around measurements and with a fixed single pixel 5See section 2.5 Chapter 3. Iterative Reduction of SHADES Data 40 Csd7 P i x 7 M e a s 2 -1 0 0 0 0 -1 -1 2 -1 0 0 0 0 0 -1 2 -1 0 0 0 0 0 -1 2 -1 0 0 0 0 0 -1 2 -1 0 0 0 0 0 -1 2 -1 V -1 0 0 0 0 -1 2 Table 3.1: Correlation Mat r ix : For a seven pixel, seven measurement, single-difference simulation. chop size. The resulting pointing matrix is P = -.5 1 -.5 0 0 0 0 0 -.5 1 -.5 0 0 0 0 0 -.5 1 -.5 0 0 0 0 0 -.5 1 -.5 0 0 0 0 0 -.5 1 -.5 (3.11) The measurements are M = 0 -.5 1 -.5 0 (3.12) The double-difference measurements do not measure an offset or a constant slope, and thus there are only five measurements. Figure 3.7 shows the map estimates for the last six iterations in a simulation with 25 iterations. The initial map estimate is the zero map. This time the algorithm is not marginally stable, as in the seven pixel, six measurment single-difference measurement simulation. Again on each iteration, the pixel values change sign and grow. The average also oscillates and grows, and the standard deviation exhibits something like exponential growth. The simulation is repeated with seven extra measurements added to the pointing matrix given by equation 3.11. The additional measurements include measurements with asymmetrical chop sizes, different chop sizes and "wrap-around" measurements. Figure 3.8 shows the simulation results. The estimate maps exhibit asymmetry, Chapter 3. Iterative Reduction of SHADES Data 41 and oscillate and grow with iteration. The asymmetry is due to the asymmetric measurements amongst the seven newly added measurements to the pointing matrix in equation 3.12. The pixel values grow without bound; however, the growth is slower than the growth depicted in Figure 3.7 for the double-difference, five measurement simulation. Adding two more measurements and repeating the simulation yields the results shown in Figure 3.9. The map is, once again, asymmetric and converged to within a constant of the correct fake sky map. The convergence is accurate to approximately 5% of the maximum pixel value in the fake sky map by iteration 100. The third plot in Figure 3.9 shows the map average converging within the first five iterations to a slightly positive value, and the map standard deviation decaying approximately exponentially. The stability of the iterative algorithm depends upon the measurement strategy. It is not clear what the dependence is and if other factors are contributing. It is also evident that there is a property of the single-difference pointing matrix, P, which is different from the double-difference pointing matrix and different between the 13 and 15 measurement, seven pixel double-difference simulations. A differ-ence is the fewer unmeasured modes in the single-difference PTP matrix, resulting in more equations in the same number of unknowns. It is instructive to look at the pixel-pixel correlation matrices for the different simulations.6 Both the pixel-pixel correlation matrix in Table 3.1 and the pixel-pixel correlation matrix in Table 3.3 correspond to converging iterative algorithm simulations. CSd 7 P i x 7 M e a s exhibits many zero off-diagonal terms and dd 7 P i x i 5 M e a s has generally a larger ratio of di-agonal entries to off-diagonal entries than C d d 7 P i x i 3 M e a s - I n the direct inversion of the least-squares method of data reduction, the presence of large off-diagonal terms is indicative of correlated pixels in the data-reduced map, and is undesirable. Con-vergence depends upon measurement strategy, but the understanding is qualitative and it is not known whether it is the measurements at the edge of the map, as in the 6 T o a i d d e s c r i p t i o n , t h r o u g h o u t the thesis p i x e l - p i x e l co r re la t i on m a t r i x refers t o PTP. T h i s is accu ra te , u p to a fac to r , w h e n t he t ime l i ne noise covar iance m a t r i x is d i a g o n a l , w i t h i d e n t i c a l ent r ies. Chapter 3. Iterative Reduction of SHADES Data 42 7 P i x 13 M e a s I 2.75 -1.5 0 -0.5 -1.5 2.75 -1 0.25 0 -1 2.75 -1 -0.5 0.25 -1 2.5 0.25 -0.5 0.25 -1 0 -0.5 -0.5 0.25 0.25 -0.5 0.25 -1 2.75 -1 0.5 -0.5 -0.5 -0.75 0 -0.5 -0.5 0.25 -1 2.75 -1 -1 \ 0.5 -0.5 -0.5 -0.75 -1 3.25 J V -1 Table 3.2: Correlation Mat r ix : Seven pixel, thirteen measurement, double-difference simulation. Cdd 7 P i x 1 5 M e a s — 3 -1.5 0 -0.25 0.25 0 -1.5 -1.5 2.75 -1 0.25 -0.5 -0.5 0.5 0 -1 3 -1 0.25 -0.25 -1 -0.25 0.25 -1 2.75 -1 0.25 -1 0.25 -0.5 0.25 -1 2.75 -1 -0.75 0 -0.5 -0.25 0.25 -1 3 -1.5 -1.5 0.5 -1 -1 -0.75 -1.5 5.25 V Table 3.3: Correlation Mat r ix : Seven pixel, fifteen measurement, double-difference simulation. Chapter 3. Iterative Reduction of SHADES Data 43 single-difference simulation in this section, or the initial map estimate, or the bulk of measurements at the centre of the field which determines the iterative algorithm convergence properties. Specifically, the following questions remain to be answered: 1. How does convergence depend on initial map estimate? 7 2. How does convergence depend on the unmeasured modes? 3. How does convergence depend on the boundary conditions?8 4. Do the small field simulations sufficiently model conditions in the actual data reduction to draw valid conclusions? 5. How does convergence depend on noise? Note that the SHADES measurement strategy is fixed and cannot be manipulated to achieve convergence. To determine the answer to the above questions by simulation requires the fixing of one parameter, say the initial map estimate, and then varying the others through all parameter space, and then repeating for different values of the parameter under investigation. This is a difficult and exhaustive procedure, which the discovery of a reformulation of the iterative deconvolution algorithm made largely unnecessary. It is the subject of Section 3.4 onwards. 3.2.3 Minimum Comparable Size The pointing matrix in the 15 measurement double-difference simulation of the last section uses measurements with all physically possible chop throws, "wrap around" measurements, and measurements with asymmetrically sized chop throws of a one-dimensional field. Despite this, the pixel-pixel correlation matrix in Table 3.3 is not very diagonal, and does not constitute an acceptable measurement strategy.9 Com-pare this to measurement strategies such as that used by WMAP, which introduces very small correlations between the pixels when iteratively deconvolving, and whose 7 "How" here means "Does? And if so, in what way?" 8 "Boundary conditions" refers to the tapering in the number of measurements of a region as the observation nears the edge. 9 U p to 20% correlation between pixels. Chapter 3. Iterative Reduction of SHADES Data 44 pixel-pixel correlation matrix is highly diagonal. WMAP does not perform every possible single difference measurement, yet achieves a relatively diagonal pixel-pixel covariance matrix. How is this accomplished? The number of unique measurements that can be made of a one-dimensional field is floor ( ^ ) Nm= £ (3.13) i=l where floor (x) returns the truncated integer representation of x. For a seven pixel map measured using double-differences there are nine unique measurements, not involving wrap-around or asymmetric measurements. Thus, the number of possible measurements of a partial-sky field is expected to be O (-Np\x). It is reasonable to expect that observations of large fields can achieve negligible correlations without exhausting all measurement possibilities. The covariance matrix of the SHADES simulation is too big to be easily con-structed and compared directly to the small simulation covariance matrices. Instead, a figure of merit is computed, which is the ratio of the central pixel's diagonal entry in the pixel-pixel correlation matrix to the sum of the absolute entries in the row. [2] This figure of merit has a maximal value of one and is applicable where the mea-surement strategy is uniform in the centre of the observed field. The figure of merit of the SHADES measurement strategy for a single pointing is 0.376, for the single-difference, seven pixel, seven measurement simulation 0.5, for the double-difference, seven pixel, thirteen measurement simulation 0.417 and for the double-difference, seven pixel, fifteen measurement simulation the figure of merit is 0.423. Calculat-ing the figure of merit for the entire SHADES strategy, though requiring much less memory than the full SHADES correlation matrix, requires the full pointing matrix (minus the zero entries) which is still prohibitive. The iterative deconvolution algo-rithm applied to the full SHADES data set in simulation is slightly unstable and if the figure of merit is a good indicator of algorithm stability then we expect a figure of merit for the full SHADES simulation between 0.417 and 0.423. Chapter 3. Iterative Reduction of SHADES Data 45 3.3 Iterative Inversion As Newton-Raphson Root Finding A method of determining the roots of a non-linear function of one variable is to make a guess as to the location of the root, compute the tangent of the function at this location, and then compute where the root would be located if the function were the tangent. The estimate of the root is the value of the function evaluated at the new estimated root location. The process is repeated until convergence. The rate of convergence in this algorithm depends upon the initial guess and the degree of non-linearity of the function. Divergence typically occurs when the slope of the tangent line to the curve at the estimated root location is zero. [13] Consider the equation Px = M (3.14) where P is the pointing matrix, x is the map of pixels, and M is the corresponding vector of measurements. Form an error vector E, such that E(x) = M-Px (3.15) The best estimate map will be the map corresponding to the roots of E. That is, the position in Np dimensional space corresponding to an error vector of minimal magnitude. Consider a small deviation, Sx, from the current map estimate, x, then NpiK BE Ei (x + Sx) = Ei (x) + V ——+ O {Sx2) (3.16) j=l d X j where Ei is the 2 t h element of E. Let x* = x + Sx (3.17) We desire x' to be the root of E (x). Neglect the higher order terms in 3.16, and set E (a?) = 0. Solve for 8x and obtain 6x = -J~lE = -P (3.18) where 13 dxj (3.19) Chapter 3. Iterative Reduction of SHADES Data 46 and P is the pointing matrix. The Jacobian, J, has dimensions Nm x iVpix and is not invertible. The least-squares optimal update is Sx = {PTP)~1 [PT (M - Ps)] (3.20) A couple of notes. First, PTP is not invertible, is large, and its inversion is the original problem, ie. 8x = x - (PTP)~1 PTM (3.21) Thus, the new map estimate is £ n e w = 2?old + SX= {PTP)~1 PTM (3.22) which is the least squares optimal solution to equation 3.14, is the topic of section 2.5.3, and is whose computation was hoped to be carried out by iteratively applying equation 3.1. Secondly, the error function is linear. That is, the higher order terms in equation 3.16 are zero, and convergence is expected in a single step, independent of the initial map estimate. The approach is abandoned, as it is simply a restatement of the original problem to which iterative deconvolution is addressed. 3.4 Formulation of Iterative Inversion as a Difference Equation The following formulation is motivated by a desire to mathematically describe the map update in a single equation. 3.4.1 Formula t ion Let the error between the data and the measurement of the map estimate be En = M - Pxn (3.23) where M is a vector of measurements, P is the pointing matrix and xn is the n t h map estimate. Note that the error on the jth measurement is El = Mj - g- • xn (3.24) Chapter 3. Iterative Reduction of SHADES Data 47 where g~j is a vector of gains corresponding to the j t h measurement. Now the update for the 2 t h pixel, given by equation 3.1, is 9i • 9i + x n - 1 The numerator in the fraction of the above equation is an inner product, ^ g i E ^ g - i - E (3.25) (3.26) and hence, 9i-E g% • 9i F + xn-l (3.27) Thus, each pixel update is an inner product between a gain vector for the 2 t h pixel and a vector of errors, scaled by the inverse square magnitude of the gain vector plus the previous pixel value. In a sense, to within a scale factor, the inner product between the gain vector and the error vector is the act of measuring the error contributed to the 2 t h pixel by the difference between the actual measurements and the measurements of the current map estimate. The last equation can be extended to give the update of the entire map. 1 \ 90-90 91-91 / go ^ 91 \ 9NP J En-1 + ^ n - l (3.28) 9Np-9NP ) Noting that 9i is equal to the tranpose of the pointing matrix P, and sub-\9NP J stituting equation 3.4 for En-i, equation 3.9 can be simplified with the following Chapter 3. Iterative Reduction of SHADES Data 48 algebra. — go-go gi-gi PT • ( M - P £ „ _ i ) + f„_i(3.29) — \ §Np-gNp J (I - DPTP) • xn-i + £ > P T • M = A • f „ _ i + 5 (3.30) (3.31) where D = l go-go 9i-Si 9Wp -9Afp / A = 7 - DPTP and P = D P T • M . Equation 3.31 is a difference equation giving the n t h estimate of the map. 3.4.2 Solution The solution is standard and is achieved by mathematical induction. From equation 3.31 the first three map estimates are x\ = A • xo + B x2 = A • xi + B = A2 • x0 + AB + B xz = A • x2 + B = A (A2 • x0 + AB + B) + B = A3 • x0 + A2B + AB + B (3.32) (3.33) (3.34) By induction the solution is n - l xn = An • xo + A T B (3.35) r=0 Chapter 3. Iterative Reduction of SHADES Data 49 Equation 3.35 gives the estimate of the map produced by the n t h iteration of the generalized Wright algorithm. There are three easy tests of equation 3.35. They are: 1. Solution behaviour with a correct initial map estimate. 2. Solution behaviour with a diagonal pixel-pixel correlation matrix. 3. Solution behaviour as the number of iterations tends to infinity. The checks are presented in the following three sections. 3.4.3 Solut ion Check #1 If XQ is correct, then The n t h map estimate should equal the initial map estimate if this is correct. The n t h map estimate becomes: P-x0 = M (3.36) n - 1 (3.37) and DPTP = I-A (3.38) then n - 1 (3.39) n - 1 (3.40) An • x0 + (J - A ) • x0 (3.41) x0 (3.42) as expected. Chapter 3. Iterative Reduction of SHADES Data 50 3.4.4 Solution Check #2 In section 3.3.2, it is noted that the stability of the generalized iterative data re-duction algorithm depends, at least in some way, upon the measurement strategy. A perfect measurement strategy results in a diagonal pixel-pixel correlation matrix. Furthermore, the diagonal elements of D in the generalized Wright algorithm scheme are equal to the inverse of the diagonal elements of the pixel-pixel correlation matrix. Thus, with a perfect measurement strategy, DPTP is equal to the identity matrix, and A equals zero. With A set to zero in equation 3.35, a single term survives: xn = DPT • M (3.43) Noting that DPTP = I: DPT = P+ (3.44) Where P+ is the pseudo-inverse10 of P. Thus, xn = P+ • M (3.45) When the sky is measured perfectly, the iterative algorithm gives the solution of P • x = M on any iteration up to a linear combination of the singular eigenvectors of P. 3.4.5 Solution Check #3 Consider the limit as n -> oo. It is desired that l i m ^ o o X n = P+ • M (3.46) Starting with the summation term in equation 3.35, note that n n+1 ^ A r - ^ A r = /-A" + 1 (3.47) r=0 r = l n Ar [I - A] = I - An+1 (3.48) r=0 n - 1 and Y,Ar = [I- An][I-A]+ (3.49) r=0 °See, for i ns tance , [21]. Chapter 3. Iterative Reduction of SHADES Data 51 Equation 3.35 becomes xn = An • x0 + [I - An] [I - A}+ (3.50) Diagonalizing the Np by Np matrix A obtain A = vXv'1 (3.51) and A2 = vXv^vXv'1 = vX2v~1 (3.52) By induction A n = vXnv~l (3.53) where v is a matrix composed of the eigenvectors of A, and A is a diagonal matrix whose diagonal entries are the eigenvalues of A. Note that the diagonalization is impossible if A does not have iVpix linearly independent eigenvectors. A is a nearly symmetric, band-diagonal, invertible matrix, which has zero entries along the diagonal. A rigorous analysis can be performed using the Jordan Form. The Jordan form analysis yields the same results with a slight modification of the algorithm convergence when the absolute value of the eigenvalues approachs one. 1 1 Now xn = vX^'1 • x0 + [I- vXnv~l] [I - A}+ B (3.54) and lim ^ o o f n = vX°°v-1 -x0+ [I- vX^v'1] [I - A}+ B (3.55) If the absolute value of the eigenvalues of A are less than one, then we see that there is no dependence on the intial map estimate So and l i n w ^ = [I-A]+B (3.56) = [DPTP]+ DPT • M (3.57) = P+ (DPT)+ DPT • M (3.58) l i m ^ o o f n = P+-M (3.59) which is the desired solution up to a linear combination of singular eigenvectors of DPTP. If any of the absolute values of the eigenvalues of A are greater than one then xn will grow without bound with increasing n, and is unstable. 1 1 See Appendix A. Chapter 3. Iterative Reduction of SHADES Data 52 -1.73 -1.33 -0.346 0.500 0.909 1.00 1.00 Table 3.4: Eigenvalues of the characteristic matrix for the seven pixel, five double-difference measurement simulation. -0.522 0.298 -0.308 0.295 -0.308 0.298 -0.522 -0.606 0.303 -0.202 0.000 0.202 -0.303 0.606 0.617 -0.177 -0.140 0.369 -0.140 -0.177 0.617 -0.632 0.000 0.316 0.000 -0.316 0.000 0.632 0.627 0.215 -0.141 -0.285 -0.141 0.215 0.627 0.378 0.378 0.378 0.378 0.378 0.378 0.378 0.378 0.378 0.378 0.378 0.378 0.378 0.378 Table 3.5: Eigenvectors of the characteristic matrix in the seven pixel, five double-difference measurement simulation. Each row represents a separate eigen-vector. The first row is the eigenvector corresponding to the left-most eigenvalue in Table 3.4. For clarity, the eigenvector components have been rounded to three significant figures. The eigenvectors are not lin-early independent. 3.4.6 Difference Equa t ion Solut ion S imula t ion The generalized Wright algorithm is simulated in a seven pixel simulation with five double-difference measurements. The 25th iteration map estimate is compared with the map estimate computed using equation 3.35. The maps are identical up to machine precision. The eigenvalues and eigenvectors of the characteristic matrix A are given in Tables 3.4 and 3.5, respectively. There are two absolute eigenvalues greater than one, and two eigenvalues that are identically one. The eigenvectors corresponding to the eigenvalues of one are the same and correspond to the offset of the map to within floating point precision. The eigenvectors are not linearly independent and the characteristic matrix is not diagonalizable. Table 3.6 and 3.7 list the eigenvalues and eigenvectors for a seven pixel, 15 double-difference measurement simulation. All of the absolute eigenvalues are less Chapter 3. Iterative Reduction of SHADES Data 53 1.00 -0.987 -0.685 0.447 0.050 -0.060 0.235 Table 3.6: Eigenvalues of the characteristic matrix for the seven pixel, 15 double-difference measurement simulation. 0.378 0.378 0.378 0.378 0.378 0.378 0.378 0.491 -0.573 0.401 -0.252 0.333 0.080 -0.298 0.174 -0.079 -0.195 0.553 -0.444 0.563 -0.325 0.509 0.606 0.105 -0.367 -0.415 -0.174 -0.159 0.438 -0.500 -0.293 -0.153 -0.454 0.006 0.493 0.403 0.064 -0.578 0.354 0.428 -0.426 -0.100 -0.003 0.159 -0.487 -0.550 0.329 0.571 -0.014 Table 3.7: Eigenvectors of the characteristic matrix for the seven pixel, 15 double-difference measurement simulation. Each row is a separate eigenvector; the eigenvector in the top row corresponds to the left most eigenvalue in Table 3.6. Chapter 3. Iterative Reduction of SHADES Data 54 than or equal to one, and the simulation is stable (See Figure 3.8). The solution is correct up to a map offset which is the eigenvector corresponding to the eigenvalue of one. Examining the last plot in Figure 3.8, we see that the average of the map exhibits some dynamics at early iterations and settles to a positive constant value. Keeping this in mind we look at equation 3.35 in the context of the seven pixel, 15 double-difference measurement simulation. In this simulation, the initial map estimate is the zero map and equation 3.35 becomes n - 1 sn = J2ATB (3-6°) r = 0 ie. xn = P+M - AnP+M (3.61) The first term in equation 3.61 is the desired solution and the second term we expect to die off with iteration. Unfortunately the eigenvectors in Table 3.7 are not orthogonal (as expected, since A is not symmetric) and are not linearly independent, and thus A is not diagonalizable. The behaviour of the map pixel average can be computed using the Jordan form, but it is preferable to analyze a simpler case. Nonetheless, qualitatively, with iteration the contributions to the map and to the mean of the map due to the eigenvectors with absolute eigenvalues less than one decrease, and we see the map mean fluctuate at early times and then stabilize. The reader is again referred to Appendix A for a discussion of the Jordan form. 3.4.7 Stab i l i ty As previously mentioned, the absolute eigenvalues of the A matrix in equation 3.35 are dominant in determining the stability of the iterative algorithm. What type of measurement strategy results in absolute eigenvalues greater than one? The answer is measurement strategies that result in a wider band in the band diagonal PTP matrix. The A matrix, in equation 3.35 is zero along its diagonal. Heuristically, the larger off-diagonal bands in A force the eigenvectors of A to contain greater structure. The physical characteristics of measurement strategies resulting in narrow band-diagonal A matrices is addressed in the section on measurement strategy. In the next chapter the rate of convergence of simulations with the same mea-surement strategy will vary depending upon the fake sky that is measured. This is Chapter 3. Iterative Reduction of SHADES Data 55 expected; consider the difference between the xn and xn-i maps: = ^2 ArDPTM - ^2 ArDPTM r=0 r = 0 Ax„ = An~lDPTM Let A be diagonalizable so that then A = vXv 1 (3.62) (3.63) (3.64) (3.65) (3.66) Axn = v\n-1v-1DPTM Convergence will vary depending on the degree of orthogonality between DPTM and the inverted eigenvectors corresponding to large absolute eigenvalues of A. Thus, one expects convergence to depend upon the measurements to a small extent. It might be possible that absolute eigenvalues greater than one correspond to attempting to invert a non-invertible matrix without constraining the unmeasured modes of PTP. Fixing a pixel, in the context of the iterative algorithm, means to never update a pixel involved in measurements. In this way the pixel is constant with iteration and stays fixed. A six pixel fake sky map is measured with the following pointing matrix. / -.5 1 -.5 0 0 0 -.5 1 -.5 0 P = 0 0 \ -.5 -.5 0 -.5 (3.67) 0 0 \ 0 0 0 -.5 1 The pixel fixing can be performed by modifing the PT matrix such that the first and last rows are zero, i.e. / 0 1 -.5 0 0 V o P1 = 0 0 0 \ -.5 0 0 1 -.5 0 -.5 1 -.5 0 -.5 1 0 0 (3.68) 0 / Chapter 3. Iterative Reduction of SHADES Data 56 The eigenvalues of the characteristic matrix A — I - DPTP are listed in Table 3.8. There is an absolute eigenvalue greater than one. -1.3 0.64 1.0 -0.31 0.97 1.0 Table 3.8: Eigenvalues of the characteristic matrix for the six pixel, four double-difference measurement simulation with two fixed pixels. There is an absolute eigenvalue greater than one. Perhaps setting, otherwise unmeasured pixels to a value will make the algorithm stable. In this case the pointing matrix, corresponding to two absolute measure-ments and two double-difference measurements, is P = (3.69) f \ 0 0 -.5 0 0 \0 0 PT really equals the transpose of the pointing matrix and the eigenvalues of the characteristic A matrix are listed in Table 3.9. There are absolute eigenvalues greater 0 0 1 -.5 -.5 1 0 0 0 0 N 0 0 -.5 0 0 1 0 1 -1.8 -0.2 1 0 Table 3.9: Eigenvalues of the characteristic matrix in the simulation with six pix-els, two double-difference measurements and two absolute measurements. There is an absolute eigenvalue greater than one. than one. Note that the fixed pixels in this case are not measured by the double-difference measurements. Next add two more double-difference measurements linking the absolute mea-Chapter 3. Iterative Reduction of SHADES Data 57 surements. The resulting pointing matrix is: / 1 0 0 0 0 0 -.5 1 -.5 0 0 0 0 -.5 1 -.5 0 0 0 0 -.5 1 -.5 0 0 0 0 -.5 1 -.5 \ 0 0 0 0 0 1 ) (3.70) Again, PT really equals the transpose of the pointing matrix and the eigenvalues of the characteristic A matrix are listed in Table 3.10. There is an absolute eigenvalue -1.4 -.62 .72 .31 .03 0.98 Table 3.10: Eigenvalues of the characteristic matrix for the six pixel, four double-difference measurements and two absolute measurements simulation. There is an absolute eigenvalue greater than one. greater than one. As a prelude to the discussion of the central update algorithm in the section on algorithm stabilization and in chapter four, consider the following modified PT matrix: / 0 0 0 0 \ 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 \0 0 0 0/ The resulting eigenvalues are listed in Table 3.11. Pl = (3.71) 3.4.8 Noise Proper t ies If the prior probability distribrution is uniformly distributed and non-zero over a reasonable range of pixel values, and the measurement noise is Gaussian , the width Chapter 3. Iterative Reduction of SHADES Data 58 -1.0 0.81 -0.81 0.31 -0.31 1.0 Table 3.11: Eigenvalues of the characteristic matrix for the six pixel, four double-difference measurement and two absolute measurement. There are no absolute eigenvalues greater than one. of the likelihood function is given by the pixel-pixel noise covariance matrix. The pixel-pixel noise covariance matrix, C, is computed as follows: C = {xnxl) - (xn) {xn)T (3.72) To compute this, first the expected map estimate is determined using equation 3.33. Let the absolute eigenvalues of A be less than one, then the nth expected map estimate is, 1 2 <£„> = ((DPTP)+DPT -(M0 + n)) (3.73) (xn) = (P+{DPT)+DPT -(M0 + n)) (3.74) {£„) = (P+• M 0 + P +• n) (3.75) (xn) = P+-M0 (3.76) Substituting equation 3.75 for the second term in equation 3.72, C becomes C = ( P + ( M 0 + n)(M0 + n ) T ( P + ) T ) - P+M0M^P+T (3.77) since the expectation of n is zero, and the last term cancels the first term in the product, the pixel-pixel noise covariance matrix is C = {P+rtnTP+T) = P+NP+T (3.78) C = P+NP+T (3.79) C = (PTN-1P)+ (3.80) If PTP is invertible, then this is equivalent to, C = (PTN-1P)~1 (3.81) 2 T h i s is cor rec t u p to the n o n - i n v e r t e d subspace of the ma t r i ces b e i n g pseudo - i nve r t ed . Chapter 3. Iterative Reduction of SHADES Data 59 where N is the timestream noise covariance matrix. Note that C is equivalent to the pixel-pixel noise covariance matrix obtained from the generalized least-squares formulation. [9] 3.5 Algorithm Stabilization The measurement strategy of SHADES is fixed, and if no modifications of the itera-tive algorithm (equation 3.35) are made, the algorithm does not converge. The goal is to modify the algorithm in a fashion that makes the algorithm stable while not ap-preciably affecting the quality of the solution. The following algorithm modifications are attempted: 1. Unstable eigenvector identification and removal. 2. Modifying the pixel size at the boundaries of the map. 3. Decoupling of measurements and pixel updates. 3.5.1 Unstab le Eigenvector Remova l The idea is to determine the eigenvectors and associated eigenvalues of A in equa-tion 3.35 which have absolute eigenvalues greater than one, and subtract them, appropriately scaled, from the final map. The trick is to find the eigenvalues and eigenvectors without forming the matrix. Even if the matrix is formed, computation of the eigenvalues and eigenvectors is an 0(n3) operation.[13] A potential probe of unstable eigenvectors is obtained by setting the measurement vector, M, in equation 3.35 to zero. Then, only the first term in equation 3.35 is computed. After many iterations, xn = An • f 0 = v^vTl • fo (3-82) i where Vi is a matrix whose columns are the eigenvectors with absolute eigenvalues greater than one. If SQ is guessed to be an eigenvector, then the nth estimate map is xn = A£a?o (3.83) Chapter 3. Iterative Reduction of SHADES Data 60 where Ae is the eigenvalue associated with an unstable eigenvector of A, namely the guessed eigenvector xo-ls The ratio of the output maps from two successive iterations gives the associated eigenvalue. Assuming knowledge of the unstable eigenvectors and eigenvalues, the least-squares optimal map can be recovered. Let u = [I - A]'1 DPTM (3.84) Equation 3.35, with a zero-map initial guess becomes xn = (l- An+l) u (3.85) Consider xn - xn-i = (-An+1 + An) u (3.86) Diagonalize A in the limit of large n: An+1 = veA2+1vT (3.87) where v is a matrix whose columns are populated with the unbounded eigenvectors of A, A is a diagonal matrix populated with eigenvalues and the subscript e stands for explode. That is, all of the eigenvectors, corresponding to absolute eigenvalues less than one, are negligible in the decomposition.14 Let ci = vie-u (3.88) be the the inner-product of u with the ith unstable eigenvector. Then, xn - = £ 4 (A* ) n * (1 - XI) (3.89) i and the least-squares optimal map is n xn = ^ 2 ArDPTM + ^ 4K+lci (3-90) r = 0 i Thus, the task is to determine the unstable eigenvectors, their associated eigen-values and the inner-products of u with the unstable eigenvectors. If these are 1 3 Guessing an eigenvector is somewhat fortunate. 1 4 Again, all of this assumes linearly independent eigenvectors. The formulation can be generalized using the Jordan form. \ Chapter 3. Iterative Reduction of SHADES Data 61 -1.73 -1.33 -0.346 0.500 0.909 1.00 1.00 Table 3.12: Eigenvalues of the characteristic matrix for the seven pixel, five double-difference measurement simulation. known, the least-squares optimal map is given by equation 3.90. The characteristic matrix, A, whose eigenvectors we wish to determine is iVpixel x Npjxe\. Determining the eigenvectors is typically an O [n3] computation, with the number of operations going to O [n2] per eigenvalue15 for a Hessenberg16 matrix. Thus, if there were a way to determine which eigenvalues one would determine prior to applying the algorithm, it might be possible to determine the eigenvectors of A. [13] Future work might investigate non-standard techniques with sparse matrices; however, with the convergence of the central update algorithm given below, this method is abandoned. 3.5.2 Boundary P i x e l Mod i f i ca t ion Changing the size of pixels at the edge of the map as a means of stabilizing the algorithm is investigated. Consider a seven pixel, one dimensional simulation of the iterative data reduction algorithm with the edge pixels enlarged to twice the size of the central pixels. The normal pointing matrix is / -.5 1 -.5 0 0 0 0 0 -.5 1 -.5 0 0 0 0 0 -.5 1 -.5 0 0 0 0 0 -.5 1 -.5 0 V 0 0 0 0 -.5 1 -.5 / (3.91) The associated eigenvalues of the characteristic matrix are given in Table 3.10, and reproduced here for convenience. With doubly enlarged edge pixels, the pointing 1 5 I f a l l e igenvalues are des i red , t h i s is aga in , a n O [ra3] c o m p u t a t i o n . 1 6 A Hessenbe rg m a t r i x c a n be f o r m e d ou t o f a n a r b i t r a r y m a t r i x t h r o u g h a series o f H o u s e h o l d e r t r ans fo rma t i ons . Chapter 3. Iterative Reduction of SHADES Data 62 -1.57 -0.797 0.406 1.00 0.953 Table 3.13: Characteristic eigenvalues for the edge-pixel enlarged pointing ma-trix. matrix becomes .5 -.5 0 0 0 \ -.5 1 -.5 0 0 0 -.5 1 -.5 0 0 0 -.5 1 -.5 V 0 0 0 -.5 .5 / The eigenvalues of the characteristic matrix A are given in Table 3.10. By com-paring the eigenvalues of the original characteristic matrix with the eigenvalues of the edge-pixel size increased characteristic matrix, it is apparent that the absolute value of the largest absolute eigenvalue in the two cases has decreased between the original characteristic matrix and the one with larger edge pixels, causing the itera-tive algorithm to become more stable. It is possible that algorithm stability might be accomplished by increasing the size of the edge pixels. In a sense, this reduces the number of unknowns while maintaining the same number of equations. Besides stability, enlarging the edge pixels might yield an improved estimate of the map by increasing the number of measurements involving an edge pixel. This is limited, however, since the resolution at the edges deteriorates as edge pixel size increases. That said, it is possible that a modest edge pixel enlargement would stabilize the full SHADES simulation of section 3.2.1. In light of the stabilizing method of the next section this method is not pursued further. 3.5.3 Measurement A n d P i x e l Upda te Decoupl ing A stable modification to Wright's algorithm was found by Colin BorysPersonal-Comm. The modification involves the D and PT matrices in equation 3.35 All oc-curences of —0.5 in PT are changed to zero, and any entries in D where, = oo, are set to zero. This modification will henceforth be called the "central update algo-rithm". The pointing matrix is left intact. Thus, the pointing matrix and the pixel Chapter 3. Iterative Reduction of SHADES Data 63 -0.866 1.00 0.866 0.000 -0.500 0.500 1.00 Table 3.14: Characteristic matrix eigenvalues for the central update algorithm in the seven pixel, five double-difference measurement simulation. All of the absolute eigenvalues are less than or equal to one. update matrix are decoupled. The result is a band diagonal A with a thinner band and smaller absolute eigenvalues. The following eigenvalue computation of A for a seven pixel double-difference pointing matrix with modified transpose demonstrates the effect. This method works because the n t h map estimate given by equation 3.55 does not depend on the D or the PT matrices. Thus, when considering only stability, they can be modified in any fashion, so long as the absolute eigenvalues of the resulting A matrix are less than one. The central update algorithm does not modify the pixel values at the edge of a map. An edge pixel that is not involved in the central part of a double-difference measurement is not updated. Thus, the initial map estimate at the edges is fixed, and any sky flux at the edges manifests itself as erroneous pixel values near the edge of the map. The effect decreases with distance from the edge and increases with the amount of flux at the edges of the field. In the next chapter performance of the central update algorithm when reducing SHADES data is characterized by simulation. Chapter 3. Iterative Reduction of SHADES Data 6 4 Figure 3.3: Top: Fake sky map. Middle: Fake sky map estimates for iteration 1, 11 and 21. Bottom row: Residual maps for iteration 1, 11 and 21. The colour bar applies to the residual maps only. Chapter 3. Iterative Reduction of SHADES Data 65 I t e r a t i o n 1 i \ \ I I t e r a t i o n 15 l1 I I t e r a t i o n 25 200 400 600 400 600 i ... 200 400 600 -i 1 1 -0 .1 - 0 . 0 5 0 0.05 0.1 Figure 3.4: Residual maps for iteration 1, 15 and 25. Note the reduction of high-frequency spatial structure between iteration 1 and 15 and the re-appearance of high-frequency structure in iteration 25. The algorithm slowly diverges. Chapter 3. Iterative Reduction of SHADES Data 66 Figure 3.5: Single-Difference: 7 Pixel Simulation, 6 Measurements Chapter 3. Iterative Reduction of SHADES Data 67 Figure 3.6: Single-Difference: 7 Pixel Simulation, 7 Measurements Chapter 3. Iterative Reduction of SHADES Data 68 Figure 3.7: Double-Difference: 7 Pixel Simulation, 5 Measurements Chapter 3. Iterative Reduction of SHADES Data 69 Figure 3.8: Double-Difference: 7 Pixel Simulation, 13 Measurements Chapter 3. Iterative Reduction of SHADES Data 70 Fake Sky Map 1.2 - 1 • • . . . . _ 1.0 E- + 0.8 r -z 0.8 0.4 z~ -_ 0.2 - -0.0 + + + + + + -j -0.2 - -0 1 2 3 Map: Iteration 20 To Iterotion 25 4 6 1.5 -1.0 =- -0.5 0.0 — 111111111 -0.5 S -1.0 + Iter. 20 - Iter 21 . Iter. 22 <>!ter. 23 .A.lter 24 Iter. 25 z 0 i 2 3 PiK*i Map: Iteration 95 To iterotion 100 4 5 6 1.0 0.5 -E 0.0 A — • • 9t — —-B-— • 0.5 1.0 + Iter. 95 • Iter 96 . Iter. 97 0 Iter. 98 .A l t e r 99 Iter. 100 : 0 1 2 3 4 5 6 Map Mean and Mop Standard Deviotion vs. Iteration Number + Map Standard Deviation L * Map Average + 0 .200000 H^*M^-H-H*H*H-H^i -H-H-H-H4i^^ H 1 1 H IH I H1111111111 B Iteration Number Figure 3.9: Double-Difference: 7 Pixel Simulation, 15 Measurements Chapter 4. Central Update Algorithm 71 Chapter 4 Central Update Algorithm The central update algorithm was introduced in section 3.5.3. In this chapter the edge effects, the algorithm behaviour when reducing measurements with additive Gaussian white noise, and the algorithm behaviour when reducing measurements with more realistic noise, are investigated by simulation with the SHADES mea-surement strategy. The application of the algorithm to real SHADES data, and subsequent source detection, are described in section 5.3. 4.1 Edge Effects To investigate the edge effects, a fake sky with a single, centred source is measured with the full SHADES measurement strategy and reduced using the central update algorithm. No noise has been added to the measurements. Figure 4.1 is a plot of the standard deviation of the residual map pixels as a function of iteration. Only pixels 150" from the map edge are included in the calculation. The standard deviation of the residual map pixels is 0.0008 mjy at iteration 101. At iteration 1001 the standard deviation of the residual map pixels is less than 0.0001 mjy Because there is no noise in this simulation, it is expected that the standard deviation of the residual pixels will asymptotically approach zero. Figure 4.2 shows residual maps for various, arbitrarily chosen, iterations. The central update algorithm converges and the edge effect, if there is one, is negligible by iteration 75. The simulation is repeated with a realistic fake sky map whose periphery is populated with sources drawn from a realistic source count model. The periphery is 150" wide (50 pixels). This is the same statistical source number count model used for various simulations of the generalized algorithm in Chapter 3. The purpose Chapter 4. Central Update Algorithm 72 Residual Map Standard Deviation 0.005 0.004 0.003 •5 0.002 0.001 0.000 400 600 Iteration Number 1000 Figure 4.1: Standard deviation of the residual map pixels as a function of iteration for a central update algorithm data reduction simulation. The simu-lated data consist of measurements of a single, 10 mjy, centred source, using the full SHADES measurement strategy. The computation ne-glects pixels within 150" from the map edge. The standard deviation of the residual map pixels is 0.0008 mjy at iteration 101. At iteration 1001 the standard deviation of the residual map pixels is less than 0.0001 mjy. Chapter 4. Central Update Algorithm 73 R e s i d u a l 1 R e s i d u a l 2 R e s i d u a l 5 I B •11 •1 1 •11 M l • 1 ! 600 -I 400 -J 200 • 200 400 600 200 400 600 200 400 600 J I L - 0 . 4 - 0 . 2 0 0 . 2 0 . 4 R e s i d u a l 11 R e s i d u a l 21 R e s i d u a l 75 600 — •ft a n 400 i : 200 — * 1 •I • n • 1 200 400 600 T 200 400 600 _ i _ i i r 200 400 600 "1 1 - 0 . 0 1 - 0 . 0 0 5 0 0 . 0 0 5 0 . 0 1 Figure 4.2: Residual maps estimated on various iterations of a simulation of the central update algorithm using the full SHADES measurement strategy. The fake sky map is a single source, centred in the map, with a peak flux of ten. The pixel values are in mJy and the contrast has been stretched on the lower three panels. The first iteration map is a measurement map. Note the convergence and the absence of a visible effect caused by the edges of the map. The maximum absolute pixel value of the residual map on the 75th iteration is less than a tenth of a percent of the source amplitude. Chapter 4. Central Update Algorithm 74 of this simulation is to determine the edge-effect produced by the central update algorithm in a realistic SHADES simulation. Inside the periphery the central portion of the fake sky map is set to the null map. Figure 4.3 is a plot of the standard deviation of the residual map pixels as a function of iteration. Only pixels inside the residual map periphery are used in the computation of the standard deviation of the residual map pixels. At iteration 401, the residual map pixel standard deviation is within 0.02 mJy of the converged value of 0.07mJy. The steps in the curve between iterations 401 and 801 are due to limited precision when reading and writing data to file.1 Note the maximum residual map pixel standard deviation of 0.185 mJy near iteration 60. At iteration one, the initial map estimate is the null map, and all pixels used in the computation of the standard deviation are equal to the fake sky map. Thus the standard deviation of the residual map pixels is zero, and rises with iteration as pixels inside the periphery are modified due to flux in the chop positions of the measurements centred upon pixels at the edge of the periphery. Note that convergence occurs much faster with measurements of the single centred source than with measurements of the statistically plausible edge-only fake sky, even though the measurement strategy, and thus the eigenvectors and eigenvalues of the characteristic matrix, are identical. This is due to the coupling described in section 3.4.7 between the measurements and the characteristic matrix. Figure 4.4 shows the fake sky map and residual maps for several iterations. With iteration, edge-effect flux propagates into the centre of the residual map when re-ducing the SHADES data with the central update algorithm. Note that the residual map mean grows increasingly negative with iteration and reduces the residual map pixel standard deviation in the higher iterations with respect to the maximum of the residual pixel standard deviation near iteration 60. At iteration 101, near the peak of the residual map standard deviation, the flux residual due to the edge ef-fect ranges from 0.25 mJy 360" from the map edge, to 0.5 mJy, 80" from the map edge. Figure 4.5 is the 801s* residual map. At iteration 801 the flux residual is less than 0.15 mJy 150" from the map edge. The 50 pixels, 150", map edge is ap-proximately the size of the outermost part of the map which is not measured in the 1 Author oversight. Chapter 4. Central Update Algorithm 75 Residuol Mop Stondord Deviotion Vs Iteration 0.201 1 1 ' I ' I ' 1 ' I 0.05 [- -0.001 i i i I ) I i i i I i i i 0 200 400 600 SCO Iteration Number Figure 4.3: Standard deviation of the residual map pixels as a function of iteration during a central update algorithm simulation. The simulated measure-ments are of an edge only, statistically plausible fake sky map and axe noiseless. Pixels within 150" of the map edge are not used in the com-putation. At iteration 401, the residual map pixel standard deviation is within 0.02 mJy of the converged value of 0.07 mJy. The steps in the curve between iterations 401 and 801 are due to limited precision when reading and writing data to file. Note the maximum residual map pixel standard deviation of 0.185 mJy near iteration 60. At iteration one, the initial map estimate is the null map, and all pixels used in the compu-tation of the standard deviation are equal to the fake sky map. Thus the standard deviation of the residual map pixels is zero, and rises with iteration as pixels inside the periphery are modified due to flux in the chop positions of the measurements centred upon pixels at the edge of the periphery. Chapter 4. Central Update Algorithm 76 Fake Sky Map Residual 21 200 400 600 200 400 600 Residual 301 Figure 4.4: An edge-only fake sky map and residual maps for a noiseless, full SHADES simulation of the central update algorithm. The fake sky map is populated with sources drawn from a realistic source count model. The edges are 150" wide. The pixel values are in mJy. Note the growth of the edge error. The maximum residual pixel standard deviation oc-curs around iteration 60. Note how the mean of the residual map pixels grows increasingly negative, reducing the standard deviation. At itera-tion 101, the absolute pixel values are less than two and a half percent of a 10 mJy source, 360" from the map edge. The white lines are spaced 200 pixels apart, each pixel is 3" wide. Chapter 4. Central Update Algorithm 77 R e s i d u a l Map 8 0 1 J I L -0.0007 -0.00065 -0.0006 -0.00055 -0.0005 -0.00045 (Janskies) Figure 4.5: The 801st residual map for the statistically plausible, edge only, central update algorithm simulation. The residual due to the edge effect is less than 0.15 mJy, 150" from the map edge. Chapter 4. Central Update Algorithm 78 SHADES measurement strategy. Each panel has 721 columns and 721 rows, which is 2163" wide. These results are encouraging. Provided noise is not greatly affected by the edge-effect, the edge-effect with a statistically plausible fake sky with the full SHADES measurement strategy is negligible compared to the expected statistical residuals. The simulation is repeated with the previous fake sky map except the null map centre is replaced with the centre of the statistically plausible fake sky map. The purpose of this simulation is to demonstrate the relatively fast disappearance of high spatial frequency structure with iteration, and the dominance of the edge effect residual in the converged residual map. (Figures 4.6 and 4.7) Plots of the residual map pixel standard deviation as a function of iteration are shown in Figure 4.6. Convergence to within 0.005 mJy of the converged standard deviation occurs by iteration 50. The converged standard deviation of 0.068 mJy is 0.002 mJy from the the converged standard deviation of the residual pixels in the edge-only simulation (Figure 4.3). The minimum of 0.0655 mJy in the standard deviation of the residual pixels occurs near iteration 40, and is consistent with the rising edge of the curve in Figure 4.3 at iterations less than 25, though the slope of the rise in Figure 4.3 must be less steep. This is purely heuristic, as the presence of sources in the fake map centre modifies the behaviour of the residual map mean with iteration. 2 The standard deviation in Figure 4.3 at iteration 30 is 0.18 mJy. The standard deviation in Figure 4.6 begins high, and then falls as the sky estimate map correctly estimates the sources outside of the map edge. With more iterations, the edge-effect residual propagates into the map and the standard deviation rises from its minimum to the converged value. Residual maps for arbitrary iterations between iteration one and 301 are depicted in Figure 4.7. 300" from the estimate map edges, the sky estimate is correct to within approximately 0.05 mJy, if the map average is estimated at —0.6 mJy, at iteration 90. The residual increases to approximately 0.1 mJy, 150" from the map edge and further increases with proximity to the map edge. Expected SHADES noise is around eight millijanskies per pixel and two millijanskies per beam. Thus, 2 In general, standard deviations of two different random variables do not add. Chapter 4. Central Update Algorithm 79 Residual Map Standard Deviation 0.09 0.08 0.07 0.06 400 600 Iteration Number 1000 Figure 4.6: Residual map standard deviation as a function of iteration during a sim-ulation of data reduction using the central update algorithm. The data consist of measurements of a statistically plausible fake sky obtained with the full SHADES measurement strategy. Convergence to within 0.005 mJy of the converged standard deviation occurs by iteration 50. The converged standard deviation of 0.068 mJy is 0.002 mJy from the the converged standard deviation of the residual pixels in the edge-only simulation (Figure 4.3). The minimum of 0.0655 mJy in the standard deviation of the residual pixels occurs near iteration 40. Chapter 4. Central Update Algorithm 80 ajgitfi » : • • • 1 -o.fl -o.6 -o.a Residual 40 Residual 50 mmm I H U N 200 400 600 Residual 80 Residual 90 _1 I I— 200 400 600 200 400 -0.7 -0.65 -0.6 -0.55 -0.5 Figure 4.7: Residual maps for a full SHADES simulation of the central update al-gorithm. The fake sky map is populated with sources drawn from a realistic source count model. Between iterations 50 and 90 the variance of the residual pixels reduces in the centre of the map and slightly in-creases toward the map edges. The pixel values are in mJy. Notice the high spatial frequency residuals, visible at iteration one, are largely gone by iteration 40. Chapter 4. Central Update Algorithm 81 in the SHADES survey, a 20% error on the measured flux emanating from a 10 mJy source is expected on 68% of the detected sources. At iteration 90 the algorithm has not converged and the dominant source of residual, over a 0.5 square degree area, is statistical. 4.2 Gaussian Noise Simulation The effect of additive Gaussian noise is determined by simulation. Random, white Gaussian noise is added to each measurement in a full SHADES simulation. The standard deviation of the noise, a, is chosen such that a/VN is around 8 mJy, where N is the number of times that a pixel is centred upon in a measurement. 8mJy is the pixel standard deviation due to noise expected from the SHADES measurement strategy. In this simulation a fake sky map with no sources is measured. Figure 4.8 shows a plot of the residual map standard deviation as a function of iteration. Pixels used in the computation of the residual map pixel standard deviation are located further than 150" from the map edge. The converged residual map pixel standard deviation is 9.95 mJy, 1.95 mJy higher than the expected value. The computation is repeated on the 1000th iteration residual map, using only pixels further than 600" from the map edge. The resulting residual pixel standard deviation is 9.06 mJy. The reduction in standard deviation is due either to large scale structure in the map or due to noisier pixels at the map periphery. Figure 4.9 shows the first and 101st iteration maps and their power spectrum maps. The power spectrum maps are invariant to 180° rotation. The lowest spatial frequency is centred in the power spectrum maps, the spectral resolution is 0.617 x 10~3 (arcsec-1) per pixel and the Nyquist frequency is | of a cycle per arcsecond. The 101st power spectrum map contains a bright cross centred in the map, as well as a fainter, checkered pattern, corresponding to the chop directions. Along the right half of the middle row there are a number of bright peaks. The first four of these bright peaks, moving right from the centre of the 101st power spectrum map, correspond to spatial wavelengths at 34.5", 22.8", 15.9", and 14.0", respectively. The first three of these wavelengths are roughly one-half of the chop throw sizes. The Chapter 4. Central Update Algorithm 82 10.0 Residual Map Standard Deviation 400 600 800 Iteration Number Figure 4.8: Residual map average and standard deviation as a function of iteration during a simulation of data reduction using the central update algorithm. Gaussian random noise is added to simulated measurements of a fake sky containing no flux. At iteration 100, the residual pixel standard deviation is 1.95 mJy higher than the expected (7measurement/v^/V The central-update algorithm converges when reducing measurements with Gaussian additive noise. Chapter 4. Central Update Algorithm 83 200 200 400 •i 0 0 1 L 1 1 1 - 0 . 0 2 - 0 . 0 1 0 0 . 0 1 0 . 0 2 PS I t e r a t i o n 1 PS I t e r a t i o n 101 • m • • • i i 1 u • o i l m r»* • ji$jk j 4' * si • • 11 ! HE • • m a • 100 200 300 400 500 100 200 300 400 500 I I "1 1 1 l e - 0 9 2 e - 0 9 3 e - 0 9 4 e - 0 9 5 e - 0 9 Figure 4.9: Measurements with Gaussian noise are iteratively reduced using the cen-tral update algorithm. The maps are, clockwise from upper left, the 1 s t iteration map, the 101st iteration map, the power spectrum of the 1 s t it-eration map and the power spectrum of the 101st and 1 s t iteration map. Note the structure that appears in the 101st iteration maps. Chapter 4. Central Update Algorithm 84 timestream noise covariance matrix is diagonal in this simulation, and the pixel-pixel correlations are due entirely to the measurement strategy. The power spectrum is the Fourier transform of the correlation function; structure in the power spectrum map is due to pixel-pixel correlations. Figure 4.4, by juxtaposing the first and 101st iteration power spectrum maps, illustrates the growth of pixel-pixel correlations formed during the convergence of the iterative deconvolution algorithm. Examining the sky estimate maps we see a cross-hatch pattern and large scale fluctuations present in the 101st iteration map that are not present in the first it-eration sky estimate map. As mentioned in the previous paragraph, this structure appears in the power spectrum maps and is correlated pixel-pixel noise. An intuitive understanding of this structure is possible. The SCUBA double-differencing mea-surement strategy, for scales much larger than the chop distance, is approximately a double spatial derivative of the sky flux. A second derivative is a high-pass filter. Thus, deconvolving the chops requires amplifying the low spatial frequencies to undo the filter, resulting in large spatial scale noise amplification. The amplification of the large scale noise will increase the residual map pixel standard deviation. The simulation is repeated with a different fake sky, this time containing four, 30mJy sources and a single 50mJy source. Figure 4.10 shows the portion of the 1st and 101st sky estimate maps, and their corresponding power spectrum maps, which contain the large sources. Again, the 101st iteration map exhibits cross-hatching and large scale fluctuations. The regions surrounding the sources exhibit no correlations not exhibited by regions far from the sources in the same map. Though not quantified, the pixel-pixel noise correlations in the iterated map are not visibly affected by the presence of a 50 mJy source. The simulation is repeated a final time, this time using the statistical model of SCUBA sources to create the fake sky. The sky estimate maps and power spectrum maps for iteration one and 101 are shown in Figure 4.11. The iteration one power spectrum map contains new structure around the centre of the map. This structure is at wavelengths greater than 34.5" and is due to the presence of many sources in the map. The 101st power spectrum map appears very similar to the previous 101st power spectrum maps, the cross-hatching, and the bright peaks are all present, in Chapter 4. Central Update Algorithm 85 I t e r a t i o n 1 I t e r a t i o n 101 300 350 400 450 300 350 400 450 Figure 4.10: No source streaking, other than the existing cross-hatch pattern, is ap-parent, in a simulated reduction of measurements of four large sources containing Gaussian additive noise using the central update algorithm. The pixel values are in mJy. Chapter 4. Central Update Algorithm 86 Iteration 1 Iteration 101 200 400 600 200 400 600 K i l l -0.02 -0.01 0 0.01 0.02 Power Spectrum 1 Power Spectrum 101 0 l e - 0 9 2 e - 0 9 3e -09 4 e - 0 9 5 e - 0 9 Figure 4.11: Measurements of a statistically plausible fake sky with additive Gaus-sian noise are iteratively reduced using the central update algorithm. The maps are, clockwise from upper left, the 1st iteration map, the 101st iteration map, the power spectrum of the 1st iteration map and the power spectrum of the 101st iteration map. The presence of sources appears as low frequency power in the 1st iteration power spectrum map. Chapter 4. Central Update Algorithm 87 the same locations, at the same values. 4.3 More Realistic Noise Simulation To test the central update algorithm in a more realistic application, a simulation of the algorithm is performed with more realistic noise. Realistic noise is created by extending SCUBA data taken during three observations of a blank field to the 4.2 x 107 measurements required in the SHADES simulation. The extension is performed as follows: the existing blank field measurements are divided into equally sized groups of contiguous measurement samples. The original blank field measurements are then extended by adding a randomly selected contiguous piece of data, drawn from the group of contiguous, existing data subsets. This process is continued, until a sufficient number of data samples is obtained. Figure 4.12 shows example timestream noise samples for four arbitrarily chosen bolometers and the magnitude of their Fourier transforms. Each bolometer timestream is pure noise and has 21,120 samples, requiring an extension factor of 22. The extension factor is arbitrarily chosen. The time streams and Fourier transform magnitudes have maximum values consistent with the timestreams and Fourier transform magnitudes of real, pre-processed SHADES data plotted in Figure 2.5. The randomly extended timestreams lack the unstationarity3 of the actual SHADES data and there are spikes in the Fourier transform magnitudes not present in the real SHADES data. These spikes are likely due to the fixed size of the timestream subsets used in the random extension of the timestreams. The randomly extended noise is added to simulations of the fake sky map with a statistically plausible fake sky in a full SHADES simulation of the central update algorithm. Plots of the residual map standard deviation are shown in Figure 4.13. At iteration 100, the standard deviation of the residual map pixels is approximately 0.35 mJy from the converged value of 11.55 mJy. Repeating the calculation for the 1000th residual map omitting the outer 600" results in a residual pixel standard deviation of 8.89 mJy. This value is much nearer the expected value of 8mJy. The 3The variance of each bolometer timestream in Figure 2.5 changes with time. That is, the data samples are drawn from time-dependent probability density functions. Chapter 4. Central Update Algorithm 88 Fourier Transform MoqnitudJ Fourier Tranitorm Mognituflc r™j«nej (tu) Figure 4.12: First four plots are of the randomly extended noise sequences for bolometers 1, 2, 9, and 35. The bolometer selection is arbitrary. Last four plots are the corresponding Fourier transform magnitudes. Note the lack of j noise in the spectrums. There is a bump at approx. 0.007 Hz in the bolometer 35 spectrum, possibly the Fourier transform spike described in chapter 2. Chapter 4. Central Update Algorithm 89 Residual Map Standard Deviation 200 400 600 Iteration Number 800 1000 Figure 4.13: Residual map standard deviation as a function of iteration in a central update algorithm simulation with a statistically plausible fake sky map and the full SHADES measurement strategy. Chapter 4. Central Update Algorithm 90 larger decrease of the residual pixel standard deviation between this simulation and the simulated reduction of Gaussian distributed measurements of a blank sky is due either to greater contributions to the standard deviation from the map edges or due to a greater contribution to the pixel variance from large scale structure. By looking at the 101st iteration sky estimate map in Figure 4.14 and comparing to the 101st iteration sky estimate map in Figure 4.9, we see that there is greater large scale fluctuation in Figure 4.14, especially at the edges. This is consistent with the greater reduction in the residual pixel standard deviation with decreasing edge contributions to the calculation. Figure 4.14 depicts the first and 101st iteration maps with their corresponding power spectrum maps. The maps in Figure 4.14 closely resemble those produced in the simulations of the central update algorithm with measurements containing Gaussian additive noise. In this simulation, it is unlikely that the timestream noise covariance matrix is diagonal, and one would expect to see extra pixel-pixel corre-lations. From Figure 4.15 this is not evident, and if present, has a relatively small effect. From iteration 100 to convergence, there is a 0.4 mJy change in the stan-dard deviation of the residual map pixels (See Figure 4.13). Though this change is small, any outlier pixels are greatly suppressed in the computation due to the normalization by a large number of pixels. It is known that the algorithm will converge, neglecting edge-effects, to an optimal solution. In noiseless simulations, by iteration 90 residuals span hundreds of arcseconds and are below 0.1 mJy (See Figure 4.7). One wishes to cease iteration when improvements due to further con-vergence are negligible. Iteration amplifies the large scale noise and correlates the pixels, while removing the negative side-lobes present around sources in the first iteration sky estimate map. In principal, the converged map is optimal; however, subsequent optimal use of the converged map will require accounting for the pixel-pixel correlations. The degree to which continued iterations between iteration 200 and iteration 1000 affect the science product will only be determined by simulation involving the method of accounting for the pixel-pixel correlations, and will depend on the objective of the analysis. The pixel-pixel correlations are complicated and it is not intuitive how their presence or enhancement will benefit, for instance, source Chapter 4. Central Update Algorithm 91 Iteration 1 Iteration 101 -0.02 -o.oi o 0.01 0.02 Power Spectrum 1 Power Spectrum 101 100 200 300 400 500 100 200 300 400 500 Figure 4.14: Clockwise from upper left: Sky estimate maps for iteration one, 101, power spectrum map for iteration 101, power spectrum map for it-eration one. Note the similarity between these maps and the maps produced in simulations with Gaussian additive noise. (Figure 4.12) The off-diagonal terms of the timestream noise-covariance matrix play a minor role in the pixel-pixel correlations present in the 101st itera-tion map. The pixel-pixel correlations are predominantly due to the off-diagonal terms in the PTP matrix for the SHADES measurement strategy. The PTP matrix was introduced in chapter two. Chapter 4. Central Update Algorithm 92 detection. The next section demonstrates the effects of measurement strategy on pixel-pixel correlations by modifying the SHADES chopping scheme in simulation. 4.4 Chopping Effects On Pixel-Pixel Correlation In this section the effect on the pixel-pixel correlations is investigated through simu-lation by increasing the number of chops to 12. Specifically, one half of the measure-ments are performed along chops rotated 45 degrees from the orignal chop directions, adding six new chops to the measurement strategy. Figure 4.15 shows the first, and 101st sky estimate maps and their corresponding power spectrum maps. The total number of measurements and the additive noise are the same as in the randomly extended noise simulations described in section 4.3. The fake sky is statistically plausible. The 101st sky estimate map in Figure 4.15 exhibits less cross-hatching than the 101st sky estimate map in Figure 4.11. The large scale structure, with a half-wavelength aroung 600", seems unaffected. Examining the corresponding 101st power spectrum maps, we see a reduction in power along the very bright cross cen-tred in Figure 4.11, and a change of the cross-hatch pattern to a star pattern. For both power spectra, the pattern corresponds to the direction of the chops. The power is distributed evenly in all directions in an inner square 100 pixels wide centred in the 101st power spectrum map in Figure 4.15. This distribution of correlations along new directions is consistent with the arguement presented in Section 6.3 regarding the reduction of pixel-pixel correlations by involving every pixel in measurements involving many different pixels. 4.5 Chapter Results The central-update algorithm, as applied to SHADES, is evaluated has been eval-uated by simulation. The residual due to the edge effect is found to be negligible compared to the statistical noise in the SHADES survey for most of the central half-square degree of the sky estimate map (Section 4.1). The presence of strong sources in the measured fake sky map, does not correlate the noise. (Section 4.2) For Chapter 4. Central Update Algorithm 93 a statistically plausible fake sky with randomly-extended realistic noise, the stan-dard deviation of the residuals upon convergence is 11.56 mJy, 0.4 mJy higher than the standard deviation of the residual map pixels at iteration 101 (Section 4.3). In Section 4.4 the effect of changing the SHADES chopping pattern is investigated and the pixel-pixel correlations found to be significantly reduced. Besides further simulations, to use the iteratively deconvolved sky estimate map, either the pixel-pixel noise covariance matrix must be used with subsequent source detection, or the pixel correlations must be reduced. In the chapter on Discussion, the effect of follow-up observations aimed at reducing the pixel-pixel correlation is addressed. As well, the use of priors to constrain the correlation is briefly mentioned in the section on Noise & Correlation Removal. Chapter 4. Central Update Algorithm 94 I t e r a t i o n 1 Iteration 101 J I I I 1 L -0.02 -0.01 0 0.01 0.02 (Jansky) 100 200 300 400 500 100 200 300 400 500 Figure 4.15: Simulation of data reduced using the central update algorithm. The data are simulated measurements of a statistically plausible fake sky with blank field extended noise samples. One half of the chops are rotated 45 degrees from the east-west and north-south chop directions. Much of the cross-hatching visibile in the 101st iteration map in Fig-ure 4.14 has disappeared. The strong cross-shape in the 101st power spectrum map in Figure 4.14 has been reduced. Large scale structure with half-wavelengths around 600" seem unaffected. Chapter 5. SHADES Subaru XMM-Deep Field 95 Chapter 5 SHADES Subaru XMM-Deep Field A list of sources detected in a measurement map of the Subaru XMM-Deep Field is presented. The reduced data includes all SHADES observations up to the 1st of Febuarary, 2004. 5.1 Data Reduction The data is pre-processed using the techniques described in Section 2.4, with the exception that the flux conversion factors are set to a nominal value of 225Jy V - 1 Figure 5.1 shows the measurement map of the Subaru XMM-Deep Field. Each pixel has a standard deviation of approximately 8 mJy, and is fully measured for most of the map. Figure 5.2 shows the hit map for the Subaru XMM-Deep Field. The value of each pixel in the hit map is the number of times each pixel is centred upon in a measurement and the noise maps are computed assuming a diagonal timestream noise covariance matrix. That is, the pith pixel in the noise map is where Nm is the number of measurements centred upon the ith pixel and on is the standard deviation of the set of measurements containing the n t h measurement. The set of measurements is grouped by bolometer and observation file, and nominally contains 320 measurements. (5.1) Chapter 5. SHADES Subaru XMM-Deep Field 96 Subaru XMM-Deep F i e l d Measurement Map 500 400 1 a» ii 300 200 • - ' 200 300 400 500 -0 .03 -0 .02 - 0 . 0 1 I 0 (Jansky) I 0.01 0.02 0.03 Figure 5.1: The Subaru XMM-Deep Field measurement map. Units are in Jy, pixels are 3" square. The horizontal and vertical axes are specified in units of pixels. Chapter 5. SHADES Subaru XMM-Deep Field 97 Subaru XMM Deep F i e l d : Noise Map ••fWlaiiiSiJi J L L fa - h J i r n 1 1 r 200 300 400 500 0.006 0.008 0.01 0.012 0.014 Subaru XMM Deep Fie l d : Hit Map r k * mm iwiJI 200 300 400 500 0 50 100 150 Figure 5.2: Subaru XMM-Deep Field noise map and hit map, top and bottom, re-spectively. The noise map units are in Jy and the noise per pixel is between 7mJy and 11 mJy for the portions of the map fully measured by the SHADES strategy measurement strategy. The SHADES tripos pointing pattern is clearly visible. Chapter 5. SHADES Subaru XMM-Deep Field 98 5.2 Source Detection Source detection can be posed as the task of determining the model which will provide the best fit to the pixel values. Ideally one would determine the most likely model by computing the probabilities of different models by marginalizing over all model parameters. Proper computation of these probabilities will intrinisically incorporate an Occam's penalty, and will select the number of sources warranted by the data.[? , BayesTutorialjhe candidate models are no sources present, one source present, two sources present, etc. Because the practical solution to this problem is difficult, a more simple approach is taken.[5] [19] [18] 1 The method is as follows: for every pixel in the map, assume that there is a source centred upon the pixel under consideration, assume that the prior probabilities are uniformly distributed , the pixel noise is Gaussian and the source full-width half-maximum is given by the beam of the telescope. Compute the most probable amplitude for the source by finding the source amplitude which maximizes the likelihood. That is, given that there is a source centred upon the pixel under consideration, the likelihood is P(d| S)ocexp(-^) (5.2) where P (d\s) is the probability of the data, d, given the model, s, and Xi = NfhzMt (5.3) i 1 where j is the j t h pixel under consideration (centred-upon), pi, is the ?th pixel, /3 is the source amplitude, fij, is the value of the source model centred upon pixel j, for the ith pixel, and o is related to the telescope beam and determines the beam full-width half-maximum. The most likely amplitude, J3, is obtained by finding the maximum of x2 and, for the jth pixel, is given by A = ^ ^fr (5-4) 1 T h e c o m p u t a t i o n requ i res for every m o d e l cons ide red a n in teg ra l of t he l i k e l i h o o d over a l l source centres a n d a l l source a m p l i t u d e s a l l owed by the pr io rs . Chapter 5. SHADES Subaru XMM-Deep Field 99 in practice, will be negligible for i's significantly different than j, and the sums involving i over all of the map pixels can be limited to the contributing terms. Once the most probable amplitude is computed, the next thing to do is to com-pute the uncertainty in the amplitdue estimate. This is accomplished by taylor expanding the x2 about the most probable amplitude, and taking the second, and final, term to be the uncertainty. Because the taylor expansion is about the maxi-mum, the first order term is zero. The second order term gives an estimate of the width of the Gaussian likelihood. The uncertainty, is given by -0.5 (5.5) The source model, fij, used in the above computation is the measurement map formed from measuring a fake sky map containing a unitless source with amplitude one with the SHADES strategy. Figure 5.3 is a map of the mask. The next step in the algorithm is to compute the amplitude uncertainty using equation 5.4, and compute the ratio between the most probable amplitude and the amplitude uncertainty. This process is repeated for every pixel in the map. The most-probable amplitude map, uncertainty map, and the ratio map for the Subaru XMM-Deep Field are shown in Figures 5.4. Once the ratio map is created, the maximum pixel in the ratio map is found. A source, with an amplitude obtained from the corresponding location in the most probable amplitude map, is subtracted from the measurement map. The subtracted source is considered detected at the pixel value in the ratio map, and is recorded. The most probable amplitude, uncertainty, and ratio maps are re-computed, and the process repeated until there are no pixels greater than three sigma left in the ratio map. The list of recorded sources will be all sources greater than three standard deviations from the map mean of zero. These sources are considered detected at greater than three sigma. To check the source detection a histogram of the pixels in the uncleaned most probable amplitude map is compared with a histogram of the pixels in the source cleaned most probable amplitude map. The histograms are plotted in Figure 5.5, along with a scaled probability density function corresponding to the probability of the union of the pixel probabilities. The scaled probability density function Chapter 5. SHADES Subaru XMM-Deep Field 100 Figure 5.3: Mask used in source detection to fit the amplitude of a point centred on each pixel. Note the negative sidelobes. The central Gaussian has a full-width half-maximum of 14.7". Each pixel is 3" square. The mask is 150" square. The mask mean pixel intensity is zero. Chapter 5. SHADES Subaru XMM-Deep Field 101 Most P r o b a b l e A m p l i t u d e Map A m p l i t u d e U n c e r t a i n t y Map >.002 0.003 n Figure 5.4: Subaru X M M Deep Field: Clockwise from upper-left, most-probable amplitude map, amplitude uncertainty map, and ratio map. The source detection acts as a statistically weighted low-pass filter. The funny pat-tern along the edges of both the most probable amplitude map and the ratio map is due to the averaging of only a few noisy values in the as-signment of the pixel values near the edges. The saturated, mask sized regions in the upper right and lower left of the most probable amplitude map and the ratio map are due to very large valued pixels along the edge of the Subaru measurement map. Chapter 5. SHADES Subaru XMM-Deep Field 102 represents the expected histogram made by histogramming the pixels of a pure noise most probable amplitude map. The probability density function is constructed by summing the probability density functions for each of the pixels. Each pixel probability density function is given by a Gaussian whose width is specified by the corresponding pixel value in the amplitude uncertainty map. The sum of the varying width Gaussians gives the probability density function for the union of the most probable amplitude map pixels. The union probability density function is scaled by eye. Near the edge of the map the noise properties of the map have not been properly investigated. The source detection mask (Figure 5.3) is a measurement map of a single source and contains negative sidelobes, the largest of which is centred 68" from the source centre. Because of this, the mask is large and much data near the edges of the map would be lost if the convolution were performed by truncating the edge of the most probable amplitude map by half the mask size. Instead, the convolution is performed as an average, such that the mask does not have to lie completely within the map prior to getting a non-zero value. When the mask is only partly within the map, the computation is more susceptible to a few, very noisy pixels. This is clear in the edges of the most probable amplitude map depicted in Figure 5.5. If histograms of the noise were produced of the edge pixels and compared with a noise model constructed from the uncertainty map, then detections at map edge would be included. This test has not been done and the author chooses to remove sources with a chop substantially off the edge of the map from the list of detected sources. To check this removal, a map of the Subaru XMM-Deep Field with a cross drawn over the source positions is created. Figure 5.6 depicts the map. The lists of detected sources are given in Tables 5.1, 5.2 and 5.3, in three parts, for the Subaru XMM-Deep Field. The columns and rows of the sources corresponding to the locations in the maps shown in this chapter are included in the tables. The correspondence of the Subaru XMM-Deep Field source list with the list given by the SHADES community as of April 12, 2004 is demonstrated in Figure 5.7. Considerable non-correspondence exists between the sources detected by the two groups. The expected non-correspondence is small as the data prior to reduction Chapter 5. SHADES Subaru XMM-Deep Field 103 1OOOO F -0 .02 -0.01 0.00 0.01 0.02 Pixel Value (mjy) Figure 5.5: Subaru XMM Deep Field most probable amplitude pixel histograms with noise model. Top, histograms produced from pixels with noise less than 4mJy. Bottom, histograms made from pixels with noise less than 3mJy. There is an excess of positive, and negative flux in the uncleaned most probable amplitude map, due to the presence of sources. The apparent lack of false detections in the top panel is due to the greater pixel noise populating the top histogram's positive tail. Chapter 5. SHADES Subaru XMM-Deep Field 104 s r c O n R a t i o B l a c k . f i t s _ 0 _ _ l I ( S N R ) Figure 5.6: Subaru XMM Deep Field source detections. The black crosses mark the location of detected sources. The closest detections are around 70" from the map edge. Chapter 5. SHADES Subaru XMM-Deep Field 105 Srcs Detected With SHADES Data Up To Feb 1, 2004 5 1 - 5 • 1 1 , | , , 1 + , 1 1 1 1 UBC 0 0 Rest of SHADES . + - -o 0 - o + + -o O i * + * + o + 0 " o + + o < _ + o - «\ Hi o + o o o o < o 0 + o o * + + o o + + + -• + + O o -T •r 0 o 1 1 + , 1 , 1 + 1 1 o 1 , , , , - 1 0 - 5 0 5 10 RA Offset (arcmins) Figure 5.7: Subaru XMM-Deep Field source list correspondence. There are 26 un-paired UBC detections, 25 unpaired detections from the SHADES com-munity, and 35 matches. Chapter 5. SHADES Subaru XMM-Deep Field 106 ^ 4.4 in 4.2 4.0 0 100 200 300 400 Iteration Figure 5.8: Sky estimate map pixel standard deviation as a function of iteration. The standard deviation is within 0.01 mJy of its converged value at iteration 100. are identical. That is, in principle, the source lists should only differ by differences in the reduction and source detection algorithms. The scatter may be due to improper calibration, or to an unknown observational problem. Different source detection algorithms are used between the groups, and mild disagreement on weakly detected sources is plausible. 5.3 Deconvolved Map Source Detection As a means of quantifying the usefulness of the deconvolved map, the central update algorithm is used to reduce the Subaru XMM-Deep Field data, and subsequent, non-optimal source detection is performed. Plots of the standard deviation of the map pixels as a function of iteration are shown in Figure 5.12. The deconvolved map pixel standard deviation is within 0.01 mJy of its converged value at iteration 100. Figure 5.13 is the 401st iteration sky estimate map. The large scale structure and cross-hatching seen in simulation in chapter four is evident. Measurements centred upon pixels involved in the centre of measurements less than thirty times Chapter 5. SHADES Subaru XMM-Deep Field 107 Figure 5.9: 401st sky estimate map. The pixel-pixel correlations are apparent. Mea-surements centred upon pixels involved in the centre of measurements less than thirty times are discarded. Source detection on this map will have to deal with the large scale structure, and if it is to be optimal, the cross-hatching. Chapter 5. SHADES Subaru XMM-Deep Field 108 are discarded. The poorly measured edges in this map affect the stability of the centre update algorithm and increase the pixel-pixel correlations. Discarding poorly measured pixels from the data reduction, which means discarding the measurements of the poorly measured pixels, increases the stability of the algorithm and decreases the pixel-pixel correlations present in the resulting map. Source detection is performed on the 401st iteration map. The choice is arbitrary, any map after iteration 100 would suffice. To remove some of the bias due to the large scale structure, the source detection is performed with the same mask used in the source detection performed on the 1st iteration map. This mask has a zero mean and is evenly symmetrical. Odd functions and constant offsets will not be measured; however, fluctuations possessing even symmetry or no symmetry on the scale of the mask size will result in non-zero response and affect source detection. The source detection is performed using the algorithm described in section 5.2. Figure 5.14 shows the correspondence between the sources detected on the decon-volved map and the sources detected in section 5.2 for the Subaru XMM-Deep Field as well as the correspondence with the source list provided by the rest of the SHADES community. The source list resulting from the source detection on the de-convolved map has not been pruned to reject detections at the edge of the map. The correspondence between any of the two source lists is equal. There is no indication which list is more reliable. The optimal source detection, taking into account the full pixel-pixel correlation matrix, is expected to out-perform the source detection applied in section 5.2 to the measurement map. The correlation matrix will need to be large enough to account for the dominant correlations, and will need to be inverted. It is possible that the matrix will take a form amenable to fast inversion; however, the scrambling of the time series will likely result in a pixel- pixel correlation matrix lacking nice symmetries. Nonetheless, the investigation has not been done and the author is speculating. Chapter 5. SHADES Subaru XMM-Deep Field 109 o U J Q 10 + o * A A & A + -F A + A -5h A - 1 0 1 o A O A A O A + A -k> - | i i i r + Iteration 1 O Iteration 4 0 ' A Rest of Shadejs A A & o A _ i i i i _ 10 RA (arcsecs) Figure 5.10: Deconvoluted source detection correspondence with the sources de-tected on the measurement map (1st iteration map) and with the source list provided by the rest of the SHADES community using the data ac-quired up to February 1st, 2004. There is no indication which source detection method is more reliable. Chapter 5. SHADES Subaru XMM-Deep Field 110 Subaru X M M Deep Field Source Detections - Part I RA (°) DEC (°) Significance Amplitude (mJy) Uncertainty (mJy) Col Row 34.3779 -4.9934 7.04 14.2 2.0 506 368 34.5151 -4.9242 5.66 11.9 2.1 342 451 34.4791 -4.8834 4.65 10.5 2.3 385 500 34.4130 -5.0892 4.54 9.1 2.0 464 253 34.4356 -4.9316 4.52 10.1 2.3 437 442 34.4255 -4.9417 4.56 9.6 2.1 449 430 34.5653 -4.9008 4.45 9.1 2.1 282 479 34.4983 -5.0842 4.16 8.5 2.0 362 259 34.4682 -5.0825 4.12 9.0 2.2 398 261 34.4289 -5.0742 4.08 9.1 2.2 445 271 34.5009 -5.1283 3.85 9.1 2.4 359 206 34.3762 -4.9975 3.76 6.9 1.8 508 363 34.4029 -5.0766 3.75 6.9 1.9 476 268 34.5787 -5.0467 3.69 9.2 2.5 266 304 34.4832 -5.1067 3.69 8.1 2.2 380 232 34.3762 -5.0167 3.67 6.7 1.8 508 340 34.4289 -5.0967 3.67 7.5 2.0 445 244 34.3695 -4.9850 3.61 7.4 2.1 516 378 34.4339 -5.0375 3.61 7.8 2.2 439 315 34.4113 -5.0608 3.60 6.8 1.9 466 287 34.3528 -4.9791 3.56 7.7 2.2 536 385 34.6054 -4.9541 3.52 9.2 2.6 234 415 34.5217 -4.9875 3.49 10.1 2.9 334 375 Table 5.1: Subaru X M M Deep Field Source Detection List - Part I Chapter 5. SHADES Subaru XMM-Deep Field 111 Subaru X M M Deep Field Source Detections - Part II RA (°) DEC (°) Significance Amplitude (mJy) Uncertainty (mJy) Col Row 34.5310 -4.9875 3.82 11.5 3.0 323 375 34.5636 -5.0308 3.45 7.8 2.3 284 323 34.4866 -4.8775 3.43 9.5 2.8 376 507 34.4022 -4.9325 3.42 7.7 2.3 477 441 34.4063 -5.0658 3.39 6.6 2.0 472 281 34.5326 -5.0675 3.39 7.6 2.3 321 279 34.3878 -5.0333 3.38 6.2 1.8 494 320 34.4189 -5.0217 3.36 7.7 2.3 457 334 34.6346 -4.9951 3.34 8.6 2.6 199 366 34.5267 -4.9258 3.33 6.7 2.0 328 449 34.5577 -4.9616 3.29 6.3 1.9 291 406 34.4415 -4.9650 3.28 6.5 2.0 430 402 34.4648 -5.0467 3.27 6.6 2.0 402 304 34.5686 -4.9200 3.25 7.2 2.2 278 456 34.3870 -4.9708 3.25 6.6 2.0 495 395 34.3879 -4.9525 3.25 7.0 2.2 494 417 34.4615 -5.1092 3.22 6.9 2.1 406 229 34.5151 -5.0708 3.22 7.3 2.3 342 275 34.4950 -4.9083 3.21 6.3 2.0 366 470 34.5041 -4.8867 3.19 7.2 2.3 355 496 34.3553 -4.9942 3.18 6.6 2.1 533 367 34.4515 -5.0008 3.15 8.1 2.6 418 359 34.5469 -5.0467 3.15 6.6 2.1 304 304 Table 5.2: Subaru X M M Deep Field Source Detection List - Part II Chapter 5. SHADES Subaru XMM-Deep Field 112 Subaru X M M Deep Field Source Detections - Part III RA (°) DEC (°) Significance Amplitude (mJy) Uncertainty (mJy) Col Row 34.3737 -4.9875 3.13 6.5 2.1 511 375 34.3921 -4.9825 3.15 5.6 1.8 489 381 34.5652 -4.9316 3.12 6.4 2.1 282 442 34.5845 -4.9458 3.11 6.9 2.2 259 425 34.5619 -4.9891 3.10 6.4 2.1 286 373 34.4239 -4.9758 3.10 6.2 2.0 451 389 34.4106 -4.8783 3.10 7.9 2.6 467 506 34.5452 -5.1292 3.09 10.5 3.4 306 205 34.4607 -4.9275 3.09 6.0 2.0 407 447 34.3712 -5.0458 3.09 7.3 2.4 514 305 34.4255 -4.9033 3.03 6.3 2.1 449 476 34.4214 -5.1341 3.03 9.1 3.0 454 199 34.4539 -5.0866 3.01 7.1 2.3 415 256 34.3946 -5.0775 3.00 6.7 2.2 486 267 34.4456 -5.0750 3.00 6.6 2.2 425 270 34.4673 -5.0783 3.04 7.0 2.3 399 266 Table 5.3: Subaru X M M Deep Field Source Detection List - Part III Chapter 6. Discussion 113 Chapter 6 Discussion This chapter summarizes the main results of the thesis, discusses follow-up ob-servation aimed at reducing the pixel-pixel correlations in the deconvolved map, discusses general issues when choosing a measurement strategy, poses the problem of determining an optimal measurement strategy, and discusses future theoretical work regarding the reduction of pixel-pixel correlations in the deconvolved map. 6.1 Main Results Data calibration using a flux conversion factor computed from calibration obser-vations of CRL618, Uranus and Mars is fournd to have a 10% difference in the correspondence of detected sources compared to the sources detected with the same data except using a flux conversion factor equalling 225 JyV - 1. The chop sizes of the differenct calibration observations is ignored (Section 2.4.1). The iterative deconvolution algorithm is generalized and formulated as the so-lution to a difference equation which gives the map for any iteration (Equation 3.35). The generalized deconvolution algorithm is shown to diverge in SHADES simulations. The convergence properties of the generalized iterative algorithm are determined by the eigenvalues of a characteristic matrix, with small dependence upon the measurements. The formulation is checked theoretically and by simula-tion. The iterative algorithm is theortically shown to converge to the least-squares optimal solution. (Section 3.59) A known work-around, the central-update algorithm, is explained by the theory (Section 3.5.3). The central-update algorithm, as applied to SHADES, is evaluated by simulation (Chapter 4). The residual due to the edge effect is found to be negligible compared to the statistical noise in the SHADES survey for most of the Chapter 6. Discussion 114 central half-square degree of the sky estimate map (Section 4.1). For a statistically plausible fake sky with randomly-extended realistic noise, the standard deviation of the residuals upon convergence is 11.56 mJy, 0.4 mJy higher than the standard deviation of the residual map pixels at iteration 101 (Section 4.3). In Section 4.4 the effect of changing the SHADES chopping pattern is investigated and the pixel-pixel correlations found to be significantly reduced. In chapter 5 a source list is presented for source detections in the Subaru XMM-Deep Field using the first iteration, or measurement map. The source list is com-pared with a list of sources provided by the SHADES community. Source corre-spondence is around 50%. 6.2 Follow-Up Observation The effect of follow-up observation on the iterated map pixel-pixel correlations is investigated by modifying the SHADES chopping strategy in simulation. Two sim-ulations are performed, the first with one extra chop, constituting an increase in the number of SHADES measurements by ^ t h the original number, and a second simulation with two mutually orthogonal chop directions of the same size, increasing the number of SHADES measurements by |th. Figure 6.1 shows the first and 101s* sky estimate maps and their corresponding power spectra, for a simulation with one extra chop per pointing. The chop is oriented at a 45 degree angle to the existing chop directions and has a size of 44". The extra chop per pointing constitutes an extra j j r t h of the SHADES observations. The fake sky is statistically plausible and the noise consist of randomly extended blank sky samples described in section 4.3. The power along the dominant cross in the 101st power spectrum map in Figure 4.14 is somewhat reduced, and the cross-hatch pattern modified. The maps are less noisy as there are more measurements. Despite this, the pixel-pixel correlations in the 101st power spectrum map are larger than in the 12 chop simulation resulting in the 101st power spectrum map shown in Figure 4.15. The effect of two extra chops is investigated. The simulation conditions are the Chapter 6. Discussion 115 Iteration 1 Iteration 101 -0.02 -0.01 0 0.01 0.02 (Jansky) Power Spectrum 1 Power Spectrum 101 100 200 300 400 500 100 200 300 400 500 Figure 6.1: Simulation of data reduced using the central update algorithm. The data are simulated measurements of a statistically plausible fake sky with blank field extended noise samples. An extra chop's worth of data has been added to measurements made with the full SHADES measure-ment strategy. The chop has a size of 44" and is rotated 45 degrees with respect to the original chop directions. Clockwise from upper left, first iteration sky estimate map, 101st iteration sky estimate map, 101st iteration power spectrum map, first iteration power specrum map. The power along the dominant cross in the 101st power spectrum map in Figure 4.14 is somewhat reduced, though not as much as in Figure 6.1, and the cross-hatch pattern modified. Chapter 6. Discussion 116 same as in the above simulation except two extra chops for each pointing are added instead of one. The extra chops are rotated 90 degrees with respect to each other and are rotated 45 degrees from the east-west, north-south chop directions. Figure 6.2 shows the sky estimate maps and their corresponding power spectra for the first and 101st iteration maps. Power along the prominant cross in the 101st power spectrum maps in Figures 6.1 and Figure 4.14 is further reduced. At low frequencies the power has been concentrated to narrow lines along the chop directions. The cross-hatch pattern throughout the 101st power spectrum map is changed and reduced. It is clear that follow-up observation with different chops reduces the pixel-pixel correlations present in the data. This effect should be two-fold; both a reduction in the pixel standard deviation by y/~N and the effect of the extra cross-linking described in the next section. However, the enhancement has a negligible effect on the large scale fluctuations in the deconvolved map (Compare the 101st sky estimate maps in Figures 4.14, 6.1, and 6.2) and modestly improves the pixel-pixel correlations, while the required observation time is long. 6.3 Measurement Strategy The SHADES measurement strategy is not optimal when producing deconvolved maps. To produce a diagonal pixel-pixel noise covariance matrix, both the timeline noise covariance matrix and the P1P matrix must be diagonal. Absolute power measurements produce a diagonal PTP matrix but often have a non-diagonal noise covariance matrix. Many sources of noise correlation in the timelines are reduced by performing differential measurements; however, the reduction is traded for a diago-nal PTP. Experiments with good measurement strategies, such as WMAP, achieve nearly diagonal pixel-pixel noise covariance matrices by taking differential measure-ments involving every pixel with as many other pixels in the map as possible, and by repeating the entire experiment many times. Each differential measurement corre-lates all of the pixels involved in the measurement, while increasing the appropriate diagonal term in PTP by, for instance, one in the central update algorithm scheme with noiseless measurements, and the off-diagonal terms by a small amount. By Chapter 6. Discussion 117 Iteration 1 Iteration 101 ~\ r -0 .02 -0 .01 0 0.01 0.02 (Jansky) Power Spectrum 1 Power Spectrum 101 ~1 1 1 1 l e - 0 9 2<e-09 3e-09 4e-09 Se-09 Figure 6.2: Simulation of data reduced using the central update algorithm. The data are simulated measurements of a statistically plausible fake sky with blank field extended noise samples. Two extra chop's worth of data has been added to measurements made with the full SHADES measurement strategy. The chops have a size of 44" are rotated 45 degrees with respect to the original chop directions and are mutually orthogonal. Clockwise from upper left, first iteration sky estimate map, 101st iteration sky estimate map, 101st iteration power spectrum map, first iteration power specrum map. Power along the prominant cross in Figures 6.1 and Figure 4.14 is further reduced. At low frequencies the power has been concentrated to narrow lines along the direction of the prominant cross. Chapter 6. Discussion 118 involving many different pixels in many measurements of the same pixel, the diag-onal terms grow with respect to the average off-diagonal term. Correlating every pixel with every other pixel, a map relatively free of pixel-pixel correlations can be constructed. An additional consideration, is that experiment repetition allows for the dis-crimination between unchanging astrophysical sources and j noise. The rate of experimental repitition will determine the frequencies at which this separation is possible, ie. The less the slow correlated noise changes between experimental re-peats, the more of it that can be determined and removed. 6.4 Optimal Measurement Strategy The problem of determining the optimal measurement strategy can be stated as the task of determining the linear combination of absolute power measurements which will result in a maximally diagonal pixel-pixel covariance matrix. Consider the noise covariance matrix formed from the noise timestream of an absolute power experiment, N-1 = (n-nT) (6.1) Form differential data from the absolute power measurements by convolving with a measurement function, h: TV" 1 = (^(h®nj (h®njT^ (6.2) The new pixel-pixel noise covariance matrix, is C = PTN~1P (6.3) where P is the new pointing matrix, containing the elements of h placed in the locations corresponding to the pixels pointed to in the absolute power measurements. —* The problem is to determine both h and where spatially the measurements that are weighted by the elements of h should be located. This is specified in the pointing matrix, P. The solution depends on N and, in general, seems difficult. Chapter 6. Discussion 119 6.5 Noise & Correlation Removal The streaking, cross-hatching and large-scale structure present in the final iteration map are indicative of a non-diagonal pixel-pixel noise covariance matrix. The final iteration map is the least-squares optimal solution to Px = M (6.4) To improve upon the pixel correlation properties in the final iteration map, it is necessary to change the type of optimal solution that is desired. One such approach is through the use of priors in Bayes' equation. An example of a prior is a uniform probablility distribution for the pixels which is zero for negative pixel values. It is unclear which priors will provide an appropriate optimal map solution, or even how to determine them. Perhaps a more difficult alternative is to re-measure the measurements with a strategy which yields a more diagonal pixel-pixel noise covariance matrix. Consider m = Ps + n (6.5) where in is the existing vector of double-difference measurements, P is the unmod-ified pointing matrix, s the sky and ft is the noise. Let Nm Nm rhj = cijiPi • s + ajini (6-6) i = i i = i where aji is the entry on the jth row and the ith column of A, a re-measurement matrix. Pi is the ith row of the pointing matrix, and rhj is the jth linear combination of existing measurements. Let P] = J2ajiPi = (AP)j (6.7) i=l and AT, n i = i then M = Ps + fi (6.9) Chapter 6. Discussion 120 Forming the new pixel-pixel covariance matrix C = PTN-1P (6.10) Expanding equation 6.10 gives C = PTAT [ANAT] 1 AP (6.11) Thus, the task is to determine A such that C is diagonal. The first conclusion to draw, is that if A is invertible, C = C. Solving for A looks tricky; however, we know of one solution. Consider the pointing matrix of an absolute power experiment. Every row of P contains a single one, and, for the sake of discussion, let N contain off diagonal terms consistent with slowly varying structure, or j noise. We perform a single-difference experiment by taking differences between the absolute power measurements. This removes the constant offset in the absolute-power measurements, and reduces the slowly varying structure, reducing the off-diagonal terms in N. At the same time, the new PTP matrix becomes more band-diagonal. Note that in this example, A is a timeline differencing matrix, and is singular.1 If the measurement strategy is very good, then there are many measurements involving a single pixel with many other pixels and the extra diagonalization of N is accomplished with a nearly diagonal PTP. The determination of A in this example, is a restatement of the problem posed in the previous section. Thus, the task is to determine the singular matrix A that will form linear com-binations of the existing SHADES data to produce a maximally diagonal C matrix. 1A is an JVm x Nm square matrix. Appendix A. Jordan Form 121 Appendix A Jordan Form For any square matrix, A, there exists a matrix M such that J = M~lAM (A.l) where the columns of M consist of linearly independent "generalized eigenvectors" whose construction is given, for example, in [21]. The result is that the decomposi-tion: A = MJM~l (A.2) exists for any square matrix A and the matrix J consists of the eigenvalues of the A along the diagonal and ones and zeros just above the diagonal. The rest of the matrix is zero. Let A equal the characteristic matrix introduced in Section 3.4.1. We are inter-ested in the behaviour of An = (MJAT 1)" (A.3) = MJnM~l (A.4) with large n. Consider the case where J in the decomposition of A in equation A.2 is / Ao 1 0 0 \ 0 Ao 0 0 0 0 Ai 1 V 0 0 0 Ai / (A.5) Appendix A. Jordan Form 122 where An and Ai are the two eigenvalues of A. Now, J = J ( \ 2 A 0 2A0 0 0 \ 0 \ 2 A 0 0 0 0 0 A? 2Ai \ 0 0 0 A? / / 3A2, 0 0 \ 0 Ag 0 0 0 0 A? 3A2 V 0 0 0 A? / (A.6) (A.7) By induction the n t h power of J is, Jn = ( \ n A o nxr1 0 0 \ 0 0 0 0 0 A" n A ^ 1 \ 0 0 0 A" / (A.8) For the elements of A to converge, the absolute eigenvalues of A must be less than one which is the same stability criterion developed in Chapter 3. Taylor expanding the function f(x) = (l-x)n (A.9) about x = 1, and neglecting higher order terms, obtain, f(x) = l n 2 (A.10) Thus, for absolute eigenvalues less than, but close to one, the non-zero off-diagonal terms in A.8 decay with increasing n in the same fashion as the diagonal terms. For absolute eigenvalues closer to zero, the decay of the off-diagonal terms with increasing n is slightly slower than the decay of the diagonal terms with increasing n. Bibliography 123 Bibliography [1] E. N. Archibald, T. Jenness, W. S. Holland, I. M. Coulson, N. E. Jessop, J. A. Stevens, E. I. Robson, R. P. J. Tilanus, W. D. Duncan, and J. F. Lightfoot. On the atmospheric limitations of ground-based submillimetre astronomy using array receivers. MNRAS, 336:1-13, October 2002. [2] R. Arendt, D. Fixen, and H. Moseley. Dithering strategies for efficient self-calibration of imaging arrays. ApJ, 2000. [3] I. Aretxaga, D. H. Hughes, E. L. Chapin, E. Gaztahaga, J. S. Dunlop, and R. J. Ivison. Breaking the 'redshift deadlock'- II. The redshift distribution for the submillimetre population of galaxies. MNRAS, 342:759-801, July 2003. [4] Andrew W. Blain, Ian Smail, R. J. Ivison, J.-P. Kneib, and David T. Prayer. Submillimeter galaxies, astro-ph, 2002. [5] C. Borys. Phd thesis: A sub-millimetre survey of dust enshrouded galaxies in the hubble deep field region. 2002. [6] C. Borys, S. Chapman, M. Halpern, and D. Scott. MNRAS, 2002. [7] James Dunlop. Shades proposal: The scuba/blast wide-field extragalactic sur-vey. 2002. [8] D. T. Emerson, U. Klein, and C. G. T. Haslam. A multiple beam technique for overcoming atmospheric limitations to single-dish observations of extended radio sources. A&A, 76:92-105, June 1979. [9] G. Hinshaw, C. Barnes, C. L. Bennett, M. R. Greason, M. Halpern, R. S. Hill, N. Jarosik, A. Kogut, M. Limon, S. S. Meyer, N. Odegard, L. Page, D. N. Bibliography 124 Spergel, G. S. Tucker, J. L. Weiland, E. Wollack, and E. L. Wright. First-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Data Processing Methods and Systematic Error Limits. ApJS, 148:63-95, September 2003. [10] D. H. Hughes, I. Aretxaga, E. L. Chapin, E. Gaztanaga, J. S. Dunlop, M. J. De-vlin, M. Halpern, J. Gundersen, J. Klein, C. B. Netterfield, L. Olmi, D. Scott, and G. Tucker. Breaking the 'redshift deadlock'- I. Constraining the star for-mation history of galaxies with submillimetre photometric redshifts. MNRAS, 335:871-882, October 2002. [11] D. Johnstone, C. D. Wilson, G. Moriarty-Schieven, J. Giannakopoulou-Creighton, and E. Gregersen. Large-Area Mapping at 850 Microns. I. Optimum Image Reconstruction from Chop Measurements. ApJS, 131:505-518, Decem-ber 2000. [12] M. R. Nolta, E. L. Wright, L. Page, C. L. Bennett, M. Halpern, G. Hinshaw, N. Jarosik, A. Kogut, M. Limon, S. S. Meyer, D. N. Spergel, G. S. Tucker, and E. Wollack. First Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Dark Energy Induced Correlation with Radio Sources. ApJ, May 2003. [13] Press, W. and Teukolsky, S. and Vetterling, W. and Flannery, S. Numerical Recipes in C. Cambridge University Press, second edition, 1992, cl988. [14] Ross, S. Introduction to Probability Models. Academic Press, seventh edition, 2000, cl972. [15] M. Rowan-Robinson. The Star Formation History of the Universe: An Infrared Perspective. ApJ, 549:745-758, March 2001. [16] M. Rowan-Robinson. Photometric redshifts in the Hubble Deep Fields: evolu-tion of extinction and the star formation rate. MNRAS, 345:819-833, November 2003. Bibliography 125 [17] S. Scott. Observing the scuba half degree extragalactic survey (shades). July 2003. [18] S. E. Scott, M. J. Fox, J. S. Dunlop, S. Serjeant, J. A. Peacock, R. J. Ivison, S. Oliver, R. G. Mann, A. Lawrence, A. Efstathiou, M. Rowan-Robinson, D. H. Hughes, E. N. Archibald, A. Blain, and M. Longair. The SCUBA 8-mJy survey - I. Submillimetre maps, sources and number counts. MNRAS, 331:817-838, April 2002. [19] S. Serjeant, J. S. Dunlop, R. G. Mann, M. Rowan-Robinson, D. Hughes, A. Ef-stathiou, A. Blain, M. Fox, R. J. Ivison, T. Jenness, A. Lawrence, M. Longair, S. Oliver, and J. A. Peacock. Submillimetre observations of the Hubble Deep Field and Flanking Fields. MNRAS, 344:887-904, September 2003. [20] D. Sivia. Data Analysis A Bayesian Tutorial. Oxford University Press, 2003, cl996. [21] G. Strang. Linear Algebra And Its Applications. Harcourt Brace Jovanovich College Publishers, third edition, 1988, cl976. [22] T. Webb, S. Eales, S. J. Lilly, D. L. Clements, L. Dunne, W. Gear, K. Adel-berger, M. Brodwin, H. Flores, S. Foucaud, O. Le Fevre, H. McCracken, A. Shapley, C. Steidel, and M. Yun. Dusty Star Forming Galaxies at High-Redshift: The Canada-UK Deep Submillimeter Survey. Bulletin of the Ameri-can Astronomical Society, 34:571—h, December 2001. [23] N. Wright. ApJ, 1996. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0091468/manifest

Comment

Related Items