Table Of ContentCoaxing Cosmic 21cm Fluctuations from the Polarized Sky using m-mode Analysis
J. Richard Shaw,1,∗ Kris Sigurdson,2 Michael Sitwell,2 Albert Stebbins,3 and Ue-Li Pen1
1Canadian Institute for Theoretical Astrophysics,
60 St. George St., Toronto, ON M5S 3H8, Canada
2Department of Physics and Astronomy, University of British Columbia, Vancouver, BC V6T 1Z1, Canada
3Theoretical Astrophysics Group, Fermi National Accelerator Laboratory, Batavia, IL 60510, USA
Inthispaperwecontinuetodevelopthem-modeformalism,atechniqueforefficientandoptimal
analysisofwide-fieldtransitradiotelescopes,targetedat21cmcosmology. Weextendthisformalism
to give an accurate treatment of the polarised sky, fully accounting for the effects of polarisation
leakage and cross-polarisation. We use the geometry of the measured set of visibilities to project
down to pure temperature modes on the sky, serving as a significant compression, and an effective
first filter of polarised contaminants. As in our previous work, we use the m-mode formalism
4 with the Karhunen-Lo`eve transform to give a highly efficient method for foreground cleaning, and
1 demonstrate its success in cleaning realistic polarised skies observed with an instrument suffering
0 from substantial off axis polarisation leakage. We develop an optimal quadratic estimator in the
2 m-mode formalism, which can be efficiently calculated using a Monte-Carlo technique. This is
used to assess the implications of foreground removal for power spectrum constraints where we
n
findthatourmethodcancleanforegroundswellbelowtheforegroundwedge, renderingonlyscales
a
J k(cid:107) < 0.02 hMpc−1 inaccessible. As this approach assumes perfect knowledge of the telescope,
we perform a conservative test of how essential this is by simulating and analysing datasets with
9
deviationsaboutourassumedtelescope. Assumingnoothertechniquestomitigatebiasareapplied,
we find we recover unbiased power spectra when the per-feed beam width to be measured to 0.1%,
]
O and amplifier gains to be known to 1% within each minute. Finally, as an example application,
we extend our forecasts to a wideband 400–800MHz cosmological observation and consider the
C
implications for probing dark energy, finding a pathfinder-scale medium-sized cylinder telescope
.
h improves the DETF Figure of Merit by around 70% over Planck and Stage II experiments alone.
p
-
o I. INTRODUCTION with GMRT [5], HERA [6], LOFAR [7], MWA [8], MI-
r
t TEoR[9]andPAPER[10]targetingtheEpochofReion-
as Recent years have seen a surge in excitement at the isationandBAOBAB[11],BAORadio[12],BINGO[13],
[ promise of radio cosmology. By using low frequency ob- CHIME[14],EMBRACE/EMMA[15],Ooty[16],Parkes
servations of the 21cm line we can survey the distribu- [17] and Tianlai [18] aiming at the low redshift intensity
1
tionofneutralhydrogenthroughoutlargevolumesofthe mapping era.
v
5 Universe. Radiointerferometersprovideanefficient,and In this paper we will focus on the low redshift, in-
9 cost effective method for doing this. tensity mapping epoch, though most of the results and
0 This transformation of radio interferometers into sur- techniques we describe apply equally well at higher red-
2 vey instruments has been driven by recent technological shift. Observationsattheselowredshiftstargetthesame
1. advances, particularly in the cheap low-noise amplifiers science as spectroscopic galaxy redshift surveys such as
0 required for mobile phones, and the constant progress [19, 20]: measuring Baryon Acoustic Oscillations and
4 ofMoore’slawmakinglarge, highbandwidthcorrelators through them probing the expansion history of the Uni-
1 economical. By correlating a large number of low cost verse [21–23]. However, they are very complimentary,
v: feeds in a compact area we can produce a telescope ide- withradioobservationsprobingalargervolumeathigher
i ally suited for wide-field surveys. redshift, with a completely different set of systematics.
X
There are three main epochs we can observe: low red- To make effective use of this new generation of radio
r shift (z < 4), where we observe the large scale emission interferometers, we must develop new methods of inter-
a
from un∼resolved galaxies, a technique termed intensity preting and analysing their data. Progress has acceler-
mapping [1, 2]; the Epoch of Reionisation (z 6–10) ated in recent years with many developments [24–31].
wheretheneutralIGMiseatenawaybythefirst∼ionising In a previous paper [32] we developed a new techique
sources [3]; andperhapseven theprimordial structure in fortheanalysisofdatafromtheseexperimentscalledthe
the dark ages (z > 30), though observations at these m-mode formalism. This method departs from the usual
very low frequenci∼es (ν < 50MHz) will be extremely interferometric analysis—making no flat-sky or small
challenging [4]. These eras are of huge cosmological im- field approximations—at the expense of being limited to
portance, a fact reflected in the large number of current transit telescopes for which it is an exact treatment. It
and planned experiments targeting 21cm observations, also brings computational advantages by allowing us to
breakthedataintouncorrelatedm-modes,makingitfea-
sible to treat the full statistics of the data. This opens
up the possibility of performing optimal map making,
∗ [email protected] foreground subtraction, and power spectrum estimation,
2
which would be extremely difficult otherwise. and demonstrate its effectiveness on simulated polarised
Perhapstheforemostchallengefacing21cmcosmology skies. In Section VIII we construct an optimal power
is the presence of bright astrophysical radio sources at spectrum estimator in the m-mode formalism, which we
frequenciesbelow1.4GHzwhicharearoundsixordersof use to study the performance of our foreground filter
magnitude brighter than the 21cm signal. This emission (Section IX). In Section X we use this estimator to show
comes mainly from synchrotron radiation, which is spec- how instrumental uncertainties give rise to power spec-
trally smooth, and in principle this allows it to be sepa- trum biases. Finally we forecast the performance of our
ratedfromthe21cmasitisdescribedbyasmallnumber example telescope at measuring the expansion history of
of modes [33] and these can simply be removed. The re- the Universe and constraining the nature of dark energy
mainingmodes,whichhavesignificantspectralstructure (Section XI).
are assumed to be free of contamination. Unfortunately,
thispictureiscomplicatedbytherealitiesofradioobser-
vation: II. POLARISED TRANSIT TELESCOPE
ANALYSIS
Frequency dependent beams lead to mixing of an-
•
gular structure into spectral structure which con-
In this section we develop a fully polarised version of
taminates the foreground clean modes [34]. This
problem, known as mode mixing, means that look- them-modeformalism,anewmethodforanalysingtran-
ing at only the frequency direction of our data is sitinterferometersthatwasfirstintroducedinaprevious
paper [32]. That treatment encapsulates all the essential
insufficient to separate these two signals.
ideasbutavoidstheaddedcomplexityoftrackingthepo-
Synchrotronemissionfromourgalaxyishighlypo- larisation, and is a useful introduction to the full treat-
•
larised, and though the totally intensity is spec- ment given here. Polarised descriptions of full-sky inter-
trally smooth, Faraday rotation by the magnetic ferometry have been given elsewhere (notably [36, 37]),
interstellarmediummeansthatthepolarisedemis- but here we develop the transit telescope limit.
sion is not. Unfortunately, the complicated polari- Any transit telescope can be viewed as a collection of
sationresponseofrealtelescopesirreversiblymixes feeds, fixed relative to the ground frame. Each feed, F
i
some fraction of the polarised sky, introducing sig- measuresacombinationoftheelectricfieldE (nˆ)coming
a
nificant frequency fluctuations [35]. As the emis- fromvariousdirectionsonthesky. Inordertoaccurately
sioncomesfromarangeofFaradaydepths,wecan- treat the polarisation when the response varies over the
not simply de-rotate the emission. sky, we need to be able to keep track of the contribution
from each direction to the electric field at a point. In
Fundamentallytherearestillthesamenumberoflarge
order to do this we define ε as the electric field density
foregroundmodes,modemixingonlymakesthemharder
in a frequency interval dν and solid angle d2nˆ by
to identify. In [32] we developed a foreground removal
techniquebasedontheKarhunen-Lo`eve(KL)transform.
dE =(µ c)1/2ε(nˆ,ν)d2nˆdν . (1)
Thisusesthefullcovariancestatisticsofthecontaminat- 0
ing foregrounds to find an optimal separation from the
With this definition the Poynting flux is conveniently
21cm signal, fully accounting for this mode mixing ef-
written as
fect. However, the technique presented there was lim-
ited in two important ways: no attempt was made to 1
S = E H
address the problem of polarised foregrounds; and it as- p µ c ×
0
sumed that we have full knowledge of the properties of (cid:90)
our instrument, including the full polarised response of = d2nˆd2nˆ(cid:48)dνdν(cid:48)nˆ(cid:10)ε(nˆ) ε(nˆ(cid:48))(cid:11) . (2)
·
theprimarybeam,andanyper-feedamplitudegainsand
phaseshiftsintroducedinthereceiversystem. Inthispa-
Radio emission from the sky is generally incoherent and
per we continue to develop both the m-mode formalism so we can write the correlations of ε explicitly in terms
and KL transform for foreground cleaning, with particu-
of the Stokes parameters
lar emphasis on investigating these two limitations.
a fWulel tsrteaarttmbeyntexotfenpdolianrgistahteiomn-(mSeocdteiofnorImI)a,laisnmd dtoiscguivses (cid:10)εa(nˆ,ν)ε∗b(nˆ(cid:48),ν)(cid:11)= 2λk2Bδ(nˆ −nˆ(cid:48))δ(ν−ν(cid:48))
how the unpolarised approach of [32] is a limiting case (cid:104) (cid:105)
(Section III)). The example telescope we use throughout × PaTbT(nˆ)+PaQbQ(nˆ)+PaUbU(nˆ)+PaVbV(nˆ) , (3)
is described in Section IV, and its harmonic space sen-
sitivity is examined in Section V. Next we take a care- wheretheindicesareoverbasisvectorstransversetothe
ful look at the geometry of the measured m-modes (Sec- line of sight. As in the unpolarised case we are more in-
tion VI), leading us to a technique which both efficiently terested in the brightness temperature on the sky, and
compresses the data and effectively removes polarised so we have written Equation (3) to make that explicit
contamination. We give an overview of the Karhunen- (thus Q, U, and V are polarisation brightness tempera-
Lo`eve scheme for foreground removal in Section VII, tures). The polarisation tensors X are related to the
Pab
3
Pauli matrices (in an orthonormal basis), specifically n . In this work we will assume that the noise from
ij
different antennas, and frequencies are uncorrelated (we
(cid:18) (cid:19) (cid:18) (cid:19)
I = 1 1 0 , Q = 1 1 0 , discuss the statistics in more detail in Appendix A).
Pab 2 0 1 Pab 2 0 1 Wenormaliseourvisibilitiessothattheyarethecorre-
−
(cid:18) (cid:19) (cid:18) (cid:19) lated antenna temperature in the noiseless limit (in par-
1 0 1 1 0 i
PaUb = 2 1 0 , PaVb = 2 i −0 . (4) ticulartheauto-correlationis theantennatemperature).
The factor of two in the definition of the transfer func-
The standard basis to use in spherical geometry are the tion ensures that for an unpolarised sky with uniform
polar and azimuthal directions, θˆ and φˆ, as these allow brightness Tb, the measured autocorrelation Vii = Tb.
Note that in Equation (7) and onwards the symbol V
spin spherical harmonics to be used straightforwardly to
denotes two different quantities, the visibility V (φ) and
decompose the polarisation field. ij
theStokesVskyfieldV(nˆ). Thedistinctionwillbeclear
Any feed on the telescope measures a weighted com-
from the context.
bination of the electric field coming from each direction
The above Equations (7) and (8) are completely ex-
on the sky. In particular we need to keep track of the
act. The general approach to interferometric analysis
antenna’s sensitivity to the orientation of the incoming
is to approximate the above to a 2D Fourier transform,
electric field. We’ll write the measured signal at the i-th
which is valid for small fields of view. For wide-field ob-
feed, as F which is given by
i
servations we can attempt to relax this with techniques
(cid:90) such as mosaicing [38] and w-projection [39] though this
F (φ)= d2nˆAa(nˆ;φ)ε (nˆ)e2πinˆ·ui(φ) , (5)
i i a quickly becomes complicated. In our case we are inter-
estedinaspecificclassoftransit interferometersintended
and is directly proportional to the voltage induced in forsurveys. However,astheseinstrumentsareextremely
the circuit. Here, and onwards, we will implicitly sum wide-field, this approach is limiting. Instead we will try
over the polarisation index a. The antenna reception a different route, restricting our domain to transit tele-
pattern Aa is a vector quantity describing the electric scopes, but otherwise attempting to keep the analysis
i
field response in a given direction. The response A exact.
∝
leff, the effective antenna length (choosing them equal To continue, we decompose into spherical harmonics,
wouldmakeFi betheantennavoltage). WenormaliseA as they are a natural way of representing fluctuations
suchthatthenormalisedantennapowerpatternPn(nˆ)= on the sky. As polarisation is not a scalar field we must
A(nˆ)2 ensuring the solid angle of the beam is expandQandU inspin-2harmonicsY(±2)(nˆ)(theStokes
lm
| | V field transforms as a scalar). This yields
(cid:90)
Ωi = d2nˆ|A(nˆ)|2 . (6) T(nˆ)=(cid:88)aT Y (nˆ), (9)
lm lm
lm
In Equation (5) we have also included an exponential (cid:88)
Q(nˆ)+iU(nˆ)= a(+2)Y(+2)(nˆ), (10)
factor which gives the phase relative to an arbitrary ref- lm lm
erence point. As both this, and the antenna orientation lm
(cid:88)
change with the Earth’s rotation relative to the sky, we Q(nˆ) iU(nˆ)= a(−2)Y(−2)(nˆ), (11)
− lm lm
write them as functions of φ, the rotation angle. lm
The fundamental quantity in radio-interferometry is V(nˆ)=(cid:88)aV Y (nˆ). (12)
the cross correlation between two feeds, the visibility lm lm
V =(cid:10)F F∗(cid:11). Using Equations (3) and (5) we can write lm
ij i j
down exactly what a visibility measures, explicitly keep- The polarised beam transfer matrices also transform as
ing track of the different sky polarisations to give spin fields, and so we decompose them in the same way,
with
(cid:90) (cid:104)
Vij(φ)= BiTj(nˆ;φ)T(nˆ)+BiQj(nˆ;φ)Q(nˆ) BT(nˆ;φ)=(cid:88)BT (φ)Y∗ (nˆ), (13)
ij ij,lm lm
(cid:105) lm
+BiUj(nˆ;φ)U(nˆ)+BiVj(nˆ;φ)V(nˆ) d2nˆ+nij(φ) (7) BQ(nˆ;φ) iBU(nˆ;φ)=(cid:88)B(+2) (φ)Y(+2)∗(nˆ), (14)
ij − ij ij;lm lm
where the beam transfer functions BX encode all the lm
informationabouttheopticsandgeomiejtryoftheinstru- BiQj(nˆ;φ)+iBiUj(nˆ;φ)=(cid:88)Bi(j−;2lm) (φ)Yl(m−2)∗(nˆ), (15)
ment. They are given by lm
(cid:88)
BV(nˆ;φ)= BV (φ)Y∗ (nˆ). (16)
2 ij ij,lm lm
BiXj(nˆ;φ)= Ω Aai(nˆ;φ)Abj∗(nˆ;φ)PaXbe2πinˆ·uij(φ) , (8) lm
ij
Note that we have decomposed with the complex con-
(cid:112)
where Ω = Ω Ω . The measured visibilities are cor- jugates of the spin-harmonics. This allows us to use the
ij i j
rupted by noise, which we include as an additional term orthogonalityofthe(spin)sphericalharmonicstorewrite
4
the visibility equation Equation (7) as productsoftheindividualFouriermodes(witharemain-
ing summation over the l index). This equation fully de-
(cid:88)(cid:104) 1 scribes how the measured visibilities are related to the
V (φ)= BT (φ)aT + B(+2) (φ)a(+2)
ij ij;lm lm 2 ij;lm lm polarised sky that we are observing.
lm Itisworththinkingaboutwhatwearemeasuring. The
1 (cid:105)
+ B(−2) (φ)a(−2)+BV (φ)aV +n (φ). (17) visibility we see is a complex time series, which roughly
2 ij;lm lm ij;lm lm ij corresponds to the signal from the sky modulated by a
complex Fourier mode. In our case the time variable is
Though this has completely transformed the problem
φ, the Earth’s rotation. Taking the Fourier transform
into harmonic space, it will be more convenient if we
of a visibility splits the time series into right and left
change into the conventional E and B mode decomposi-
moving waves (positive and negative m respectively). A
tion as they are real scalar fields [40]. This can be done
correlated beam pointing south of the north pole only
by making the standard substitutions
producesmodesmovinginonedirection(asthebeamon
a(+2) = (cid:0)aE +iaB (cid:1) , (18) the sky is a Fourier mode), however pointing the same
lm − lm lm beam beyond the north pole (that is north of it as de-
a(−2) = (cid:0)aE iaB (cid:1) (19)
lm − lm− lm fined in the ground frame), produces the other modes
as the Fourier mode on the sky moves in the opposite
aswellasthecorrespondingchangesforthebeammatri-
directionwithrespecttotheEarthrotation. Oneimpor-
ces
tant consequence of this is that if we use the freedom to
B(+2) = (cid:0)BE iBB (cid:1) , (20) choosetheorderoftheourfeedpairssuchthatthebase-
ij;lm − ij;lm− ij;lm linevectorspointtowardstheeast,positivem-modesare
B(−2) = (cid:0)BE +iBB (cid:1) , (21) produced below the pole, and negative m-modes come
ij;lm − ij;lm ij;lm
from above. If the primary beam does not extend over
leaving the visibility as the pole only positive m’s are produced, though a small
amount of negative m’s are seen because of the effect of
V (φ)=(cid:88)(cid:104)BT (φ)aT +BE (φ)aE the primary beam.
ij ij;lm lm ij;lm lm
Infact, whilstthepositiveandnegativem-modesmay
lm
(cid:105) be independent measurements they are still observations
+BiBj;lm(φ)aBlm+BiVj;lm(φ)aVlm +nij(φ). (22) of the same sky — for a real field alm = a∗l,−m and thus
both Vij and Vij∗ measure the same harmonics on the
In the above the harmonic coefficients are now all the m −m
sky. Itwillbeusefultochangeournotationtomakethis
transformsofrealscalarfields(theBX arethecomplex
ij;lm fact transparent.
conjugates of the spherical harmonic coefficients). Let us separate out the positive and negative m parts
Giventheperiodicityofthesysteminφ,Fouriertrans-
by defining
forming Equation (22) is an obvious next step
BX,+ =BX n+ =n (25)
(cid:90) dφ ij;lm ij;lm ij;m ij;m
Vij;m = 2πVij(φ)e−imφ . (23) BiXj;,l−m =(−1)mBiXj;∗l,−m n−ij;m =n∗ij;−m (26)
which is valid for m 0. Additionally to prevent double
WecalltheseFouriercoefficients,m-modes,andtheywill countingthem=0m≥easurementweneedtosetBX− =
ij;l0
becomethekeyquantityinouranalysis. Asthevisibility n− = 0. For brevity of notation, we will introduce a
is a complex timestream, the positive and negative m’s ij;0
label α which indexes both the positive and negative m
are independent measurements.
partsofallincludedfeedpairsij,suchthatanyparticular
As the φ dependence simply rotates the functions
α specifies exactly the values of ij, (exactly how α is
about the polar axis the transfer function is trivially ±
packed is unimportant). This gives the final form of the
BX (φ) = BX (φ = 0)eimφ. The integral over the
ij;lm ij;lm m-mode visibility equation that we use as the basis of
exponential factors generates the Kroenecker delta δmm(cid:48) this work,
and removes the summation over m entirely, and we can
write the m-modes as V =(cid:88)(cid:104)BT aT +BE aE
α;m α;lm lm α;lm lm
(cid:88)(cid:104) l
V = BT aT +BE aE (cid:105)
ij;m ij;lm lm ij;lm lm +BB aB +BV aV +n . (27)
l α;lm lm α;lm lm α;m
(cid:105)
+BB aB +BV aV +n . (24) As in [32] we can write this equation in an explicit
ij;lm lm ij;lm lm ij;m
matrixformwhichwillallowustosimplifythenotation.
Though slightly hidden, this a property of the convo- The beam transfer matrices above can be written in an
lution theorem. For a transit telescope the visibility explicit matrix notation
timestream is an azimuthal convolution of the beam and (cid:16) (cid:17)
sky. This means its Fourier conjugate, the m-modes, are BXm =BαX;,mνδνν(cid:48) (28)
(αν)(lν(cid:48))
5
wheretherowindexlabelsallbaseline(α)andfrequency
N
combinations (ν), whereas the column index is over all
W
multipole(l)andfrequencies(ν(cid:48)). Similarlywecandefine
vectors for the visibilities and harmonic coefficients
(v ) =Vν (cid:0)aX(cid:1) =aXν . (29)
m (αν) α;m m (lν) lm
D
To keep track of the different polarisation states we
define the block matrix and vector
aT
B=BT BE BB BV , a= aaBE (30)
a
V
such that
v =Ba+n. (31) FIG.1. Aschematicofacylindertelescope,consistingoftwo
cylinders aligned North-South on the ground. Each cylinder
Thisistheessenceofthem-modeformalism: asimple, is of width W, and has Nfeeds regularly spaced a distance D
apart. In this paper we will only consider cylinders which
linear matrix relation that exactly describes the whole
are touching, making the total width of the array 2W. The
measurement process for a transit interferometer. As we
cylinders are assumed to be long enough that there are no
will discuss in Section V both the number of m-modes
opticaldifferencesbetweenfeedsattheedgeandinthecentre
and the dimensionality of the B matrices is bounded by
of the array.
the physical size of the instrument. This means that we
caneasilyapplyallthestandardtoolsofstatisticalsignal
processing without even remembering that we’re dealing
Under these constraints there is only one relevant linear
with an interferometer. In the following sections we do combination of the four XX, YY, XY and YX visibili-
this with gusto.
ties that is sensitive to the total intensity, the average of
Despite this being an interferometry paper, the uv- the XX and YY visibilities
plane has not been mentioned at all so far. Though it is
prevalent in many interferometric applications, as both 1
V = (V +V ) (34)
an extremely useful aid for physical understanding and u 2 XX YY
for computational efficiency (by virtue of the FFT), the 1 (cid:90) 1
= d2nˆA2(nˆ)e2πinˆ·uT(nˆ)+ (n +n ) .
m-mode formalism does not make use of it. Eschewing Ω 2 XX YY
the uv-plane is part of its power, helping it to work triv-
ially for wide-field analysis, and focusing us on only the Because of the properties of the polarisation matrices,
measured degrees of freedom. However, it comes at the this is not sensitive to the different polarisation modes
costofmakingitdifficulttohaveconcretephysicalinter- Q, U and V, whilst all the orthogonal combinations are
pretations of the process. insensitive to the total intensity T.
In this limit, the combination V is equivalent to the
u
unpolarisedformalismgivenin[32],ifwerelabelthenoise
III. UNPOLARISED LIMIT terms such that n=(n +n )/2. Provided the noise
XX YY
terms are uncorrelated this reduces the power spectrum
In[32]wedevelopedanunpolarisedformalismbecause downbyafactoroftwo—thatistheunpolarisedsystem
it gives a simpler problem to analyse, both conceptually temperature is T =T /√2.
sys,u sys,p
and computationally. However, under certain assump-
tions it is directly equivalent to the full polarised case.
For a telescope with dual polarised antennas, let us IV. CYLINDER TELESCOPES
suppose that we can engineer our telescope optics such
that the field patterns of the two feeds (labelled X and
Cylinder telescopes are interferometric arrays consist-
Y) obey two constraints. First, that their normalised
ing of one or more parabolic cylindrical reflectors. They
power patterns are equal everywhere
have a long history in radio astronomy, with well known
A 2 = A 2 =A2 , (32) facilities like the Molongo Synthesis Telescope [41], and
| X| | Y| the Ooty Radio Telescope [42]. Though advances in am-
and second that their polarisation orientations are or- plifiertechnologymeanttheysteadilylostfavourtodish-
thogonal all over the sky basedinterferometers, interestinthemhasrecentlybeen
revived. Reasons are twofold: the development of cheap,
A A =0. (33) room temperature, low noise amplifiers has dramatically
X · Y
6
antenna beam. Rather than try to accurately solve for
TABLE I. Parameters of the example cylinder telescope.
the beam, we try to capture these two effects. We will
Parameters Value break the model down into the product of two 1D func-
Number of cylinders 2 tions: a function for the E-W direction, calculated by
Cylinder width [m] 20 illuminatingthecylinderwiththedipolebeam,andsolv-
Feeds per cylinder 64 (dual-pol) ing for the diffraction in the Fraunhofer limit; and a N-S
Feed spacing [m] 0.3
function which is just the reflected feed amplitude in the
T [K] 50
sys N-S direction. We will also model the polarisation direc-
Bandwidth [MHz] 400–800
tionasbeingthesameasthatofanunfocuseddipole(in
Channel width [MHz] 2.5
spherical co-ordinates, for a dipole along the polar axis,
Number of Channels 160 (in groups of 40)
Telescope Latitude 45◦ the polarisation direction is θˆ).
First we model the beam amplitude for the unfocused
dipole in the E-plane and H-plane as taking the form
improved sensitivity; and 21cm intensity mapping has (cid:18) ln2 tan2θ (cid:19)
providedanapplicationforwhichtheyareideallysuited. AD(θ; θW)=exp − 2 tan2θ , (35)
W
Intensity mapping requires a large collecting area in
a compact region to achieve high brightness sensitivity, where θ is the full width at half-power of the beam.
W
whichcylinderscanprovidedcheaply. Additionallycylin- For a horizontal dipole mounted a distance λ/4 over a
der telescopes are a cost effective way of surveying large conducting ground plane (see [46, section 4.7]), we can
amounts of sky at high speed [43]. And whilst arrays exactly calculate the widths in the H-plane (θ =2π/3)
H
of dipoles provided a bigger instantaneous field of view, and E-plane (θ 0.675θ ). We use these value for our
E H
≈
the large number of elements required at a fixed angular fiducial beam model, though we will vary them later in
resolution makes the receiver and correlation hardware this paper.
increasingly expensive. In the E-W direction we are solving the Frauhofer
Each cylinder has a parabolic cross section such that diffraction problem of a cylinder feed illuminating an
they focus only in one direction. In the layout we as- aperture of a finite width. This has the solution
sume(seeFigure1), thisgivesalongandandthinbeam
othnethNeorstkhy-,Seoxuttehnddiinrgecntieoanrlybufrtowmhhicohrizisonontloyhaorroizuonndin2 AF(θ;θW,W)∝(cid:90)−WW2 AD(2tan−1(2Wx);θW)e−ikxsinθdx
degrees wide East-West. Feeds are spaced along the axis 2
of each cylinder — when correlated these provide resolu- (cid:90) 1 − ln2 u2 −iπWusinθ
e tan2θW 1−u2 λ du (36)
tion in the N-S direction. As the telescope operates as a ∝
−1
transit telescope this means that the entire visible sky is
observed once per sidereal day. where we have used the fact that for a cylinder with
Inthispaperweillustratethem-modeformalismusing an f-ratio of 1/4 a ray striking a distance x from the
cylindercentrereflectsbyanangleθ =2tan−1(2x/W)=
amediumsizedcylindertelescope,similartotheCHIME
2tan−1u where W is the cylinder width.
Pathfinder. Table I lists the parameters of this example
Putting these components together, our overall beam
instrument.
model can be written as the product of three functions.
For the X feed
A. Beam model
AX(nˆ)=A (sin−1(nˆ xˆ); θ ,W)
a F · E
In the m-mode formalism knowledge of the primary AD(sin−1(nˆ yˆ); θH)pa(nˆ; xˆ) (37)
× ·
beams of our instrument is crucial. In our model we as-
and for the Y feed
sume an arrangement such that at each location there
are two perpendicular dipoles: the X feed where the
dipole is aligned across the cylinder (pointing East), and AYa(nˆ)=AF(sin−1(nˆ ·xˆ); θH,W)
Y feedwherethedipoleisalignedalongtheaxis(pointing A (sin−1(nˆ yˆ); θ )p (nˆ; yˆ) (38)
D E a
North). Inbothcasesthefeedshangbelowaconducting × ·
ground plane which stops the beam spilling above the where the vectors xˆ is a unit vector transverse to the
cylinder (which is assumed to have an f-ratio of 1/4). cylinder, pointing East, and yˆ is along the cylinder,
Solving for the beam on the sky for a feed placed in pointing North. The function pa gives the unit vector
a parabolic cylinder is a complex problem (for one ap- polarisationdirection on the sky for adipole in direction
proach see [44, 45]). Crudely the cylinder acts in two dˆ
ways: in the parabolic direction it focuses the antenna
beam to a diffraction limited beam on the sky; in the pˆ (nˆ;dˆ)= 1 (cid:2)dˆ (nˆ dˆ)nˆ(cid:3) . (39)
orthogonal direction it acts like a mirror, inverting the a (cid:0)1 (nˆ dˆ)2(cid:1)1/2 − · a
− ·
7
where the part in square brackets is the coefficient in a
0.001 0.01 0.1 1 0.0003 0.003 0.03 0.3
spherical harmonic expansion of the plane wave. The
90 amplitude of a spherical harmonic function can be con-
I I P I venientlywrittenintermsofintegralsofBesselfunctions
→ →
60 [47, section 5.4]. For large l 1 we find
(cid:29)
2l+1(cid:90) ∞(cid:20) (cid:18)tsinθ(cid:19)(cid:21)2
30 Y (θ,φ)2 = J J (t)dt
| lm | 4π m 2 2l+1
0
g
NS/de 0 ≈ 2lπJm(lsinθ)2 (43)
30
− where we have used the approximation that
lim J (x) = δ(x n). Combining this with
n→∞ n
−60 Equation (42) shows that−the magnitude of the spherical
harmonic coefficients of a plane wave are
90
− −10 −5 EW/0deg 5 10 −10 −5 EW/0deg 5 10 |alm|2 =8πljl(2π|u|)2Jm(lsinθ)2 . (44)
In particular this shows that the coefficients are effec-
FIG. 2. The primary beam of the cylinder telescope forms
tivelyboundedinatrianglebyl<2π u and m <lsinθ
a long strip on the sky from North to South. This figure | | | |
because of the exponential decay of the Bessel functions
illustrates the transfer into an instrumental Stokes I (XX+
for large order.
YY polarisations), from the total intensity on the sky (left
Thehighestfrequencyfouriermodemeasuredonasky
panel), and from the polarised sky only (right panel). The
red contour in the top panel, marks the half power point of by an individual baseline comes from the maximum dis-
the beam. tance between illuminated areas on the correlated an-
tennas (this is the largest distance from the origin in the
uv-plane). FollowingthroughfromEquation(44),weex-
In Figure 2 we illustrate the on-sky beam for this our
pecttherangeofmeasureablemodestobel<2πd /λ
max
example telescope. We plot the response of an ‘instru-
and m < 2πd /λ, where d is the largest distance
E-W max
mental Stokes I’, constructed from the combination of
assoicated with baselines and d the largest in the E-
XX+YY polarisation,toStokesIandpolarisatedemis- E-W
W direction.
sion on the sky. The reponse to Stokes I on the sky is
Let us consider our cylinder (see Figure 1). A feed
given by
on the cylinder effectively illuminates the whole width
R =(cid:0)Aa Ab +AaAb (cid:1) I . (40) of the cylinder, but a very short distance along its axis.
I→I X X Y Y Pab
This makes the largest E-W distance of all feed pairs
As a measure of the response to polarised radiation we
N W,andthelargestN-SdistanceN D. Intermsof
cyl feeds
use
sphericalharmonicscoefficientsonthesky,wearelimited
R2 = (cid:88) (cid:2)(cid:0)Aa Ab +AaAb (cid:1) P(cid:3)2 . (41) to
P→I X X Y Y Pab
P∈{Q,U,V} 2π(cid:113)
l< (N W)2+(N D)2 , (45)
For a beam with no polarisation leakage, this response λ cyl feeds
is zero. Though our example has no leakage on-axis, 2π
m< N W . (46)
Figure 2 clearly shows that there is significant pickup of λ cyl
polarisation away from the beam centre.
Though this result is correct for a cylinder telescope,
foraninterferometerwithacompactfieldofview,point-
ing away from the celestial equator, it needs modifying.
V. SENSITIVITY IN HARMONIC SPACE
As before the resolution in the E-W direction is deter-
mined by the maximum distance d , however, if the
The geometry of an interferometer on the ground de- E-W
primary beam does not cross the equator this resolution
terminesitsangularsensitivityonthesky,withthetotal
corresponds to a larger fraction of the circle of constant
sizeoftheopticalsystemdeterminingthesmallestscales
declinationatthatpoint. Asthem-modecorrespondsto
that can be measured. This limits the number of har-
the Fourier mode in the azimuthal direction, this means
monic modes on the sky that we are able to measure,
thatthelimitonmisinfactm<2πcosδd /λ,where
reducing Equation (27) to finite sums. E-W
δ is the declination of the point in the primary beam
The set of spherical harmonics that a given baseline is
closest to the celestial equator.
sensitive to can be found by expanding a plane wave on
Tolookatthesensitivityofthetelescopeinmoredetail
the sky
wecancalculatetheFishermatrixofthea coefficients
lm
e2πinˆ·u =(cid:88)(cid:2)4πilj (2π u)Y∗ (uˆ)(cid:3)Y (nˆ) (42) (we discuss the interpretation of Fisher matrices in de-
l | | lm lm tail in Section VIII). For Gaussian noise the likelihood
lm
8
function for the a ’s is
lm Sensitivityin(l,m)space—log10[F−1/2/µK]
8
(cid:18) (cid:19)
1 6
L(a;v)∝exp −2(v−Ba)†N−1(v−Ba) . (47) 400 4 7
2
6
From this we can calculate the Fisher matrix for a par-
ticular m 300 1
5
(cid:28) ∂2 (cid:29)
Fll(cid:48) =− ∂aT∂aT lnL l 200 4
(cid:104) l l(cid:48)(cid:105) 0.5 dth dge 3
Weexpectthatinge=neraBl†TthNis−m1BatTrixll(cid:48)w.illbesingular(a4n8d) 100 0.5 cylinderwi telescopee 12
hencewecannotfindthecovariancematrixofthea co-
lm
efficients by finding F−1. One obvious source of this is 0 0
0 100 200 300 400 500
thattheinterferometerdoesnotseethewholesky—any- m
thingdeclinationlessthanδ = 45◦ isbelowthehorizon
—andthismanifestsitselfasc−orrelatedcombinationsof FIG.3. Sensitivityofthearraytotemperature,derivedfrom
a ’s that we cannot separate. Additionally the angular the inverse of the diagonal elements of the Fisher matrix
lm (F )−1. The plot above shows the log of the sen-
resolution falls off towards the horizon meaning that we (lm)(lm) 10
sitivity in units of µK. The sensitivity to the three remain-
do not have uniform sensitivity across the sky.
ingStokesparametersarelargelyidentical. Thedashedblack
Thisconsequenceofthisisobviousfromsimplycount-
linesmarkthemcorrespondingtotheseparationbetweenthe
ing the degrees of freedom involved. For the example
cylinders, and the total width of the cylinders. As we would
telescope, there are 762 unique baselines each of which expect the sensitivity peaks in m at the zero separation, and
gives a noisy complex measurement of the sky. However, the single cylinder separation. It then falls off rapidly at the
we are sensitive up to lmax 400 for each polarisation, edge of the telescope.
∼
giving 4(l m) complex degrees of freedom on the
max
−
sky,sotheremustbesomecombinationsaboutwhichwe
have no information. V
3
As in general we cannot determine the covariance ma-
trixofthea , wewillusetheFishermatrixitselftode-
lm
scribe the sensitivity. In Figure 3 we show the diagonal P
elements of the Fisher matrix at each m for the exam-
ple telescope, this gives an illustration of the amount of
V
informationwehaveaboutanysphericalharmonicmode. 2
T
T’
VI. SVD PROJECTION
V
For21cmcosmologyweareonlyinterestedinderiving 1
Visibility Space Sky Subspace
real properties of the unpolarised sky. As we shall see
thisisusuallyofmuchlowerdimensionthanthespaceof
measurementsmadebyaninterferometer,leavingalarge
FIG. 4. The information about the sky does not spread
number of redundant degrees of freedom which are just
throughout thespace of visibilities butis containedin a sub-
filledbytheinstrumentalnoise. Eliminatingthesewould
space, a linear combination of the measured signals which
allow us to significantly compress the data space, with- does not span the whole visibility space. Directions orthog-
out losing useful information. In Figure 4 we illustrate onal to this subspace are excited only by the instrumental
the geometry of the measured visibilities. The matrix B noise, and contain no information about the sky. The left
wholly describes the mapping between the sky and the panelillustratesthegeometryofthefullvisibilityspace,show-
measured visibilities, and understanding its structure is ing the sky subspace as a plane. In the right panel we show
the key to isolating the important degrees of freedom. only the sky plane. Within this sky subspace, there are yet
lowerdimensionalsubspacesthatthetotalintensity(labelled
To start with let us concentrate on how to reduce to
T) and polarised (P) signals get mapped to. However, they
only the degrees of freedom on the sky (ignoring their
need not be orthogonal, an effect we must take into account.
polarisations for now). The matrix B tells us how a sub-
One way of treating this is to project onto the space orthog-
space of the spherical harmonics a map into a subspace onal to polarisation (labelled T(cid:48)), this eliminates polarised
in visibility space v. This visibility subspace (shown by contamination at the expense of some sensitivity to total in-
the plane in the Figure 4), is termed the image of B. tensity. This is discussed in detail in the text.
The subspace of visibilities orthogonal to the image, is
9
called the cokernel. The cokernel has no mapping to the As a first attempt we might consider projecting onto
sky, and so measuring this subspace yields no useful in- the image of B , the total intensity transfer matrix,
T
formation. By projecting our data onto the image, we rather than the full B. In Figure 4 this would corre-
remove the cokernel and compress our data by retain- spond to projecting straight onto the T vector, rather
ingonlytherelevantdegreesoffreedom. InFigure4this than just the plane.
correspondstoprojectingontotheplane,eliminatingthe Unfortunately as illustrated in Figure 4, the image of
perpendicular dimensions. thetotalintensityneednotbeorthogonaltothesubspace
The number of retained degrees of freedom is given by containingthepolarisedimage. Thisisamanifestationof
the dimensionality of the image — that is, the rank of polarisationleakage. Inthiscasebydoingthiswelosethe
B — and cannot exceed the number of measured modes abilitytodifferentiatebetweenpolarisedandunpolarised
on the sky. For a single frequency and m, the rank is signals from the sky, resulting in catastrophic leakage of
guaranteed to be less than the total number of spherical polarised foregrounds.
harmonics required to describe the polarised sky, that is A resolution to this problem is to project not onto
4(lmax−m). However, inthecaseofincompleteskycov- the image of BT but to perform another projection, this
erage, we cannot measure all spherical harmonic modes time onto the polarisation cokernel. In Figure 4 this is
independently, and this coupling means that the rank is equivalent to projecting onto the vector T(cid:48). By doing
likely to be reduced to around 4fsky(lmax m), where this we ensure that there is no leakage of the polarised
−
fsky is the fraction of sky observed. skyintoourcompresseddata,attheexpenseofthrowing
Thesenumbersdependonlyonthephysicalsizeofthe awayinformationaboutthetotalintensitysignalthatlies
telescope, and not details of the feed distribution. For in the overlap between the two spaces.
compactinterferometerswithlittleredundancy,thenum- To project out the polarised signal, we first construct
ber of feed pairs rapidly exceeds the rank of the matrix, the polarisation transfer matrix
and so projecting onto the image gives a large computa-
tional saving.
To find the image of B we can use the Singular Value Bpol =BE BB BV , (51)
Decomposition (SVD). However, first we will pre-whiten
the visibilities with respect to the instrumental noise.
This transforms it to be uncorrelated with unit variance then we use this to isolate the polarisation cokernel in
and can be done by multiplying them with N−21, a ma- the sky compressed basis by performing another SVD
trix such that N−21(N−21)† = N−1. As N is positive
definite this factorisation always exists and can be found U†IN−12Bpol =UpolΣpolV†pol . (52)
by Cholesky factorisation or eigendecomposition. This
leaves Equation (31) as As before we separate into the image and cokernel of
this matrix, by dividing up U into U and U
N−12v =N−21Ba+N−21n. (49) respectively. Asbeforethesepaproaltionontpoolt,hIetwosppaocl,eNs
We then take the SVD of the whitened beam transfer isnotexact,butdonethroughanumericalthreshold. By
matrix projecting our dataspace onto the cokernel we achieve
this final compression.
N−12B=UΣV† . (50) Overall we have applied three transformations to our
data:
The matrix U defines the image and cokernel, given by
columnsofUcorrespondingtonon-zeroandzerosingular Whiten the instrumental noise by applying N−12
values respectively. In practice many singular values are •
numerically small but not precisely zero, giving modes Project onto the sky subspace by using U†
whichareeithernon-zerobecauseofnumericalprecision, • I
orsimplycarryaverysmallbutnon-zeroamountofinfor- Project out the polarised sky using U†
mationaboutthesky. Inthiscaseweseparatetheimage • pol,N
and cokernel using a numerical threshold. We partition Combined these define a new basis in which to consider
the columns of the matrix U into two matrices UI and our data. One which strives to preserve as much of the
UN which give the image and cokernel respectively. To relevant information as possible, whilst vastly reducing
compress our data we simply filter with the matrix UI the number of degrees of freedom we must consider. We
to give v(cid:48) =U†N−21v. define our filtered visibility data as
I
While this filtering can yield a large compression, we
should note that it preserves all the information about v¯ =U† U†N−21v . (53)
pol,N I
the sky. However, the cosmological signal we are in-
terested in is purely unpolarised and requires only We can write a modified version of the measurement
∼
fsky(lmax m) modes per frequency and m to describe equation (31) which relates this to the sky signal
−
it. This suggests that we should be able to improve our
compression by around another factor of four. v¯ =B¯ a+n¯ (54)
10
where we have defined SVDspectrum(log10Σ/Σmax)
0
B¯ =U† U†N−12B, (55) 140
pol,N I
n¯ =U† U†N−21n. (56) 120
pol,N I
100 2
As the columns of UI and Upol,N are orthonormal, de −
the instrumental noise still has the identity covariance mo 80
(cid:104)on¯fn¯th†e(cid:105)=neNw¯ =maIp.pIinngFimguarteri5xwB¯e,schloewarltyheillsuinstgrualtairngvatlhuaest SVD 60 -3 4
we are only sensitive to a small number of modes on the 40 -5 −
sky. -2
20
In order to visualise our data we will want to make -1
mapsfromourfiltereddataset. ForGaussiandistributed 0 6
0 50 100 150 200 250 300 350 400 −
instrumental noise it is straightforward to make maxi-
m
mum likelihood maps of the sky as discussed in [32]. As
we have whitened the instrumental noise, our data has a FIG. 5. The singular values of B¯ for the 400MHz channel
likelihood function after removal of the polarised modes. Large singular values
(cid:18) (cid:19) represent modes on the sky that are well measured. We see
(a;v˜) exp 1(cid:12)(cid:12)v¯ B¯a(cid:12)(cid:12)2 (57) that at each m there are less than 100 measured degrees of
L ∝ −2 − freedom from the sky, with the spectrum dropping off very
steeply beyond this. This is a significant saving, before com-
and thus we can solve for the maximum-likelihood solu-
pression there are twice as many modes as there are unique
tion using the Moore-Penrose pseudo-inverse, giving our
baselines, including positive and negative m. In our example
best estimate of the spherical harmonics simply as thereare762uniquebaselines(withoutauto-correlations),so
there would be ∼1500 modes.
aˆ =B¯+v¯ . (58)
As in [32], to make a full map of the sky, we simply use
skywillvarywiththeobservedfrequency,drivenby
this estimator on a per-m and per frequency basis and
theopticaleffectsofusingafixedphysicalaperture
collate the estimates. We can then perform an inverse
or feed spacing. Even if the angular fluctuations
Spherical Harmonic Transform to produce sky maps at
on the sky were frequency independent as we scan
each frequency. As we have projected onto the polarisa-
through in frequency the beam structure changes,
tion cokernel, the data does not contain any information
andthisintroducesvariationsofourmeasurements
about the polarised sky. Combined with the minimum
with frequency.
powerpropertyoftheMoore-Penrosepseudoinversethis
meansthatthepolarisedsphericalharmonicswillbezero. Model uncertainties:
Astrophysical foregrounds are poorly constrained
at the small angular and frequency scales that will
VII. FOREGROUND REMOVAL WITH THE be probed by upcoming 21cm intensity mapping
KARHUNEN-LOE`VE TRANSFORM
experiments. Whilst there exist theoretical and
phenomenological models of this regime, a success-
The foremost challenge for any 21cm intensity map- fulforegroundremovalmethodshouldberobustto
ping experiment is separating the cosmological signal uncertainties in the foreground statistics. Though
from astrophysical contaminants which are around 104– most effort has focused on the uncertainties in the
106 times larger. The primary sources are the diffuse two-point correlations, we must also ensure that
synchrotron emission from our own galaxy and emission higher order moments do not impair our analysis.
from extra-galactic point sources [48]. All significant
foregrounds are expected to be spectrally smooth [33], Given these complications, we would prefer a fore-
however, the 21cm signal decorrelates quickly as each ground removal method to be conservative, throwing
frequency corresponds to a different spatial slice. This away potentially useful information in order to be ro-
gives an opportunity to separate the two. bust to them. It is better to be cautiously correct than
Conceptuallyforegroundremovalissimple—wejustre- precisely wrong.
move the smooth frequency component from our obser- Acceptingthatwemayprefertoloseinformationabout
vations. Unfortunatelytherealityisfarfromstraightfor- the21cmsignalinordertobeunbiasedbyresidualfore-
ward. The large dynamic range between the amplitude grounds, we would still like to perform the best job we
of the foregrounds and the 21cm signal makes several can, requiring that we are
effects extremely problematic.
Statistically Optimal:
Mode mixing: Whateverspacetheforegroundsareremovedinwe
In a real experiment the shape of the beam on the must be able to keep track of the statistics of both
