Table Of ContentDraft version November 13, 2015
PreprinttypesetusingLATEXstyleemulateapjv.05/12/14
MODELING ATMOSPHERIC EMISSION FOR CMB GROUND-BASED OBSERVATIONS
J. Errard30,3, P.A.R. Ade27, Y. Akiba16, K. Arnold14, M. Atlas14, C. Baccigalupi17, D. Barron14, D. Boettger5,
J. Borrill3,30, S. Chapman9, Y. Chinone16,13, A. Cukierman13, J. Delabrouille1, M. Dobbs24, A. Ducout11,
T. Elleflot14, G. Fabbian17, C. Feng32, S. Feeney11, A. Gilbert24, N. Goeckner-Wald13, N.W. Halverson2,6,15,
M. Hasegawa16,31, K. Hattori16, M. Hazumi16,31,19, C. Hill13, W.L. Holzapfel13, Y. Hori13, Y. Inoue16,
G.C. Jaehnig2,15, A.H. Jaffe11, O. Jeong13, N. Katayama19, J. Kaufman14, B. Keating14, Z. Kermish12,
R. Keskitalo3, T. Kisner3,30, M. Le Jeune1, A.T. Lee13,25, E.M. Leitch4,18, D. Leon14, E. Linder25, F. Matsuda14,
T. Matsumura33, N.J. Miller21, M.J. Myers13, M. Navaroli14, H. Nishino19, T. Okamura16, H. Paar14, J. Peloton1,
D. Poletti1, G. Puglisi17, G. Rebeiz7, C.L. Reichardt29, P.L. Richards13, C. Ross9, K.M. Rotermund9,
5 D.E. Schenck2,6, B.D. Sherwin13,20, P. Siritanasak14, G. Smecher24, N. Stebor14, B. Steinbach13, R. Stompor1,
1 A. Suzuki13, O. Tajima16, S. Takakura22,16, A. Tikhomirov9, T. Tomaru16, N. Whitehorn13, B. Wilson14,
0 A. Yadav14, O. Zahn25
2
Draft version November 13, 2015
v
o ABSTRACT
N Atmosphere is one of the most important noise sources for ground-based cosmic microwave back-
ground(CMB)experiments. Byincreasingopticalloadingonthedetectors,itamplifiestheireffective
2
noise, while its fluctuations introduce spatial and temporal correlations between detected signals. We
1
present a physically motivated 3d-model of the atmosphere total intensity emission in the millimeter
and sub-millimeter wavelengths. We derive a new analytical estimate for the correlation between
]
M detectors time-ordered data as a function of the instrument and survey design, as well as several
atmospheric parameters such as wind, relative humidity, temperature and turbulence characteristics.
I Using an original numerical computation, we examine the effect of each physical parameter on the
.
h correlations in the time series of a given experiment. We then use a parametric-likelihood approach
p to validate the modeling and estimate atmosphere parameters from the polarbear-i project first
- season data set. We derive a new 1.0% upper limit on the linear polarization fraction of atmospheric
o
emission. We also compare our results to previous studies and weather station measurements. The
r
t proposed model can be used for realistic simulations of future ground-based CMB observations.
s
a
[
1. INTRODUCTION therefore reach higher angular resolution, but must deal
2 with residual atmospheric effects which can be an im-
Observations of the cosmic microwave background
v portant limitation to their ultimate performance. Emis-
(CMB)areauniqueprobeofthefundamentalphysicsat
1 sionfromwatervaporanddioxygenmoleculesdominates
work in the early universe. Ongoing and upcoming high
1 atmospheric emission at millimeter and sub-millimeter
sensitivity observations of CMB temperature and polar-
9 wavelengths. This effect can be significantly minimized
izationanisotropiesprovideconstraintsontheproperties
7 by observing in so-called atmospheric windows, where
of the very early Universe, for example through observa-
0 emission levels are greatly reduced. However, even this
tionsofCMBprimordialpolarizationfromlargetointer-
.
1 mediate angular scales. In addition, through the precise residual emission inevitably increases the optical power
0 characterizationoftheimpactofgravitationallensingon incident on the detectors and therefore their (photon)
5 noiselevel, e.g., Arnold(2010). Inaddition, theinhomo-
smaller scales, the analysis of CMB polarization brings
1 geneousdistributionofwatervapormolecules, drivenby
informationaboutlargescalestructures,andcorrespond-
: complex mechanisms, which depend on the properties of
v ing constraints on neutrinos masses and species and on
theatmosphereaboveagivenobservationsite, resultsin
i the dark energy, e.g., Planck Collaboration et al. (2013);
X bothtemporalandspatialvariationsofthereceivedopti-
Das et al. (2013).
calpower. Theamplitudeandcorrelationofatmospheric
r Ground-based experimental efforts have been playing
a fluctuations depend on both the scanning strategy and
andwilcontinuetoplayaprominentroleintheexploita-
the properties of the atmosphere at the time of obser-
tion of this potential (e.g. The POLARBEAR Collabo-
vation, such as the wind direction and speed. If treated
rationetal.(2014);BICEP2Collaborationetal.(2014);
as an additional noise-like component, this atmospheric
Naess et al. (2014); Crites et al. (2014); Abazajian et al.
contamination results in an additional colored and spa-
(2013); Wu et al. (2014)). As compared to space-borne
tially correlated signal in the time stream of any specific
instruments such as WMAP1 or Planck2, or planned
detector. This would potentially affect not only the ef-
next-generation space mission such as COrE3, Lite-
fective noise level of the entire instrument, obtained as
BIRD4 and PIXIE (Kogut et al. (2011)), ground-based
a combination of the noise levels of all its detectors, but
experiments can deploy larger primary reflectors and
alsotherateatwhichthisnoiseleveldecreaseswhenthe
number of deployed detectors grows.
1 http://map.gsfc.nasa.gov/
2 http://www.esa.int/Our Activities/ Theatmosphericemissionisexpectedtobelargelyun-
Space Science/Planck polarized, e.g., Hanany & Rosenkranz (2003); Spinelli
3 http://www.core-mission.org/ et al. (2011) and therefore appropriate hardware solu-
4 http://litebird.jp/eng/
2
tions, such as the dual polarized TES (Arnold (2010)), 2. ATMOSPHERIC TRANSMISSION AND
may be used to mitigate its impact on the polarization- EMISSION
sensitive observations of the CMB. Nevertheless, even if
The atmosphere is only partly transparent at millime-
such measures are implemented, the atmospheric emis-
ter wavelengths. Water vapor is responsible for most
sion will in general contribute to the detected polarized
of the continuum absorption in the 100GHz–1THz fre-
signals due to instrumental limitations such as instru-
quencyrange. Inaddition,verystrongabsorptionoccurs
mental polarization, imperfect half-wave plates or fre-
at frequencies within the broad oxygen absorption band
quency bandpass mismatch (e.g., Shimon et al. (2008);
around60±10GHz,aswellasaroundtheoxygenlineat
Errard (2012)). Moreover, total intensity sensitivity
119 GHz and the water vapor lines at 22, 183, 325, and
is crucial for calibration, which can impact, for exam-
380GHz. This leaves only four main atmospheric win-
ple, the performance of component separation based
dows available for CMB observations from the ground:
on a parametric approach, e.g. Stompor et al. (2009).
below 50GHz, and around 95, 150 and 250GHz. Ob-
Therefore, though the focus of this work is on total in-
servations are also possible at about 340GHz, but only
tensity measurements, its conclusions are relevant for
when the atmosphere is exceptionally dry, i.e., for a low
polarization-sensitive observations alike.
Precipitable Water Vapor (PWV) value. At higher fre-
Modeling of the atmospheric effects is complex.
quencies,atmosphericopacityisincompatiblewithsensi-
Fluctuations of the atmospheric optical depth generate
tiveground-basedCMBobservations,evenfromthebest
emission with amplitude and scale which both depend
observing sites like the Atacama desert.
on the observation site (dryness, air density, etc.) and
on the time of observation (temperature, pressure,
PWV [mm] 0.25 0.5 1.0 2.0
etc.). Moreover, wind can displace atmospheric struc-
time below 5% 25% 50% 65%
tures, introducing hard-to-model non-stationary effects.
Nonetheless,ithasbeenshown,Lay&Halverson(2000);
Bussmann et al. (2005), that in the case of a given TABLE 1
experimental setup (focal plane, scanning strategy, Percentage of time below the listed amount of total
precipitable water vapor, for the Llano de Chajnantor in
etc.), and with help of some simplifying assumptions
Chile.
useful models can be developed and implemented.
These previous efforts have used the 2d frozen screen
approximation to characterize observations above the Table 1, which summarizes the percentage of time at
south pole. In contrast, we employ a 3d model, derived Atacama when the PWV is below a certain amount,
as a generalization of the model of Church (1995) and, shows that the typical amount of water vapor is about
for the first time, apply it to the data gathered at the 1 mm, the atmospheric conditions being better about
Atacama Desert. The modeling proposed here requires 50% of the time5. We use the Atmospheric Transmis-
significant computational resources, but is more general sionsatMicrowaves(ATM)code(Pardoetal.(2001))to
and permits, for instance, accounting for near- and compute the atmospheric transmission T, defined as the
far-field beam regimes. frequency-dependent ratio I/I0 between the radiation I
receivedbyadetectorandtheradiationI abovetheat-
0
The paper is organized as follows. In section 2, we mosphere. Atthezenith,fromtheALMAsitelocatedat
describe the total emission of the atmosphere, its im- 5040 m elevation in the Atacama desert in Chile (Llano
pactonthephotonnoiseofaground-baseddetectorasa deChajnantor),wegetthetransmissioncurvesshownin
function of frequency, and the typical emission law of Fig. 1.
atmospheric emission fluctuations due to water vapor
inhomogeneity. In section 3, we derive a new general
expression for the auto- and cross-correlation induced
by atmosphere between two detectors of a given focal
plane geometry. We present in section 4 several results
obtained with the model, computed using Quasi Monte
Carlo algorithms, to illustrate the impact of the various
parameters on the properties of the atmospheric signal.
In section 5, we compare the modeling predictions with
realCMBdatasetsfromthefirstseasonof polarbear-
i observations. In addition, we derive estimates of the
atmosphere model parameters, and compare them with
measurements taken by a nearby weather station. We
alsocompareourresultswithmeasurementsmadeinLay
& Halverson (2000) above the Atacama desert, and esti-
Fig. 1.—AtmospherictransmissionfromtheAtacamaplateauat
mate an upper limit for the polarization fraction of the
thezenithfordifferentamountsofprecipitablewatervapor. This
atmosphere. Finally, we conclude and discuss our re-
isobtainedusingtheATMcode,Pardoetal.(2001).
sults in section 6. As a complement, we introduce in ap-
pendix B a faster method to estimate an approximation For different zenith angles, the optical depth τ has to
of the correlated noise between detectors, encapsulated be modified for the airmass m(90◦−el), so that the at-
inafullbinnednoisecovariancematrix,andtoefficiently
simulate realistic noise time streams. 5 https://almascience.nrao.edu/documents-and-tools/
overview/about-alma/atmosphere-model
3
mospheric transmission is antenna-coupled TES (Arnold (2010)), the considered
time stream would effectively be the sum of the time-
T =exp(−τ)=exp(−m(90◦−el) τ ). (1)
0 ordered data (TOD) coming from two orthogonal anten-
The airmass m(90◦−el) ranges from 1 at zenith (el = nas confined within a focal plane pixel.
90◦) to about 40 at the horizon (el = 0◦), and τ is
0 3.1. Atmosphere contribution to antenna temperature
the optical depth at the zenith. The airmass can be
computedforanyzenithangleusingthefollowingfitting Spatial inhomogeneities in the atmosphere emission
formula by Kasten & Young (1989): can get imprinted in the TOD when the instrument line
of sight scans across them. In addition, the wind blows
1 structuresthroughthebeam. Theseinhomogeneitiesare
m(θ)∝ . (2)
cosθ+0.50572(96.07995−θ[deg])−1.6364 due to turbulences along the line-of-sight, and are some-
times assumed to be concentrated in a two-dimensional
We note that for high elevations we have m(90◦−el) (cid:39) layer: atmospheric contamination is then approximated
1/sin(el), as expected for a plane-parallel approxima- as a screen of brightness temperature, moving with the
tion. wind (Lay & Halverson (2000); Bussmann et al. (2005);
The instantaneous noise of current CMB detectors is Sayers et al. (2010)). Furthermore, these fluctuations
dominated by the statistical fluctuations in photons ar- are usually assumed to be frozen as the involved turbu-
riving at the detectors. Consequently instrument envi- lent processes are much smaller than the displacements
ronment, specifically ground or space, plays a critical imposed by the wind (Taylor (1938)).
role. On the one hand, for space-borne detectors, the ThemodelproposedbyChurch(1995)isnotrestricted
total background in the frequency range of interest for to a single turbulent layer, but instead sees the atmo-
CMB observations is the sum of the CMB blackbody sphere as a continuum three-dimensional medium that
emissionandtheemissionfromtheoptics,dominatedby depends on water vapor distribution, wind speed, tem-
the emission from telescope mirrors, perature, etc. as one moves away from the telescope,
along the line of sight. Our work consists in the imple-
I(ν)=B (T )+ε (ν)B (T ), (3)
ν CMB tel ν tel mentation and exploitation of this 3d modeling, hence
where B (T) is the emission of a thermal blackbody at significantly differing from previous studies. Based on
ν
temperature T. For a space mission such as Planck, the Kolmogorov model of turbulence (Tatarskii (1961)),
T (cid:39)40 K, and on the basis of the pre-launch measure- the power spectrum of the fluctuations in a large three-
tel
ments performed on the Planck reflectors (Tauber et al. dimensional volume is
(2010)), we can model the emissivity of a two-reflector
P(k)∝kb (6)
space-borne telescope as
where k is the wavenumber with units m−1, and b =
(cid:16) ν (cid:17)0.5
ε (ν)=0.002× . (4) −11/3. As a consequence, one can show that the con-
tel 140GHz taminated time streams have power law spectra, with
On the other hand, for ground-based experiments, the an index related to the Kolmogorov spectrum, typically
temperature of the telescope and cryostat window is of P(f)∝fb,e.g.,Marriage(2006);Errard(2012). Follow-
theorderof280Kandthetotalemissivityistypicallyof ing Church (1995), a coherence distance scale is defined
the order of a few percent. In addition, the emission of over which the Kolmogorov spectrum holds — the outer
theatmospherecontributestothetotalbackgroundwith scaleLo correspondstothetypicaldistancethatsatisfies
an additional term E(ν) of that condition.
In the following reasoning, we use the geometry de-
E(ν)=[1−T(ν)] Bν(Tatm), (5) picted in Fig. 2. An imaginary central detector of the
focal plane observes along a vectorˆr =ˆr (t), defined in
where T (cid:39) 280 K is the atmosphere temperature. s s
atm spherical coordinates by φ and θ (equivalently azimuth
Thus,at150GHz,thetotalbackgroundemissionis∼ 20
and elevation). All other detectors are assumed to ob-
K on the ground while it is ∼ 1 K in space.
serve in directions slightly off from the central line of
In addition to this background emission, the emission
sight,ˆr i(t)=ˆr (t)+δri, with δri defined by the design
of the atmosphere fluctuates as a function of time and s s
of the focal plane.
pointingdirection. Amodelofthesefluctuations, dueto
Each detector has an associated beam pointing in a
inhomogeneitiesofitsphysicalproperties,isdiscussedin
given direction ˆr i, with an effective area B(ˆr i,r) such
the next sections. s s
that, for a monochromatic detector observing at wave-
length λ:
3. ATMOSPHERE PHYSICAL MODEL
In this section, we present the modeling of the atmo- B(ˆri,r)= 2λ2|ˆris·r|2 ×exp(cid:20)−2(r2−(ˆris·r)2)(cid:21) (7)
spheric contamination occurring for ground-based ob- s πw2(ˆri ·r) w2(ˆri ·r)
s s
servations, through the derivation of an analytical ex-
where
pression based on the auto- and cross-correlation be- (cid:115)
tween detectors time streams, initially proposed by (cid:18)λˆri ·r(cid:19)2
w(ˆri ·r)=w 1+ s (8)
Church (1995). Since the atmosphere mainly radiates s 0 πw2
non-polarized light, the detectors considered through- 0
out the following sections are assumed to give measure- with w0 is the beam waist given by
ments of the sky total intensity, i.e., the Stokes param-
λ
eter I, expressed in K. For example, in the case of w ≡ , (9)
0 πθ
b
4
point, i.e.,
1 dV
dT (r)= B(ˆri,r)α(r,λ)T (r) (10)
ant λ2 s phys r2
where the transmission T(λ) mentioned in Eq. (1) is re-
lated to the absorption α(r,λ), expressed in m−1, by
T(λ)=e−τ(λ) =e−(cid:82)0rα(r(cid:48),λ)dr(cid:48). (11)
T (r) in Eq. (10) is the physical temperature of the
phys
given volume of atmosphere. It is convenient to assume
an adiabatic atmosphere, so that T depends linearly
phys
on altitude z ≡r·ˆz,
T (r)=T (z)
phys phys
Fig. 2.—Adetectorofthefocalplaneobservesalongavectorˆrsi (cid:18) z (cid:19)
(line of sight), defined in spherical coordinates φi and θi (equiva- =T 1− (12)
lently azimuth and elevation: θi =0 corresponds to the horizon). ground z
atm
ˆrs is the pointing direction of an imaginary central detector onto
the focal plane. The scan strategy, ss, can be expressed as the withT ∼280Kthetemperatureatthegroundlevel
ground
evolutionofˆrswithrespecttotime,ss=dˆrs/dt. Windisavector and z ∼104 m a typical height which depends on the
orthogonal to the z-axis and is characterized by its direction φW atm
anditsnormW. observation site. Note that Eq. (12) is not expected to
be correct above the tropopause, generally located be-
tween 10 and 20 km, e.g. Hobgood (1993). We keep
this assumption through the paper as our results do not
depend significantly on z , the atmospheric emission
atm
being also weighted by the water vapor density profile.
ThislatterdecreasesmuchmorerapidlythanT with
phys
altitude as explained in the following paragraph.
3.2. Analytical expression for the auto- and
cross-correlation between detectors
By integrating Eq. (10) over r, the total antenna tem-
peratureT asmeasuredbyadetectoricanbewritten
ant
as
T i ≡T (ˆr i(t))
ant ant s
1 (cid:90) dr
= B(ˆr i(t),r)α(r)T (r). (13)
λ2 r2 s phys
Fig. 3.— Illustration of three functions used in the modeling of Moregenerally,usingEq.(13),thecorrelationbetween
the correlation Citjt(cid:48): the water vapor column, χ2(z,z), Eq. (18); two samples t and t(cid:48), measured by two given detectors i
the temperature profile, T , Eq. (12); and the Gaussian beam and j, is defined as (cid:104)T i (t)T j (t(cid:48))(cid:105), expressed as:
phys ant ant
widthw(z),Eq.(8).
Ctt(cid:48) ≡(cid:104)T i (t)T j (t(cid:48))(cid:105)≡(cid:104)T (ˆr i(t))T (ˆr j(t(cid:48)))(cid:105)=
ij ant ant ant s ant s
1 (cid:90) dr(cid:90) dr(cid:48)
B(ˆr i(t),r)B(ˆr j(t(cid:48)),r(cid:48))
λ4 r2 r2 s s
θ beingthebeamopeningangleinradians. Asexpected,
fobr large distances from the telescope, i.e., |ˆri ·r|(cid:29)w , ×(cid:104)α(r)α(r(cid:48))(cid:105)Tphys(r)Tphys(r(cid:48)) (14)
s 0
w(ˆris ·r) has a slant asymptotic with slope w0, leading where the average (cid:104)·(cid:105) is taken over realizations of the
to B(ˆri,r) ∝ θ2. It is computationally expensive to sky. As in Church (1995), we reduce the correlation (cid:104)·(cid:105)
s b
consider a realistic behavior for the beam effective area, between two given points r and r(cid:48) in the atmosphere to
and this was usually neglected in previous works (Lay & the (cid:104)α(r)α(r(cid:48))(cid:105) term. Note also that we neglected the
Halverson (2000); Bussmann et al. (2005); Sayers et al. wavelengths in (cid:104)α(r,λ)α(r(cid:48),λ(cid:48))(cid:105). We will limit ourselves
(2010)). Fig.3showsthebehavioroftheGaussianbeam tothemonochromaticcaseinthisworkbutwillmention
width w(z) given in Eq. (8), compared to two other im- the λ(cid:54)=λ(cid:48) possibility in section 6.
portant functions used in this modeling and introduced QuantitiesinEq.(14)dependontheatmosphereprop-
below, the water vapor density (Eq. (18)) and the tem- erties(hiddeninthe(cid:104)α(r)α(r(cid:48))(cid:105)×T (r)T (r(cid:48))term)
phys phys
perature profile (Eq. (12)). and on the experimental design and the operation of the
Given the expression for B(ˆri,r), the contribution telescope(included in theB(ˆr i(t),r)B(ˆr j(t(cid:48)),r(cid:48))term).
s s s
dT to the antenna temperature of a small element dV At this point, if no wind is assumed, the only time de-
ant
of atmosphere located at a distance r is proportional to pendenceinEq.(14)isencodedinthescanningstrategy
theeffectiveareaofthetelescopebeamasseenfromthat i.e.,ˆr (t).
s
5
used for χ naturally depends on the site of observa-
2
tion. For example, water vapor density column at South
Polefallsoffquitesignificantlyataltitudesoutsidea300-
2,000 m layer (Bussmann et al. (2005)). We come back
to this difference in appendix A.
We assume in this work that wind behaves such that
it does not mix atmosphere between different altitudes.
Wind is thought to displace atmospheric structures, in
particular the fluctuations encoded in the (cid:104)α(r)α(r(cid:48))(cid:105)
correlation term. We suppose that wind affects the
χ (∆r) term of the correlation, Eq. (16), and leads to
1
anewtimedependenceinthislatterχ (∆r)→χ (∆r(cid:48))
1 1
such as
∆r(cid:48) ≡∆r(cid:48)(∆t)=∆r−W∆t (19)
with ∆t ≡ |t(cid:48)−t| and where W is a vector describing
wind direction and amplitude. W is supposed to not
Fig. 4.—ComparisonofχKolm.(∆r)(gray)givenbyEq.(17)and depend on the altitude z as in Church (1995) or Lay &
1
χ1(∆r) (orange), given in Eq. (16), normalized to 1 for ∆r = 0 Halverson (2000), i.e.,
m. Here we set Lo = 100 m, κmin = (Lo)−1 = 0.01 m−1 and
κmax=100m−1. W=Wxxˆ+Wyyˆ ≡Wwˆ, (20)
where wˆ corresponds to the unitary vector along the di-
To derive an analytic expression for the correlation rection of the wind. It is possible to refine Eqs. (19,20)
term, (cid:104)α(r)α(r(cid:48))(cid:105), we follow Church (1995) and assume and include a dependence of W on the altitude. We do
that not explore this possibility in this work and rather focus
(cid:104)α(r)α(r(cid:48))(cid:105)=χ (|r−r(cid:48)|)χ (z,z(cid:48)) (15) on a simple effective wind, described by a constant wind
1 2 norm W =|W| [m s−1] and a single direction φ [deg].
W
whereχ istheturbulentpartofthecorrelation,effective Thisapproximationturnsouttoreasonablydescribereal
1
for lengths satisfying |r−r(cid:48)| ≤ L . In this work, as in observations, as detailed in section 5.
o
Church (1995), χ is approximated by a Gaussian, i.e., AmplitudetermsinEqs.(16)and(18)ensurethatthe
1
(cid:18) ∆r2(cid:19) dimension of Citjt(cid:48) is K2. In this work, we set χ01 = 1.0
χ (∆r)=χ0exp − , (16) and χ0 = 1.0 m−2. The global amplitude of the corre-
1 1 2L2 2
o lation is parametrized by an effective ground tempera-
where ∆r ≡|r−r(cid:48)|. As depicted in Fig. 4, this is nearly ture, T0, expressed in Keff. Note that measurements of
the true ground temperature T , as measured by a
equivalent to Kolmogorov correlation, given by ground
weather station, and the estimation of T from Ctt(cid:48), cf.
1 (cid:90) κmax 0 ij
χKolm.(∆r)∝ κ−11/3sin(κ∆r)κdκ (17) section 5, can give an estimate for the conversion factor
1 ∆r κmin between Keff and K.
where κ and κ are proportional to the inverse of
min max 3.3. Comparing the 3d and 2d approaches
typical sizes of the turbulences, respectively their outer
(energy provision) and inner (energy dissipation) scales. PreviousstudiessuchasLay&Halverson(2000);Buss-
In the case depicted in Fig. 4, the relative difference mann et al. (2005); Sayers et al. (2010) used the approx-
between the integration of Eqs. (16) and (17) over ∆r imation of a 2d frozen screen. Appendix A shows that
is about 10%. In the following, we will use the Gaus- the 3d modeling introduced above is general enough to
sian correlation given in Eq. (16) for computational effi- reproduce the 2d screen predictions. In fact, parameters
ciency. While the quantitative change due to the Gaus- can be tuned so that the angular correlations from both
sian approximation is less than 10% for Ctt(cid:48), the Gaus- approachesarenearlyequivalentonsmallangularscales.
ij However, this would assume a stepped function for the
sianshapeforχ drivesthephysicalinterpretationofthe
1 water vapor column, which is not supported by observa-
typical scale L .
o tionsabovetheAtacamadesert(Giovanellietal.(2001)).
The second term in Eq. (15), χ , only depends on the
2 Unlikemeasurementsfromsouthpole,estimatedparam-
altitude and can be interpreted as the water vapor dis-
eters from a 2d screen approach might be difficult to in-
tribution. This latter is assumed to be an exponential
terpret: inparticular,thetypicalaltitudehandthickness
decreasing with altitude z,
∆h of the turbulent layer would not have clear physical
(cid:18) z+z(cid:48)(cid:19) meanings in the case of distributed atmospheric turbu-
χ (z,z(cid:48))=χ0exp − , (18)
2 2 2z lences.
0 Moregenerally,notethattherearetwosignificantrea-
which seems to be a reasonable approximation as in- sons for using the 3d formalism. First, the beams are
dicated by some observations made above the Chilean overlappingwheremuchoftheemissionisoccurring,and
desertofAtacama(Giovanellietal.(2001)). Forillustra- this can be the case for large aperture instruments. Sec-
tion, Fig. 3 shows χ (z,z) compared to the beam width ond, the scale length can be smaller than the thickness
2
w(z) and the temperature profile T (z). The model of the layer where the atmosphere is emitting.
phys
6
Despite the caveats mentioned above, we will com- changetotheintegralvalue. Inaddition,weusethespa-
pare the effective brightness of the atmospheric fluctu- tial locality of the functions involved in the integrand to
ations, as estimated from both approaches. This com- speed-up the computation, namely B(ˆr (t),r), Eq. (7),
s
parisonisreasonablesincetheobservedbrightnessofat- and χ (∆r), Eq. (16): we reduce the integration bounds
1
mospheric fluctuations, B2, should not depend on the in Eq. (14) to scales within the beam size or the turbu-
ν
assumed modeling. At small angular scales, identifying lent typical length. For example, we reduce dr(cid:48) to a 3d
Eqs. (A1) and (A6) leads in particular to T2 ∝B2, with volume of radius
0 ν
a coefficient of proportionality depending on L , h, ∆h,
o
the angular scale γ and zatm. In section 5, we will have ˆrs(t)·r ± Lo, (23)
estimates for T and L derived from real observational
0 o which is effectively the 1σ region of the correlations in-
data. Although significant uncertainties on h and ∆h
duced by the turbulences. In addition, the (θ,φ) and
above the Atacama desert make any quantitative com-
(θ(cid:48),φ(cid:48)) variables are integrated around
parison between the 2d and 3d approaches quite risky,
estimates of atmosphere brightness derived from both (cid:18) (cid:18) (cid:19) (cid:18) (cid:19)(cid:19)
L L
modelings are presented in section 5.4. θ ±max θ , o ,φ ±max θ , o . (24)
s b r s b r
4. IMPLEMENTATION AND AN ILLUSTRATIVE
Moredistantlocationsaresuppressedeitherbythebeam
EXAMPLE
term B(ˆr (t),r), or by the turbulence term, χ . Outside
s 1
In this section, we illustrate the 3d model introduced the ‘’1σ” volume defined by Eqs. (23) and (24), inte-
insection3bycomputingelementsofCtt(cid:48), Eq.(14), for grand contributes to less than 5% of the total integral
ij
various atmosphere conditions, as parametrized by the in Ctt(cid:48). Finally, the implementation takes advantage of
ij
atmospheric fluctuations typical scale Lo, the wind W, symmetry relations such as Ctt(cid:48) =Ctt(cid:48).
the typical altitudes z and z , as well as the effective ij ji
0 atm Foragivensetofphysicalparameters{L ,T ,...},the
groundleveltemperatureT . ContrarytoChurch(1995), o 0
0 implementationofthealgorithmusedforthisarticleper-
no simplifications nor approximations are imposed on
Eq. (14), and we perform the full computation of the forms the computation of a single element Citjt(cid:48) in ∼ 1
six integrals to evaluate each Ctt(cid:48) element. ms on one CPU. A basic parallelization of the algorithm
ij is easy to implement and efficient to compute the entire
4.1. Assumptions: scan strategy, focal plane layout set of Citjt(cid:48) elements.
Forillustrationpurposes,wesetascanningstrategyof
the telescope, ˆr (t), to be constant in elevation (θ (t) = 4.3. Illustration of the modeling for a detector
s s
θ )andsweepinginazimuth. Asaroughapproximation auto-correlation
s,0
of usual scanning strategies adopted by CMB ground- Assuming the specific observing strategy given in
based experiments, ˆrs(t) is assumed in this section to Eq. (21), Fig. 5 depicts Ctt(cid:48) for the central detector
follow a simple cosine function, ij
i=j =0asafunctionofbotht(cid:48) (x-axisassumingt=0)
∆az and one of the parameters involved in the atmosphere
φ (t)=φ + cos(2πf t+ψ), (21)
s s,0 2 scan conditions description (y-axis, all the other parameters
are set to their fiducial values). C0t(cid:48) corresponds to
where we set ∆az to be the angular size of the scan 00
the time-domain autocorrelation of the central detector.
(usuallyoftheorderofacoupleofdegrees)andf the
scan Moreprecisely,Fig.5depictsthequantityC0t(cid:48)/C00,fid
scan frequency, defined as ij ij
where C00,fid is computed for the chosen fiducial pa-
ss[deg s−1] ij
f [Hz] ≡ , (22) rameters { L = 300m, z = 2000m, z = 40,000m,
scan ∆az[deg] FWHM = 3o.5arcmin, ss0 = 1 deg s−1a,tm∆az = 5 deg,
where ss = |ss| is the norm of the scanning speed. For φs0 = 0 deg, θs,0 = 45 deg, ψ = 0 deg } and { W = 25
simplicity, we only consider a single detector in this sec- m s−1, φW =0 deg } in the windy cases (panels (5) and
tion, pointing along the central line of sight ˆr (t) of the (6)only). AcommonfeatureacrossallpanelsinFig.5is
s
telescope. Throughout this paper, this detector will be the presence of nearly scan-synchronous signals, induced
denoted i = j = 0. Note that in section 5, in which by atmospheric structures, roughly acting as a sky-like
the modeled atmospheric contamination is compared to signal for low wind speed and/or short period of time.
realdata,Ctt(cid:48) iscomputedusingtherealpointingofthe We describe below the effect of changing a parameter
ij individually on the auto-correlation function:
experiment {φ=az(t),θ =el(t)}.
4.2. Implementation • Turbulence typical length, Lo: panel(1)shows
that an increase of this length results in a larger
Significant computational resources are required to
width of the scan-synchronous features. In the
compute Eq. (14). We use a Quasi Monte Carlo inte-
limit of an infinitely large turbulent typical scale,
gration method to numerically estimate the six-integral
all the detectors would be 100% correlated, for an
included in the computation of a single Ctt(cid:48) element.
ij infinitely long time if there is no wind. Contrarily,
Thenumberofquasi-randomsampleshasbeenoptimized a L → 0 m reduces the nearly scan-synchronous
o
and turned out to be ∼ 102 for the computation of a features, bringing the averaged correlation to neg-
singleCtt(cid:48) element. Moresamplingsleadtosub-percent ligible levels.
ij
7
Fig. 5.— Normalized auto-correlation of a detector as a function of time, i.e., C 0t(cid:48)/C 00,fid as defined in Eq. (14), for various
00 00
atmospherephysicalparameters. Thefiducialparametersarechosenas{Lo=300m,W =0.00ms−1,φW =0.00deg,∆az=5.00deg,
ss = 1.00degs−1 }. The observing experiment is assumed to scan the sky in azimuth (here 6 back-and-forth movements in 30 s), at a
givenelevation. Adetaileddescriptionofeachpanelcanbefoundinthemaintext.
• Correlation amplitude, T : as mentioned ear- spheric signal: the lower is z , the more important
0 0
lier, thisisaneffectivegroundtemperatureandits is the contribution of atmospheric fluctuations at
squareistheglobalnormalizationofthecorrelation lower altitudes.
matrix Ctt(cid:48), Eq. (14), in the case of monochro-
ij • Wind properties: results are depicted on pan-
matic detectors. Therefore T simply scales the
0
els (5) and (6). The first one shows the effect of
entirecorrelationfunction,andisnotillustratedin
Fig. 5. wind speed on C0t(cid:48), for a constant wind direc-
00
tion φ = 0 deg, cf. the geometry illustrated in
W
• Telescope scanning properties, ss and ∆az: Fig. 6. The larger is the wind speed, the more
they impact fscan, Eq. (22), and consequently the rapidly harmonics of C0t(cid:48) are suppressed. Since
frequencyofthescansynchronousfeaturesinCtt(cid:48). windisdisplacingatmosp00hericturbulences,anyde-
ij
As shown in panel (2), in the absence of wind, the tector cannot observe the same sky features at a
largerthescanspeedis,themorescan-synchronous given azimuth pointing. From Eq. (14), having a
features are present. As depicted in panel (3), the non-zero W leads to Citjt(cid:48) ∝ e−(W|t(cid:48)−t|)2/2L2o: for-
larger ∆az is (for a fixed scan speed), the fewer gettingaboutthecross-terminthisexpression,this
scan-synchronousfeaturesappear. Alongandslow is why the correlation is suppressed for long time
(respectively short and rapid) movement of the intervals |t(cid:48)−t|(cid:29)L /W.
o
telescope modulates the atmospheric structures to
low (respectively high) time streams frequencies. Panel (6) shows the effect of a wind direction φW
change at a given wind speed, W = 25ms−1.
• Typical altitudes z and z : they have simi- The geometry of the scanning strategy with the
atm 0
lar impact on the correlation. Panel (4) of Fig. 5 wind direction is depicted in Fig. 6: as defined
shows the effect of changing z on C0t(cid:48): smaller in Eq. (21), the assumed scanning strategy is cen-
water vapor column results in0smaller00correlation teredaround(cid:104)φ(t)(cid:105)t =0 degsincewetookthefidu-
amplitude, with a mild impact due to the form of cial φs,0 = 0 deg. On one hand, pseudo-harmonic
χ ,cf. Eq.(18). Inaddition,notethattheposition structures of C0t(cid:48) are efficiently suppressed for
2 00
of the water vapor densest region would also play φ ∈ {−180, 0, +180} deg. For the chosen scan-
W
a role on the width of the nearly scan-synchronous ningamplitude∆az =5 deg,thesewinddirections
features. The χ term acts as a weight for the are roughly orthogonal to the scanning direction,
2
contribution of turbulent structures to the atmo- which means that the wind efficiently brings new,
8
Fig. 6.— Scan strategy of the telescope is a cosine function,
Eq.(21),centeredonφs,0=0. Thewindisassumedtobeparallel
tothex−y planeandhaveafiducialdirectionφW =0.
and hence uncorrelated, atmospheric fluctuations
across the telescope line of sight. On the other
hand, even if atmospheric structures are shifted
and suppressed as a function of ∆t, significant
features remain in the time streams for the case
φ = ±90deg. The physical interpretation of
W
this effect is that the telescope resamples already
scanned structures on the sky: the φ = ±90deg
W
windshiftsatmosphereturbulencesinanearlypar- Fig. 7.—Autopowerspectra,obtainedastheFouriertransform
allel direction to the telescope scans. Specifically, ofC000t(cid:48),shownhereforvariousLo(upperpanel,assumingW =25
given our assumptions for the scanning strategy, a m−1andφW =0deg),variouswindspeeds(middlepanel,assum-
φwWith=th−e9l0indeeogfwsioguhltdmmootivoenadtumroinspghtehreic1sstt,ru3crdtu,r.e..s, i(nlogwφerWpa=nel0,dasesgumanindgLWo ==23500mms−)1anadndvaLroio=us30w0inmd).directions
(2ith+1)subscanofthetelescope. Onthecontrary,
a φ =+90deg would move fluctuations with the
W
line of sight motion during the 2nd, 4th, ..., 2ith To summarize, we presented the dependence of the
subscan. Bussmann et al. (2005) gives a detailed autocorrelation regarding the turbulence typical size, on
description of this phenomenon. theobservingstrategy,onthewindpropertiesandonthe
water vapor column typical height. In the limit of a low
The remarks above hold for cross-correlations, i.e., for wind speed, the model predicts the presence of nearly
i(cid:54)=j case. However, wealwayshaveCt=t(cid:48) ≤Ct=t(cid:48) and scan synchronous features in the detector time streams.
i(cid:54)=j i=j
the wind can impact correlations amplitude and shape However, the width of these features strongly depends
between detectors, depending on their specific positions on the typical size L of the turbulence, and depends
o
on the focal plane, as well as on the scanning strategy. weakly on the typical altitudes z and z . Moreover,
0 atm
Complementary to Fig. 5 which illustrates the cor- theamplitudeofthesefeaturesusuallydecreasesrapidly
relation in the time domain, Fig. 7 shows the Fourier as a function of the elapsed time |t(cid:48) −t|, mostly driven
transform of several C0t(cid:48), i.e., the auto-spectra, for by the wind speed. The pseudo-periodicity of these
00
variouswindspeedandwinddirections. First,theupper features is modulated by the wind direction relative to
panel of Fig. 7 shows the effect of various turbulence the scanning strategy.
scales L at a fixed wind direction and speed (W = 25
o
m s−1 and φ = 0 deg): similarly to the observations The next section will quantitatively compare the
W
predictions of the modeling described above with re-
derived from C0t(cid:48), the auto spectra present features
00 cent observations of the CMB ground-based experiment
appearing at fscan = 0.2 Hz. Large Lo results in polarbear-i.
more power at low frequencies while small L leads
o
to a broader spectrum, with spread power at high
5. QUANTITATIVE COMPARISON BETWEEN
frequencies. Second, the middle panel shows the effect
THE MODELING PREDICTIONS AND REAL
of various wind speeds at a fixed wind direction and L
o CMB DATA SETS
(φ = 0 deg and L = 300 m). The larger the wind
speWed, the broader aore the peaks in frequency space. Inthissection,wefitrealPolarbear-idatasetswith
As mentioned previously, wind displaces atmospheric the modeling introduced in section 3 and derive distri-
structures and the detector does not resample these butionsofatmosphereparameterssuchastypicalturbu-
structures at a constant frequency f . Third, the lence scale and wind speed. Our approach here is quite
scan
lowerpanelshowstheeffectofvariouswinddirectionsat similar to Bussmann et al. (2005), but differs in two as-
a fixed wind speed and L : as suggested by the correla- pects: first, we use the 3d modeling of the atmosphere
o
tion behavior in time domain, wind direction modulates and second, we estimate the physical atmospheric pa-
the typical frequencies of atmospheric contamination rameters through the maximization of a parametric like-
features. lihoodbasedonthecorrelationmatricesCtt(cid:48). Aquanti-
ij
9
a new set of atmospheric parameters (except z and
0
z as explained below) to be estimated and therefore
atm
a new likelihood to be maximized.
For a given polarbear-i CES, we estimate the full
real datasets’ covariance matrix as
(cid:28)(cid:68) (cid:69) (cid:29)
Dtt(cid:48)≡ d(m,S)d(n,S) (25)
ij i j
{m,n}∈{t,t(cid:48)} S
where d(m,S) denotes the calibrated total intensity data
i
from a pixel i of the focal plane, at a given time sam-
ple m—sample which comes from a group of subscans
Fig. 8.— Drawing of the scanning strategy of ∼ 8-hour S. The average (cid:104)·(cid:105)S is taken over several periods made
polarbear-iobservations. Itiscomposedofseveralconstantele- of consecutive subscans. In each of these groups, we
vationscans(CES). consider an even number of subscans, so that left-going
(right-going) subscans are always averaged with left-
tativecomparisonwithpreviousresultsderivedfromthe going (right-going) subscans only.
2d modeling approach above the Atacama desert (Lay Given the instrument specifications, the complete
& Halverson (2000)) is also presented. This analysis Dtt(cid:48) would have a size (n ×n )2 ∼ (108)2.
involves an estimator of the full covariance of the real ij pixels samples∈S
It would therefore be challenging to quickly perform
data sets introduced in section 5.1, and a parametric
operationsonthisdataset. Tosimplifythecomputation
maximum-likelihood framework detailed in section 5.2.
of the real data covariance matrix Dtt(cid:48), we use some
Results of the fit are summarized in section 5.3. ij
properties of the atmospheric contamination. For
5.1. Analysis of POLARBEAR-I data sets reasonable scanning strategies, since the atmosphere is
evolving at rather low temporal (≤ 1 Hz) and spatial
polarbear-i is a CMB polarization imaging experi-
frequencies (≥0.5deg), a few detectors among the focal
ment6 locatedontheChajnantorplateau,at17,000feet,
plane with low frequency time streams can capture
nearto(andhencelookingthroughasimilaratmosphere
enough information about atmospheric fluctuations
as) the Atacama Cosmology Telescope7, the ALMA8 properties. We hence choose “cardinal” pixels across
and the APEX9 observatories. It started its operations the focal plane (e.g., at the center and edges), and we
in early 2012 and has performed observations of three downsample the time streams of each of these pixels in
patches of the sky at 150 GHz, leading to an effective ordertolightenDtt(cid:48) andthereforespeeduptheanalysis
coverage of 30 deg2 of the sky. The polarbear collab- ij
without degrading the parameters estimation.
oration has recently published the results of the first
season of observation (2012-2013), and demonstrated
First the compression of {m, n} into temporal bins
measurements of CMB B-modes at sub-degree angular {t, t(cid:48)}inEq.(25)correspondstothedownsamplingpro-
scales (Ade et al. (2014); POLARBEAR Collaboration cess. Second, dm in Eq. (25) is also the result of an
et al. (2013); The POLARBEAR Collaboration et al. i
averageoftheTODoversomeregionsofthefocalplane,
(2014)).
preciselyoverdetectorsclosetothe“cardinal”pixelsi.e.,
polarbear-i is composed of 1274 bolometric detectors
(637 dual-polarization pixels sensitive to total intensity) dm →(cid:104)dm(cid:105) , (26)
i k k∈neighbor pixels to pixel i
sampled at 190.73 Hz (Kermish et al. (2012)). As
schematically depicted in Fig. 8, polarbear-i CMB where the neighbor pixels k are the nearest observing
detectors to a “cardinal” pixel i. Complementary to
observations are composed of constant elevation scans
(denoted CES hereafter), each of them lasting ∼ 15 Eq. (25), the error bar on each Dtt(cid:48), denoted ∆Dtt(cid:48),
ij ij
minutes. Each CESisitself composedof ∼100subscans is taken as the standard deviation of
i.e., right-to-left and left-to-right going scans, with a (cid:68) (cid:69)
∆φ ∼ 5 deg amplitude. The motion of the telescope d(m,S)d(n,S) ,
i j
{m,n}∈{t,t(cid:48)}
is a bit more complex than the one considered in
Eq. (21), with consecutive acceleration/deceleration computed over groups of subscans S. We compute a
(turnarounds) and constant velocity periods. The scan- single D matrix for each polarbear-i CES.
ning strategy (θ (t), φ (t)) is an input to our modeling,
s s
and we use the true motion of the Polarbear-i tele- Atmospheric contamination encapsulated in Dtt(cid:48) can
ij
scope in this section. The instrument periodically scans sometimes be approximated by a stationary process,
inazimuthand,asassumedinsection4,foragivenCES meaning that Dtt(cid:48) (cid:39) D0∆t, with ∆t ≡ t(cid:48) − t. This
we are in a θ =constant case. In this work, we fit each ij ij
s pseudo-stationarity is exploited in Appendix B, where
CES independently: each of them is considered as an
anapproximationofDtt(cid:48) isintroducedforrapidcharac-
independent realization of the atmospheric fluctuations, ij
terizationofcoloredandcorrelatednoise. However,wind
6 http://bolo.berkeley.edu/polarbear/ usually breaks stationarity since it modulates the atmo-
7 http://www.princeton.edu/act/ sphericsignalperiodicity,whichdependsonthepointing
8 http://www.almaobservatory.org/ of the telescope and on the direction of the wind. For
9 http://www.apex-telescope.org/ example, in the case of a wind nearly parallel to the
10
scanning strategy (φ = ±90 deg in Fig. 6), the de- the physical parameters p = {L , T , W,...} describing
W o 0
tectorsofaleftgoingfocalplaneorofarightgoingfocal the atmosphere conditions, with the real-data matrices,
plane would not show the same correlation pattern, cf. Dtt(cid:48), estimated from the polarbear-i first season data
ij
section 4.3. sets.
The comparison of modeled and real covariances has
beentackledinotherastrophysicalcontexts,forexample,
5.2. Parametric likelihood and implementation
in the case of astrophysical foregrounds separation from
We want to construct a likelihood to quantitatively multi-frequency data sets. Following the formalism by
androbustlycomparethepredictionsoftheatmospheric Pham & Cardoso (2001), the negative log-likelihood L
contamination modeling, Ctt(cid:48) = Ctt(cid:48)(p), function of can be written as
ij ij
−2log(L(p))∝(cid:88)(cid:26)tr(cid:18)Ctt(cid:48)(p)(cid:16)Dtt(cid:48)(cid:17)−1(cid:19)−log(cid:20)det(cid:18)Ctt(cid:48)(p)(cid:16)Dtt(cid:48)(cid:17)−1(cid:19)(cid:21)−constant(cid:27) (27)
ij ij ij ij
t,t(cid:48)
(cid:39)(cid:88)(cid:26)tr(cid:20)(cid:16)Ctt(cid:48)(p)−Dtt(cid:48)(cid:17)(cid:16)Dtt(cid:48)(cid:17)−1(cid:16)Ctt(cid:48)(p)−Dtt(cid:48)(cid:17)(cid:16)Dtt(cid:48)(cid:17)−1(cid:21)(cid:27), (28)
ij ij ij ij ij ij
t,t(cid:48)
whereEq.(28)isthequadraticapproximationofEq.(27) We implement and test the algorithms on NERSC
(Delabrouilleetal.(2003)). ForthecomputationofL(p), systems10, and use a simple parallelization scheme in
note that Ctt(cid:48)(p) is given by the integration of Eq. (14) which a single process performs the analysis of one
for agiven siejt ofparameters p, cf. section4, and Dtt(cid:48) is CES; this corresponds to the computation of a Ditjt(cid:48) and
ij
the covariance matrix computed from polarbear-i real the estimation of the parameters {Lo, W, φW, T0} by
maximizing the likelihood given in Eq. (28). For this
data sets, cf. section 5.1.
work, we decide to fix z = 2000 m and z = 40,000
For the estimation of D, we consider time streams 0 atm
d m over groups of 6 consecutive subscans, defining S m as they turn out to be badly conditioned: the de-
ini Eqs. (25) and (26): averages (cid:104)·(cid:105)S are then taken over pendence of Citjt(cid:48) on these quantities is quite weak, and
all the available groups of 6 consecutive subscans within make the atmospheric model degenerate, in particular
a CES, corresponding to ∼ 35 s long time streams. As with the amplitude T (section 4.3). We therefore allow
0
mentioned earlier, note that S is composed of an even only {L , W, φ , T } to vary. In addition, scan speed,
o W 0
numberofsubscans, sothatleft-goingscansareonlyav- elevationandazimuthinvolvedinEq.(14)aresetbythe
eraged with left-going scans, cf. Bussmann et al. (2005). real observations of the telescope, i.e., we compute the
Similarly, right-going scans are only averaged over right- modeled C for the detectors location used to compute
going scans. Prior to Eqs. (25) and (26), we remove the D, using the same scanning strategy, and at the same
globalmeansignalfromalltheTODoverthecourseofa binned samples {t,t(cid:48)}.
given CES. A ground-template, obtained as the estima-
tion of an azimuth-fixed signal, is removed as described Finally, as illustrated in Fig. 9, the estimation of the
by The POLARBEAR Collaboration et al. (2014). In atmospheric modeling parameters is performed in two
addition, we choose 5 cardinal pixels (one at the cen- steps: we first compare D0∆t with precomputed tem-
ij
ter and four on the edge, forming a square shape). The plates C0∆t, estimated for various atmospheric condi-
ij
average in Eq. (26) is taken over the 5 closest working
tions and saved on disk. The minimization of −2log(L)
pixels to each of the cardinal pixels. We downsample
givesaroughestimateoftheadjustedparameters,which
the time streams from 190.73 Hz to 2 Hz, in order to
areusedasastartingpointforthesecondstep. Wemaxi-
speed up the algorithms without missing the main fea-
mizethelikelihoodexpressedinEq.(28),basedthistime
tures of the atmospheric contamination. As shown in
on the matrices Dtt(cid:48) and Ctt(cid:48), using gradient informa-
Fig. 10 and described in section 5.3, these features are ij ij
essentially pseudo scan-synchronous occurring at a fre- tion in a truncated Newton algorithm, e.g., Wright &
quency of ∼0.25 Hz in the case of polarbear-i. Nocedal (2006).
In addition, we provide the optimization algorithm
with the derivatives of the likelihood, with respect to 5.3. Results
the considered parameters, namely ∂L/∂Lo, ∂L/∂W, Fig. 10 shows the correlations between the five
∂L/∂φW and ∂L/∂T0. These derivatives are semi- “cardinal” pixels, D0∆t + constant ± ∆D0∆t (colored
analyticallycomputedusingEq.(28),andincludederiva- ij ij
points), and the corresponding estimated predictions
tives of the modeled covariance Ctt(cid:48) defined inEq. (14).
ij from the atmospheric model, C0∆t + constant (solid
Havingsemi-analyticalderivativesofthelikelihoodleads lines), with constant ≡ −minij(cid:0)D0∆t(cid:1). Only four
to a significant speed-up of the optimization. ij
10 National Energy Research Scientific Computing Center ma-
chineshttp://www.nersc.gov/
