SHARP - III: First Use Of Adaptive Optics Imaging To Constrain Cosmology With Gravitational Lens Time Delays

1.6 MB·

by  Geoff C. -F. Chen

MNRAS000,1–20(2015) Preprint16August2016 CompiledusingMNRASLATEXstylefilev3.0
SHARP - III: FIRST USE OF ADAPTIVE OPTICS
IMAGING TO CONSTRAIN COSMOLOGY WITH
GRAVITATIONAL LENS TIME DELAYS
Geoff C.-F. Chen,1,2,3 ⋆ Sherry H. Suyu,1 Kenneth C. Wong,1,4
6
1 Christopher D. Fassnacht,3 Tzihong Chiueh,2,5,6 Aleksi Halkola, I Shing Hu,7
0 Matthew W. Auger,8 Leon V. E. Koopmans,9 David J. Lagattuta,10
2
John P. McKean9,11 and Simona Vegetti12
g
u 1Institute of Astronomy and Astrophysics, Academia Sinica, P.O. Box 23-141, Taipei 10617, Taiwan
A 2Department of Physics, National Taiwan University,Taipei 10617, Taiwan
3Department of Physics, Universityof California, Davis, CA 95616, USA

4 4National Astronomical Observatoryof Japan, 2-21-1 Osawa, Mitaka, Tokyo181-8588, Japan
1 5Institute of Astrophysics, National Taiwan University,Taipei 10617, Taiwan
6Centerfor Theoretical Sciences, National Taiwan University,Taipei 10617, Taiwan

] 7Department of Mathematics, National Taiwan University,Taipei 10617, Taiwan
O
8Institute of Astronomy, Universityof Cambridge, Madingley Rd, Cambridge, CB3 0HA, UK
C 9KapteynAstronomical Institute, Universityof Groningen, P.O.Box 800, 9700 AV Groningen, The Netherlands
. 10CRAL, Observatoire de Lyon, UniversitLyon1, 9 Avenue Ch. Andr, F-69561 Saint GenisLaval Cedex, France
h 11Netherlands Institutefor Radio Astronomy (ASTRON), P.O. Box 2, 7990 AA Dwingeloo, The Netherlands
p 12Max Planck Institute forAstrophysics, Karl-Schwarzschild-Strasse 1, D-85740 Garching, Germany
-
o
r
t
s
a Accepted—.Received—;inoriginalform16August2016
[

2
ABSTRACT
v
Accurateandprecisemeasurementsofthe Hubble constantarecriticalfortesting our
1
currentstandardcosmologicalmodelandrevealingpossiblynewphysics.WithHubble
2
3 Space Telescope (HST)imaging,eachstronggravitationallens systemwith measured
1 time delays can allow one to determine the Hubble constant with an uncertainty of
0 ∼7%.Since HST willnot lastforever,we exploreadaptive-optics(AO)imaging as an
. alternative that can provide higher angular resolution than HST imaging but has a
1
less stable point spread function (PSF) due to atmospheric distortion. To make AO
0
6 imaging useful for time-delay-lens cosmography, we develop a method to extract the
1 unknownPSFdirectlyfromtheimagingofstronglylensedquasars.Inablindtestwith
: two mock data sets createdwith different PSFs, we are able to recoverthe important
v
cosmologicalparameters (time-delay distance, external shear, lens mass profile slope,
i
X andtotalEinsteinradius).OuranalysisoftheKeckAOimageofthestronglenssystem
RXJ1131−1231 shows that the important parameters for cosmography agree with
r
a thosebasedonHST imagingandmodelingwithin1-σuncertainties.Mostimportantly,
the constraint on the model time-delay distance by using AO imaging with 0.045′′
resolutionistighterby∼50%thantheconstraintoftime-delaydistancebyusingHST
imagingwith0.09′′whenapower-lawmassdistributionforthelenssystemisadopted.
Our PSF reconstruction technique is generic and applicable to data sets that have
multiple nearby point sources, enabling scientific studies that require high-precision
models of the PSF.
Key words: gravitational lensing:strong – cosmolgy:distance scale – methods:data
analysis – adaptive optics
1 INTRODUCTION
The discovery of the accelerated expansion of the Uni-
⋆ E-mail:[email protected]
verse (Perlmutteret al. 1999; Riess et al. 1998) and obser-
(cid:13)c 2015TheAuthors
2 G.C.-F. Chen et al.
vations of the Cosmic Microwave Background (CMB; e.g., point spread function (PSF) for the lens mass modeling.
Hinshawet al.2012;Planck Collaboration et al.2015)have However, HST’s lifetime is finite1, and the angular res-
established a standard cosmological paradigm where our olution is also limited by its aperture size. Given the
Universe is spatially flat and is dominated by cold dark dozens of time-delay lenses from COSMOGRAIL2 (e.g.
matter (CDM) and dark energy: the so-called flat ΛCDM Vuissozet al. 2007, 2008; Courbin et al. 2011; Tewes et al.
model, where Λ represents a constant dark energy density. 2013b,a; RathnaKumaret al. 2013; Eulaers et al. 2013),
While the CMB provides strong constraints on the param- and hundreds of new lenses to be discovered in the near
eters of this model, a relaxation of the assumptions in this future (e.g., Oguri & Marshall 2010; Agnello et al. 2015;
model,suchasspatialflatnessorconstantdarkenergyden- Chan et al. 2015; Marshall et al. 2015; More et al. 2015),
sity, leads to a strong degeneracy between the cosmological finding an alternative long-term solution for this promising
parameters,particularlythosewiththeHubbleconstantH . method is timely.
0
Therefore, independent and accurate measurements of H One alternative approach is imaging from the ground
0
provide one of the most useful complements to the obser- via adaptive optics (AO), which is a technology used to
vationsof theCMB inconstraining thespatial curvatureof improve the performance of optical systems by reducing
theUniverse,darkenergyequationofstate,andthenumber theeffectofwavefrontdistortions(e.g., Rousset et al. 1990;
of neutrino species (e.g., Hu 2005; Riess et al. 2009, 2011; Beckers 1993; Watson 1997; Brase 1998). In other words, it
Freedman et al. 2012; ?).The recent inferred valueof Hub- aims at correcting the deformations of an incoming wave-
ble constant H = 67.8 0.9kms−1Mpc−1, based on the front by deforming a mirror and thuscompensating for the
0
±
PlancksatellitedataoftheCMBandtheassumptionofthe distortion.Theadvantagesof usingAOimaging are (1) the
flat ΛCDM model, is low in comparison to several direct angular resolution obtained with telescopes that are larger
measurements including those from the Cepheids distance than HST can be higher than that of HST since a per-
ladderwithH =74.3 1.5(stat.) 2.1(sys.)kms−1Mpc−1 fect AOsystem wouldlead toadiffraction limited PSF,(2)
0
(Freedman et al. 2012)±and H =±73.8 2.4kms−1Mpc−1 ground-basedtelescopes are more accessible. Thedisadvan-
0
±
(Riess et al. 2011). If this indication of tension is not ruled tageisthatwedonothaveastablePSFmodelapriori,since
out by systematic effects, then this could indicate new the atmospheric distortion varies both temporally and spa-
physics beyond the standard flat ΛCDM model. Therefore, tiallyacross theimage.Lenstargetstypicallydonothavea
pinningdown theHubbleconstant with independentmeth- nearbybrightstarwithin 10arcseconds,andstarsat fur-
∼
ods is a keyapproach tobetter understandourUniverse. therangular distancefrom thetarget may beinsufficientin
Strong gravitational lensing with time delays provides providinganaccuratePSFmodelgiventhespatialvariation
a one-step measurement of a cosmological distance in the of thePSF across the field.
Universe.Thebackgroundsourceiscomposedofacentrally In HST imaging, we can use the lensing arcs to con-
varying source, such as an active galactic nucleus (AGN), strain thelens mass model by usingthestable PSFof HST
and its host galaxy. The time delays between the multi- toseparatethearcfrom thebrightAGN,butwecannotdo
ple images of the source, induced by the foreground lens, so in AO imaging. The contamination of the AGN light on
are given by ∆t = 1D ∆τ. Here, ∆τ is dependent on the thelensingarcsinAOimagingmakesitdifficulttoconstrain
c ∆t
geometry and the gravitational potential of the lens sys- thelensmodel,andconsequentlyH0.Onethereforeneedsto
tem; ∆τ can be tightly constrained by the spatially ex- obtainagoodPSFmodelfortheAOdata,andtherearere-
tended images (we usually call them“arcs”) of the lensed centstudiesthataimtodosodirectlyfromtheAOimaging.
background galaxy (e.g., Kochanek et al. 2001; Suyuet al. Lagattuta et al. (2010) use three Gaussian components as
2009), together with stellar kinematics of the foreground thePSFmodeltosubtracttheAGNlightwhichissufficient
lens galaxy (e.g., Treu & Koopmans 2002; Koopmans et al. tostudy thelensing galaxy and its substructures.However,
2003; Suyuet al. 2010, 2014) and studies of the lens the analytical model is not sufficient to describe the com-
environment combined with ray-tracing through numeri- plexity of the PSF (see Figure 1 of Lagattuta et al. 2010)
cal simulations (e.g., Hilbert et al. 2007, 2009; Suyuet al. which could potentially impact the cosmographic measure-
2010;Fassnacht et al.2011;Greeneet al.2013;Collett et al. ments.Rusu et al. (2015)useeitherananalyticorahybrid
2013). The stellar kinematics and lens environment stud- PSFtostudythehostgalaxiesofthelensedAGNs(seealso
iesareimportantforovercomingthemass-sheetdegeneracy Rusuet al. 2014). The hybrid PSF is built from elliptical
and source-position transformations in lensing (Falco et al. Moffat profiles (Moffat 1969) with central parts iteratively
1985;Schneider& Sluse2013,2014;Xu et al.2015).There- tunedtomatchasingleAGNimage.WhilethishybridPSF
fore,bymeasuringthetimedelaysbetweenthemultipleim- is useful for extracting properties of the AGN host galaxy,
ages and modeling the lens and line-of-sight mass distribu- thecentralpartsofthePSFmodelcouldmanifest thenoise
tions, we can constrain D , which is the so-called time- pattern in the image (see Figure B.7 of Rusu et al. 2015)
∆t
delay distance that encompasses cosmological dependences which also could potentially impact cosmographic measure-
and is particularly sensitive to the Hubble constant (e.g., ments. Agnello et al. (2015) use an iterative method to re-
Suyuet al. 2010). The time delays in combination with the construct the PSF directly from lens imaging by averaging
stellar velocity dispersion measurements of the lens galaxy thedoublylensedAGN.Thismethodisvalidonlywhenthe
further allow us to infer the angular diameter distance to lensed AGN are separated far enough from each other. For
thelens galaxy (Paraficz & Hjorth 2009;Jee et al. 2014).
? have shown that for each lens system we can mea-
sure H0 to 7% precision. Hubble Space Telescope (HST) 1 Andnoequivalentopticalspace-basedtelescopemightbeforth-
∼
imaging is imperative for this analysis because it not comingsoon.
only provides high angular resolution but also a stable 2 COSmologicalMOnitoringofGRAvItational Lenses
MNRAS000,1–20(2015)
Strong lensing AO imaging with time delays for cosmography 3
typicalquad(four-image)lenssystems,thelensedAGNsare
′′
often close in separation (within 2 ), leading to overlaps in
thewingsoftheAGNimagesthataresmearedbythePSF.
Ourgoalistoprovideageneralmethodtoovercometheun-
known PSF model problem by extracting the PSF directly
from strong lensing imaging and simultaneously modeling
thelensmassdistribution.Wetestourmethodonsimulated
AO images, and apply the method to the known gravita-
tionallensRXJ1131 1231withKeckAOimaging,apartof
−
datafromSHARP3,whichisaprojectthatfocusesonstudy-
ing known quadruple-image and Einstein ring lenses using
high-resolution AO imaging, in order to probe their mass
distributions in unprecedented detail (e.g., Lagattuta et al.
2010,2012;Vegettiet al.2012,Hsuehetal.submitted).The
gravitationallenssystemRXJ1131 1231wasdiscoveredby
−
Sluseet al. (2003) who also measured the lens and source
redshifts to be 0.295 and 0.658, respectively. The HST ob-
servationsofthesystemRXJ1131 1231havebeenmodeled
−
by ?Suyuet al. (2014) for cosmography and more recently
byBirrer et al. (2015).
The outline of the paper is as follows. In Section 2, we
describe the observation of RXJ1131 1231 with the adap-
−
tive optics imaging system at the Keck Observatory. We
briefly recap in Section 3 the basics of cosmography with
time-delay lenses. In Section 4, we describe our new proce- Figure 1. Keck AO image (K′ band) of the gravitational lens
duretoanalyzeAOimageswithoutinformationonthePSF RXJ1131-1231.ThelensedAGNimageofthespiralsourcegalaxy
are marked by A, B, C and D, and the star-forming regions in
in advance. In Section 5, we use simulated data to test and
the background spiralgalaxy formplentiful lensed features. The
verify themethod.InSection 6,we demonstratetheresults
foregroundmainlensandthesatelliteareindicatedbyGandS,
fromrealdataandprovideacomparisonbetweentheresults
respectively.
fromHST imagingandAOimaging.Finally,wesummarize
in Section 7.
3 BASIC THEORY
3.1 The Theory of Gravitational Lensing
with Time Delay
2 OBSERVATION
Inthissectionwebrieflyexplaintherelationbetweengravi-
The RXJ1131 1231 system was observed on the nights of
− tationaltimedelaysandcosmology.Whenalightraypasses
UT2012May16andMay18withtheNearInfraredCamera
near a massive object, it experiences a deflection in its tra-
2(NIRC2)on theKeck-2Telescope (e.g., Wizinowich et al.
jectory and acquires a time delay by the gravitational field
2003).ThisimagewasapartofSHARPdata.Theadaptive
with respect to the travel time without the massive object.
optics corrections were achieved through the use of a R =
Therefore,thetimedelayhastwocontributions:(1)thegeo-
15.8tip-tiltstarlocated54.5arcsecondsfromthelenssystem
metricdelay,∆t ,whichiscausedbythebenttrajectory
andalaserguidestar.Thesystemwasobservedinthe“Wide geom
′′ ′′ beinglongerthan theunbentone,and (2) thegravitational
Camera” mode, which provides a roughly 40 40 field
of view and a pixel scale of 0.0397 arcseconds.×This pixel delay, ∆tgrav, which is due to the fact that the space and
time are affected around the gravitational field, so after in-
scaleslightlyundersamplesthepointspreadfunction(PSF),
tegrating the gravitational potential along the path, a far
but the angular extent of the lens system and the distance
away observer receives the light later by a Shapiro delay
from thetip-tilt starmadetheuseof theWideCamera the
(Shapiro1964; Refsdal 1964).
preferable approach.
The combination of thetwo delays is
The observations consisted of 61 exposures, each con-
sisting of 6 coadded 10 s exposures, for a total on-source
integration time of 3660 s. The data were reduced by a ∆t= D∆t 1(θ β)2 ψ(θ) , (1)
python-based pipeline that has steps that do the flat-field c 2 − −
(cid:20) (cid:21)
correction, subtract the sky, correct for the optical distor-
whereθ,β,andψ(θ)are theimage coordinates, thesource
tions in the raw images, and combine the calibrated data
coordinates, and the lens potential respectively. The time-
frames (for details, see Augeret al. 2011). The final image
delay distance is definedas
has a pixel-scale of 0.04 arcseconds and is shown in Figure
1. D (1+z )DdDs H−1, (2)
∆t ≡ d D ∝ 0
ds
whereD ,D andD aretheangulardiameterdistancesto
d s ds
3 Strong-lensingHighAngularResolutionProgram(Fassnachtet thelens,tothesource,andbetweenthelensandthesource,
al.inprep.) respectively. Thus, we can measure D∆t via gravitational
MNRAS000,1–20(2015)
4 G.C.-F. Chen et al.
lensing with time delays. Notice that the gradient of the (2005)toachieveefficientsamplingandtotestconvergence.
term in thesquarebracketsin Equation (1) vanishesat the The pragmatic procedure for convergence is described in
positions of thelensed images and yields thelens equation Suyu& Halkola (2010). We use Bayesian language in the
following sections.
β=θ ψ(θ), (3)
−∇
which governs thedeflection of light rays.
We refer the reader to, e.g., Schneideret al. (2006), 4 METHOD: PSF RECONSTRUCTION AND
Bartelmann (2010), Treu (2010), Suyuet al. (2010), LENS MODELING
Treu & Ellis (2014) for more details.
Inthissection,wedescribeanovelproceduretoanalyzethe
AO imaging without a PSF model a priori. Readers who
3.2 Probability Theory are not planning to use this method may wish to proceed
directly toSection 5 on thescientificresults enabledby the
Ameaningfulmeasurementshouldhaveanuncertaintyasa
method.
reference and it is also thekeyto confirm or ruleout possi-
The assumption of this strategy is that the PSF does
ble models. Thus, we need to analyze our data based on a
not change significantly within several arcseconds, which
probabilitytheorythatcanpresentthisidea.Bayes’theorem
is typically valid in AO imaging (vanDam et al. 2006;
provides the conditional probability distribution, so we can
Wizinowich et al. 2006). We show an overall flow chart in
obtain the posterior probability distribution of the model
Figure 2 to illustrate how to obtain iteratively the PSF,
parametersgiventhedatafrom Bayes’rule.Forexample,if
backgroundsourceintensity,thelensmass andlight model.
weareinterestedintheposterioroftheparametersπ ofthe
In Section 4.1, we decompose the observed light from
hypothesis model H given the data d, it can be expressed
the lens system into three components (lens galaxy, lensed
as
arcs of the background source galaxy, and the lensed back-
ground AGN) and introduce the notation that we will use
likelihood prior in the subsequent discussion. In Section 4.2, we obtain the
posterior
P(dπ,H)P(πH) preliminaryglobalstructureofAGNlightmodel,whilesep-
P(πd,H)= | | , (4)
| z }|P(d{Hz)}| { aratingthelenslightandarclight.InSection4.3,weobtain
z }| { | the fine structure of the AGN light and incorporate it into
evidence(marginalizedlikelihood)
the preliminary AGN light model. This is accomplished by
where the Bayesian evidence c|an{zbe} used to rank the correcting the PSF model. In Section 4.4, we update the
modeland ourpriorbased onthedata(e.g., MacKay 1992; PSFanduse it tomodel thelens mass and sourceintensity
Hobson et al. 2002; Marshall et al. 2002) distributions.Sincethelensgalaxylightisquitesmoothand
In addition, if we are interested in the posterior of a less sensitive to the PSF model, we use the PSF built from
specificparameter,πN,theposteriordistributioncanbeob- the AGN light for the lens galaxy light model. The PSF
tained by marginalizing overother parameters updatingandlensmassmodelingarerepeateduntilthecor-
rections to the PSF become insignificant. (See the criteria
N−1 in Section 4.3.3 and Section 4.4.3.)
P(π d,H)= P(πd,H) dπ . (5)
N i
| |
Z i=1
Y 4.1 Light components of the lens system
As shown in Figure 3, our model for the observed light in
3.3 Markov chain Monte Carlo
thelenssystemon theimage planehasthreecontributions:
the lens galaxy light, the arc light (the lensed background
Obtaining the probability distribution function of the pa-
source, i.e., the host galaxy of the AGN), and the light of
rameters in a model can benon-trivial,especially when the
themultipleAGNson theimage plane. We define
number of parameters is high. It is computationally unfea-
sible to explore a high-dimensional parameter space on a d=dP+n, (6)
regulargridsincethenumberofthegridpointsforthetask
exponentiallyincreaseswiththenumberofdimensions.Due wheredisthevectorof observeddata(image pixelvalues),
to the fact that the parameter space is typically large in lenslight arclight AGNlight
strong lensing analyses, one can bypass the use of grids by
dP = Kg + KLs + Mw , (7)
obtainingsamplesinthemulti-dimensionalparameterspace
thatrepresenttheprobabilitydistribution(i.e., thenumber and n is the noise inz}t|h{e datza}c|h{aractezri}z|e{d by the covari-
densityofthesamplesisproportionaltotheprobabilityden- ancematrixC (weusesubscriptDtoindicate“data”).The
D
sity). A Markov Chain Monte Carlo (MCMC) provides an blurring matrix K accounts for the PSF convolution, the
efficientwaytodrawsamplesfromtheposteriorprobability vector g is the image pixel values of the lens galaxy light,
density function (PDF) of the lens parameters, because of thematrix L maps source intensity to the image plane, the
theapproximatelylinearrelationbetweenthecomputational vectorsdescribesthesourcesurfacebrightnessonagridof
time and thedimension of theparameter space. pixels, the matrix M is composed using the positions and
WeuseMCMCsamplingthatisimplementedinGlee, theintensitiesoftheAGNs,andwisthevectorofpixelval-
a strong lens modeling software developed by S. H. Suyu ues of the PSF grid. We refer to Treu & Koopmans (2004)
and A. Halkola (Suyu& Halkola 2010; Suyuet al. 2012b). for constructing K and L, and illustrate the effect of M in
It is based on Bayes’ theorem and follows Dunkleyet al. Figure 3.
MNRAS000,1–20(2015)
Strong lensing AO imaging with time delays for cosmography 5
Figure 2. The flow-chart describes the overall procedures in Section 4. We use the procedures to reconstruct the PSF directly from
lensimageanddothelensmodeling.Instep1,weuseanearbystar (oroneofthe lensedAGNitself)as theinitialPSF;instep2,we
sequentiallyobtain thelens light,arclight,AGNs light,andthe positionsandrelativeamplitudes of AGNs;steps 3to5formaninner
looptoaddthecorrection(finestructures)intothePSFandaccumulatethecorrectionuncertainties;instep6,weentertheouterloop
which updates the image covariance matrix, PSF of all light model, and then repeat the full procedure until the image χ2 no longer
decreases.
At first, since we do not know the AO PSF a priori, thepredictedtotallightsequentiallyintolenslight,arclight,
K and w are just the initial blurring matrix and PSF grid and AGN light. We detail this process in Section 4.2.1 to
values,respectively.Asweiterativelymodelthelightcompo- Section 4.2.3 below.
nentsandcorrect thePSF,weupdatew (andsubsequently
K).
4.2.1 Lens Light Model (Step 2)
4.2 Determining the light components
Formodelingthelightdistributionofthelensgalaxy,weuse
The goal in this section is to obtain the preliminary model
parametrized profiles, such as theelliptical S´ersic profile,
ofeachof thethreelightcomponents.Instep1ofFigure2,
weinputtheobservedimageintothelensmodelingsoftware
I (θ ,θ )
GleewithanearbystarasourinitialPSFmodel.Ifthereis S 1 2
nonearbystar,anystarinthefieldcanbeusedastheinitial
PSForwecanuseoneoftheAGNimages.Adifferentinitial
PSFdoesnot affect thefinalresults, althoughwe notethat θ2+θ2/q2 1/ns´ersic
a good initial PSF would be helpful as they would require =I exp k 1 2 L 1 , (8)
feweriterationsofPSFcorrections.Instep2,wedecompose s −  p Reff ! − 
  
MNRAS000,1–20(2015)
6 G.C.-F. Chen et al.
Figure 3. Top panel: we decompose the image into lens light, arc light, and AGN light sequentially. Bottom panel: we model the
AGN lightbyplacingthe PSFgrid(describedby vector w)at each ofthe AGN positionsand scalingeach PSFbyitsrespective AGN
amplitude.ThisprocedurecanbecharacterizedbyamatrixM,suchthattheAGNlightmodelontheimageplanecanbeexpressedas
Mw.
where I is the amplitude, k is a constant such that R is isthenormalization fortheprobability.Thecovariancema-
s eff
theeffectiveradius,q istheminor-to-majoraxisratio,and trix,C ,istheoriginalcovariancematrixwithen-
L D,mArcAGN
n is the S´ersic index (S´ersic 1968). triescorrespondingtothearcandAGNmaskregionboosted
s´ersic
In order to get a preliminary model of the lens light, (see Appendix A), and N is the number of image pixels.
d
wemaskoutthearclight andAGNlightregion;thatis,we Wedenoteηˆ asthemaximumlikelihoodparameters(which
boost the uncertainty of the region where the arc light and maximizes equation 11).
theAGNlightareapparentlydominant.Thus,inthefitting Since the initial PSF is a prototype, usually there are
region, equation (7) becomes effectively somesignificantresidualsinthelenslightcenterwhenmax-
imizingtheposterioroflenslightparameters.However,this
dP =Kg. (9)
does not affect the subsequent lens light prediction in the
By Bayes’ rule,we have arc region, because the residuals are far from the arc re-
gions. To recap, we can obtain ηˆ by masking out the arc
P(ηd) P(dη)P(η), (10)
| ∝ | light and AGNlight regions.
where η represents the parameters of lens light (such as
I ,q ,n ,R ).Weassumeuniformprioronthelenslight
s L s´ersic eff
parameters, so we want to obtain 4.2.2 Arc Light Model (Step 2)
P(dη)= exp[−ED,mArcAGN(d|η)], (11) For modeling the arc light, we describe the source inten-
| ZD,mArcAGN sity on a grid of pixels on the source plane and map the
source intensity values onto the image plane using a lens
where,
mass model (via the operation KLs in equation (7)). We
ED,mArcAGN(dη) use elliptically symmetric power-law distributions to model
|
thedimensionless surface mass density of lens galaxies,
= 12(d−Kg)TC−D,1mArcAGN(d−Kg) 3 γ′ θ γ′−1
κpl(θ1,θ2)= 1−+q θ2+Eθ2/q2! , (14)
1 2
1
= 2χ2mArcAGN, (12) whereγ′ istheradialpower-lawpslope(γ′=2corresponding
to isothermal), θ is the Einstein radius, and q is the axis
and E
ratio of the elliptical isodensity contour. In addition to the
Z =(2π)Nd/2(det C )1/2 (13) lens galaxies, we includea constant externalshear with the
D,mArcAGN D,mArcAGN
MNRAS000,1–20(2015)
Strong lensing AO imaging with time delays for cosmography 7
Figure 4.Theleftpanel istheglobalstructureofthe PSF.Themiddlepanel isthecumulativefinestructureofthePSF.Weaddthe
finestructuretoglobalstructuretogetthePSFmodelintherightpanel.
following lens potentialin polar coordinates θ and ϕ:
ψext(θ,ϕ)= 12γextθ2cos2(ϕ−φext), (15) = 12(d−Kgˆ−KLs)TC−D,1mAGN(d−Kgˆ−KLs), (20)
whereγ isshearstrengthandφ istheshearangle.The
ext ext and
◦
shear position angle of φ =0 corresponds to a shearing
ext
along θ1 whereas φext = 90◦ corresponds to shearing along ZD,mAGN=(2π)Nd/2(det CD,mAGN)1/2 (21)
θ .4
2
We model the arc light with the lens light fixed. Since is the normalization for the probability. We discuss the
theAGNlightdominatesneartheAGNimagepositions,we “mAGN”regionsinAppendixA.InthesecondtermofEqua-
mask out the region where the arc is hard to be seen; that tion (19), θAGN,i is the measured AGN image position and
is, we want to minimize the contribution to the source in- σAGN,i is the estimated positional uncertainty of AGN im-
tensityreconstructionfromtheAGNlight.Sincetheregions agei;inthethirdterm,∆ti isthemeasuredtimedelaywith
oftheAGNaremaskedout,wetemporarily5 droptheAGN uncertaintyσ∆t,i for image pair i=AB,CB, or DB. After
component,Mw, in Equation (7),which given ηˆ becomes we maximize the likelihood of the data, we obtain ζˆ, and
alsothepredictedarclightofthereconstructedbackground
dP =Kgˆ+KLs, (16) source intensity, sˆ, of the AGN host galaxy. Note that if
where gˆ=g(ηˆ). The posterior based on thearc light is there is no time-delay information, one can remove the last
term in Equation (19).
P(ζd,∆t,ηˆ) P(d,∆tηˆ,ζ)P(ζ), (17)
| ∝ |
where ζ are the parameters of the lens mass distributions
(such as γ′, θE, γext). The likelihood of the data can be 4.2.3 AGN Light Model (Step 2)
expressed as
In Equation (7), we use Mw to represent the AGN light.
P(d,∆tηˆ,ζ)= dsP(d,∆tηˆ,ζ,s)P(s), (18) Inthenextsection, wefurtherdecomposethePSF,w,into
| |
Z the global structure and the fine structure that are shown
where in Figure 4. In particular, we define
P(d,∆tηˆ,ζ,s) w=w +T δw , (22)
| [0] [0] [0]
where w is the vector of global structure, δw is the
[0] [0]
exp[ E (dηˆ,ζ,s)]
= − D,mAGN | vector of fine structure and the subscript, [0], represents
ZD,mAGN the zero-th iteration. Since, in this section, we focus on the
global structure of the PSF, we postpone the discussion of
NAGN 1 exp |θAGN,i−θpAGN,i|2 T toEquation (29) and let
· i=1 √2πσAGN,i − 2σA2GN,i ! w=w[0]. (23)
Y
1 (∆t ∆tp)2 By using ηˆ, ζˆ, and sˆ from the previous two sections and
exp i− i , (19) keeping them fixed, we model the global structure of the
·i=1√2πσ∆t,i "− 2σ∆2t,i # PSFwith Gaussian profiles, each of theform
Y
θ2+(θ2/q2)
ED,mAGN(d|ηˆ,ζ,s) IG(θ1,θ2)=Igexp − 1 2σ22 g , (24)
(cid:20) g (cid:21)
4 Our(right-handed)coordinatesystem(θ1,θ2)hasθ1alongthe whereIg istheamplitude,qg istheaxisratio,andσg isthe
East-Westdirectionandθ2 alongtheNorth-Southdirection. width.We findthat a few Gaussians (∼2−4) with a com-
5 WewillputtheAGNcomponentbackinnextsection,4.2.3 moncentroidaresufficientindescribingtheglobalstructure
MNRAS000,1–20(2015)
8 G.C.-F. Chen et al.
Figure5. ThePSFcorrectiongridsoftheiterativePSFcorrectionscheme.IntheinnerloopofthePSFcorrectionscheme(samenbut
differentm),westartwithasmallcorrectiongridδw andincreaseitsequentially.Thisaccommodates forthelargercorrectionsneeded
in the central parts of the AGN. Left panel: a small PSF correction grid is placed at each of the four AGN images A, B, C and D in
Figure 1 via the matrix M, and the values of the PSF correction grid is determined via a linear inversion to reduce the overall image
residuals. Since the AGN centroids are typically non-integral pixel values, we linearly interpolate the correction grid onto the image
plane.Rightpanel:theenlargedcorrectionsgridsafterseveraliterationsofPSFcorrection,showingoverlapbetweenthegrids.Whenthe
peripheralareaofacorrectiongridoverlapswiththecentralpartsofanotherAGNimage(e.g.,AGNimageCinthelower-rightpartsof
the correctiongridofimageA), wemaskout thecenter of theAGN regioninorder toprevent thecorrection gridsfromabsorbingthe
residualswhichcomefromthemismatchofthesharpintensityofAGNcenter (seeAppendixAformoredetails).
of the PSF.6 SubstitutingEquation (23) into Equation (7), and ZD = (2π)Nd/2(det CD)1/2. We denote νˆ and ξˆ as the
given ηˆ, ζˆ and sˆ, we obtain maximum likelihood parameters (that maximizes Equation
(26)) from which we can obtain the optimal AGN light on
dP =Kgˆ+KLˆsˆ+Mw , (25)
[0] theimage,giventheoptimizedsourceandlensmassmodels
where Lˆ =L(ζˆ), which is kept fixed at this step. Note that from theprevioussections.
theKmatrixhereisbasedontheinitialPSFmodel,before
the multi-Gaussian fitting. The posterior of the PSF and
AGN parameters is given by 4.3 Pixelated Fine Structure of AGN light
P(dηˆ,ζˆ,ν,ξ)P(ν,ξ) Inthissection,weintroducetheinnerloopwhichaimsatex-
P(ν,ξ|d,ηˆ,ζˆ)= | P(dηˆ,ζˆ) , (26) tractingthefinestructure,δw[0],inEquation (22) byusing
| a correction grid. Weshow it visually in Figure 5. The goal
whereνrepresentstheparametersoftheGaussianprofilesin of the inner loop is to incorporate most of the fine struc-
Equation(24)thatyield w[0];ξ aretheamplitudesandthe tures into the PSF model; then in the outer loop, we can
positionsoftheAGN,whicharecodedinM.Thelikelihood use the updated PSF model obtained from the inner loop
of Equation (26) is to remodel all the light components (which require a given
PSF model). Since this section is the starting point of the
P(dηˆ,ζˆ,ν,ξ)= exp[−ED(d|ηˆ,ζˆ,ν,ξ)], (27) inner loop and outer loop, the ηˆ, ζˆ, sˆ, νˆ, and ξˆ we get by
| ZD optimizing Equations (10), (17), and (26) are actually the
where zero-th outer loop iteration and the zero-th inner loop iter-
ation, which we denote byηˆ , ζˆ , sˆ , νˆ ,and ξˆ .
E (dηˆ,ζˆ,ν,ξ) [0] [0] [0] [0,0] [0,0]
D
|
1 4.3.1 PSF Correction for Each Iteration
= 2(d−Kgˆ−KLˆsˆ−Mw[0])T (Inner Loop: Step 3)
In general, given ηˆ , ζˆ , sˆ , νˆ 7, and ξˆ , where
[n] [n] [n] [n,m] [n,m]
C−1(d Kgˆ KLˆsˆ Mw ), (28) m is the iteration number of the inner loop and n is the
· D − − − [0]
6 ThedifferentGaussiancomponents canvarytheiramplitudes, 7 νˆ is onlypresent whenn=m=0, which corresponds to
[n,m]
positionanglesandaxisratios. parametersoftheGaussianprofilesinEquation(24)
MNRAS000,1–20(2015)
Strong lensing AO imaging with time delays for cosmography 9
iteration number of the outer loop, we can write out the determined yet), and N is the number of pixels of
δw,[n,m]
equation as thecorrection grid.Weusethecurvatureform forthefunc-
tion E , which is discussed in Suyuet al. (2006).
dP =K gˆ +K Lˆ sˆ δw,[n,m]
[n] [n] [n] [n] [n] Again, it is easy to understand that we want to maxi-
mize Equation (30). We obtain the most probable solution
+Mˆ (wˆ +T δw ) dP , (29)
[n,m] [n,m] [n,m] [n,m] ≡ correction δw =(F+λ H)−1(M T )TC−1 u,
where [n,m] δw,[n,m] [n,m] [n,m] D,mAc
(34)
wˆ = w[n,m](νˆ[n,m]) if n=m=0 where
[n,m] w otherwise
(cid:26) [n,m] F= E
gˆ = g(ηˆ ), Lˆ = L(ζˆ ), Mˆ = M(ξˆ ), K ∇∇ D,mAc,[n,m]
[n] [n] [n] [n] [n,m] [n,m] [n]
is the nth blurring matrix (we explain how to get the nth
blurring matrix in Section 4.4.1), T[n,m] is a matrix which =TT (MT C−1 M )T , (35)
makes δw the same length as wˆ by padding with [n,m] [n,m] D,mAc [n,m] [n,m]
[n,m] [n,m]
zeros(weshowitvisuallyinAppendixB),andδw isthe
[n,m]
finestructurewe want to obtain bytheend of this section.
The posterior of δw[n,m] is H=∇∇Eδw,[n,m], (36)
P(δw d,ηˆ ,ζˆ ,νˆ ,ξˆ ,λ ,R)
[n,m]| [n] [n] [n,m] [n,m] δw,[n,m]
u=d K gˆ K Lˆ s Mˆ wˆ , (37)
P(dδw ,ηˆ ,ζˆ ,νˆ ,ξˆ ) − [n] [n]− [n] [n] [n]− [n,m] [n,m]
= | [n,m] [n] [n] [n,m] [n,m] and
P(dλ ,ηˆ ,ζˆ ,νˆ ,ξˆ ,R)
| δw,[n,m] [n] [n] [n,m] [n,m] ∂
. (38)
∇≡ ∂δw
P(δw λ ,R), (30) [n,m]
· [n,m]| δw,[n,m]
Now,wegobacktofindtheoptimalregularizationcon-
where P(δw λ ,R) is the prior on δw given
[n,m]| δw,[n,m] [n,m] stant;that is, we want to maximize
λ ,R withRdenotingaparticularformof“regular-
{ δw,[n,m] }
ization”onδw andλ characterisingthestrength P(dR,λ )P(λ )
[n,m] δw,[n,m] P(λ d,R)= | δw[,[n,m] δw,[n,m] (39)
of the regularization. We can write the likelihood in Equa- δw,[n,m]| P(dR)
tion (30) as |
using Bayes’ rule. If we assume a flat prior in log λ ,
P(d|δw[n,m],ηˆ[n],ζˆ[n],νˆ[n,m],ξˆ[n,m]) we want to maximize P(d|R,λδw,[n,m]), whichδwis,[n,tmh]e
evidence in equation (30). Following the results from
Suyuet al. (2006),we get
exp[ E (dδw ,ηˆ ,ζˆ ,νˆ ,ξˆ )]
= − D,mAc,[n,m] | [n,m [n] [n] [n,m] [n,m] , 2λˆ E (δw )
Z δw,[n,m] δw,[n,m] [n,m]
D,mAc
(31)
where“mAc”standsfor maskAGNcenter,
ED,mAc,[n,m](d|δw[n,m],ηˆ[n],ζˆ[n],νˆ[n,m],ξˆ[n,m]) =Nδw,[n,m]−λˆδw,[n,m]Tr[(F+λˆδw,[n,m]H)−1H], (40)
whereTrdenotesthetraceandλˆ istheoptimalreg-
δw,[n,m]
= 21(d−dPcorrection)TC−D,1mAc(d−dPcorrection), (32) fiunlaerisztartuiocntucroen),stwaento.bIftawinesδewt[mn,0=].0D(uzeertooththiteersahtairopnionfttehne-
sityoftheAGNcenter,theresidualstherearemuchstronger
andZD,mAc=(2π)Nd/2(detCD,mAc)1/2 isthenormalization than theperipheral area. If we directly extract the full cor-
for the probability. We discuss the mAc (maskAGNcenter) rection grid, the regularization intends to under-regularize
regions in AppendixA on the peripheral area and over-regularize on the center.
The prior/regularization we impose in Equation (30) To avoid this problem, at first, we extract the correction
on the correction grid (fine structure of PSF) is to prevent only around the AGN center; that is, we start from small
thecorrectiongridfromabsorbingthenoiseintheobserved N (half lightradiusorsmaller) andincrease itgrad-
δw,[n,m]
image. Weexpress theprior in thefollowing form ually(around1.2timesprevioussizeeachtime)asweobtain
P(δw λ ,R) δw[n,m].WeshowtheideainFigure5(notethattheindices
[n,m]| δw,[n,m] on δw in the figure are labeling the pixels, rather than the
iteration numbers).
= exp(−λδw,[m]Eδw,[n,m](δw[n,m]|R)), (33) Since every iteration of δw[n,m] has their own fine
Z (λ ) structure (correction) uncertainty, according to Suyuet al.
δw,[n,m] δw,[n,m]
(2006), we also take as estimates of the 1σ uncertainty on
where λ is the regularization constant of cor-
δw,[n,m] eachpixelvaluethesquareroot ofthecorrespondingdiago-
rection, Zδw,[n,m](λδw,[n,m]) = dNδw,[n,m]δw[n,m] exp nalelement of thecovariance matrix given by
( λ E )isthenormalizationofthepriorprob-
a−biliδtwy,[dn,imst]ribδuwt,[ino,nm](note that theRoptimal λ is not C =(F+λˆ H)−1. (41)
δw,[n,m] δw,[n,m] δw,[n,m]
MNRAS000,1–20(2015)
10 G.C.-F. Chen et al.
4.3.2 Add Fine Structure into Global Structure and
(Inner Loop: Step 4)
ξˆ ξˆ ξˆ . (49)
[n,Ninner] ≡ [n+1,0] ≡ [n+1]
We start with the zeroth inner loop iteration, by setting
m = 0, of the global structure, w , and fine structure,
[n,0] 4.4 Lens Modeling with updated PSF
δw (which we can obtain by following the previous two
[n,0]
sections). We then add the fine structure into the global The goal of the outer loop in Figure 2 is to remodel all the
structurebydefining light components with the updated PSF; that is, we want
toobtainabetterlenslightmodelandarclightmodelwith
w =w +T δw , (42)
[n,1] [n,0] [n,0] [n,0] the new blurring matrix, and the underlying fine structure
wherew is thefirstiteration in innerloop. More gener- can then berevealed.
[n,1]
ally, we definethem+1th iteration of the PSFas
w =w +T δw . (43) 4.4.1 Update the Blurring Matrix and the Image
[n,m+1] [n,m] [n,m] [n,m]
Covariance Matrix (Outer Loop: Step 6)
We recalculate the AGN parameters every time after
getting a new w , so given the same ηˆ and ζˆ in Blurringmatrix:AfterobtainingthelastversionofthePSF
[n,m+1] [n] [n]
Equation(29),theposterioroftheAGNparametersisgiven from Section 4.3.3, we update the blurring matrix, K, in
by Equation (7). In order to accelerate modeling speed, which
highlydependsonthesizeofthePSFforconvolutionofthe
P(ξ d,ηˆ ,ζˆ )
[n,m+1]| [n] [n] extendedimages,wechoosethecentrall[n]×l[n]pixelsofthe
updatedPSFgrid(thathasN pixels)asthenew
δw,[n,Ninner]
PSFtoconstructK forthespatiallyextendedimages.8
= P(d|ηˆ[n],ζˆ[n],ξ[n,m+1])P(ξ[n,m+1]). (44) Image covarian[cne+1m] atrix : We accumulate the uncer-
P(dηˆ ,ζˆ ) tainty of the PSF pixel grid from every inner loop. The ac-
| [n] [n]
cumulated uncertaintyis
(recall that ξ represents the relative amplitudes and
the positions[no,fmt+h1e] AGNs in the nth outer loop iteration, n2 =Ninner T C δ , (50)
andm+1thinnerloopiteration).ThelikelihoodofEquation δw,[n+1],k [n,m],ki δw,[n,m],ij ij
m=0 i
(44) is X X
where T is the element at k row and i column of
[n,m],ki
T , C is the element at i row and j column of
[n,m] δw,[n,m],ij
P(d|ηˆ[n],ζˆ[n],ξ[n,m+1]) nC+δw,1[nth,mn],oiasnedvδeicjtoisrtihsedKefirnoendecakserdelta.Theelementofthe
= exp[−ED,[n,m+1](dZ|ηˆ[n],ζˆ[n],ξ[n,m+1])], (45) n[n+1],µ=sn2µ+Xk Mˆ[n+1],µkn2δw,[n+1],k, (51)
D
which is characterized by the covariance matrix C 9.
where D,[n+1]
Notethatn istheelementoftheoriginaldatanoisevector.
µ
E (dηˆ ,ζˆ ,ξ )
D,[n,m+1] | [n] [n] [n,m+1]
4.4.2 Lens Modeling with All Light Components
= 1(d Ω)TC−1(d Ω) (46) (Outer loop: Step 2)
2 − D − In general, when executingthenextiteration of outerloop,
with we can express Equation (7) as
Ω=K[n]gˆ[n]+K[n]Lˆ[n]sˆ[n]+M[n,m+1]w[n,m+1], (47) dP =K[n+1]g[n+1]+K[n+1]L[n+1]s[n+1]+M[n+1]w[n+1]
and ZD, as usual, is (2π)Nd/2(det CD)1/2. After maximiz-
ing Equation (44), we obtain ξˆ . We then replace the dP . (52)
[n,m+1] ≡ total
ξˆ from theprevious iteration with the ξˆ we just
[n,m] [n,m+1] Theposterior can be written as
obtained,and conductthenext innerloop iteration.
P(η ,ζ ,ξ d,∆t)
[n+1] [n+1] [n+1]|
4.3.3 The Criteria to Stop the Inner Loop.
P(d,∆tη ,ζ ,ξ )P(η ,ζ ,ξ ).
During every inner loop, we gradually increase the size, ∝ | [n+1] [n+1] [n+1] [n+1] [n+1] [n+1]
(53)
N , of the correction grid. Then, if (1) there is no
δw,[n,m]
residuals outside the correction grid, (2) Equation (34) has
no intensity, and (3) Equation (46) no longer decreases, we 8 Forexample,wechoosel =nl
[n] [1]
stop the inner loop. Assuming we have Ninner iterations in 9 Thepurposeofupdatingthetheimagecovariancematrixisto
theinnerloop,weobtainw andξˆ .Wethen speedupthemodelingtothefinalanswersincethecorrectionun-
define [n,Ninner] [n,Ninner] certaintythatweaddintotheimagecovariancematrixisaround
AGN;thatis,weweightthearcregionmore.However,intheend,
w[n,Ninner] ≡w[n+1,0] ≡w[n+1] (48) ifthereisno“correction”,Equation(50)isclosetozero.
MNRAS000,1–20(2015)