Interstellar Holography

0.36 MB·

by  Mark Walker

Mon.Not.R.Astron.Soc.000,000–000(0000) Printed3February2008 (MNLATEXstylefilev2.2)
Interstellar Holography
M.A. Walker1,2,3,4, L.V.E. Koopmans3, D.R. Stinebring5,6, W. van Straten4,7,8
1.ManlyAstrophysicsWorkshopPtyLtd,Unit3,22CliffStreet,Manly,NSW2095,Australia
2.SchoolofPhysics,UniversityofSydney,NSW2006,Australia
3.KapteynAstronomicalInstitute,UniversityofGroningen,P.O.Box800,9700AVGroningen,TheNetherlands
4.NetherlandsFoundationforResearchinAstronomy,P.O.Box2,7990AADwingeloo,TheNetherlands
5.OberlinCollege,DepartmentofPhysicsandAstronomy,Oberlin,OH44074,U.S.A.
6.LeidenUniversityObservatory,Leiden,TheNetherlands
7.CenterforGravitationalWaveAstronomy,UniversityofTexas,Brownsville,TX78520,U.S.A.
8.SwinburneUniversityofTechnology,CentreforAstrophysicsandSupercomputing,Hawthorn,VIC3122,Australia
8
0
0 3February2008
2

n
a ABSTRACT
J Thedynamicspectrumofaradiopulsarisanin-linedigitalhologramoftheionisedinterstellar

8 medium.Ithaspreviouslybeendemonstratedthatsuchhologramspermitimagereconstruc-
2 tion,inthesensethatonecandetermineanapproximationtothecomplexelectricfieldvalues
as a function of Doppler-shift and delay, but to date the quality of the reconstructions has

] beenpoor.Herewereportasubstantialimprovementinthemethodwhichwehaveachieved
h
bysimultaneousoptimisationofthethousandsofcoefficientsthatdescribetheelectricfield.
p
ForourtestspectrumofPSRB0834+06wefindthatthemodelprovidesanaccuraterepresen-
-
o tationofthedataoverthefull63dBdynamicrangeoftheobservations:residualdifferences
r betweenmodelanddataarenoise-like.Theadventofinterstellarholographyenablesdetailed
t
s quantitativeinvestigationoftheinterstellarradio-wavepropagationpathsforagivenpulsarat
a
eachepochofobservation;weillustratethisusingourtestdatawhichshowthescatteringma-
[
terialtobestructuredandhighlyanisotropic.Thetemporalresponseofthemediumexhibits

1 ascatteringtailouttobeyond100µsandapulsearrivaltimemeasurementatthisfrequency
v and this epoch of observation would be affected by a mean delay of 15µs due to multipath
3 propagationintheinterstellarmedium.
8
1 Keywords: techniques:interferometric—pulsars:general—pulsars:individual:B0834+06
4 —scattering—ISM:structure—turbulence
.
1
0
8
0
1 INTRODUCTION facthasbeenemphasisedbyconsiderationofthetwo-dimensional
:
v powerspectraofthedynamicspectra,whereinpowerisoftenseen
i With modern instrumentation pulsar dynamic spectra can be tobeconcentratedinparabolicarcsandinvertedarclets(Stinebring
X
recordedathighspectralandtemporalresolutionyieldingadataset etal2001).Thereisaconsensusthatthisphenomenonarisesdi-
r withalargeinformationcontentfromjustoneobservation.Forex- rectlyfromthegeometryofthescatteringprocess(Stinebringetal
a
ample:ifweobserveapulsarforonehour,samplingthespectrum 2001;Cordesetal2006;Walkeretal2004),withwavesscattered
with1kHzchannelsevery10seconds,overatotalbandwidthof throughanangleθ(cid:126)experiencingaDoppler-shiftproportionaltoone
100MHz,thenwehave∼4×107independentfluxmeasurements. componentofθ(cid:126)andadelayproportionaltoθ(cid:126)·θ(cid:126).Incaseswherethe
Andifthepulsarisbrightandthetelescopeislargetheneachof parabolae are very sharp it has been argued that the scattering is
these measurements can have signal-to-noise ratio of order unity highlyanisotropic(Walkeretal2004;Cordesetal2006).Sharply
implyingatotalinformationcontentofpotentially∼40Mbit.This definedarcs/arcletsalsorequirethatthescatteringarisesinaregion
informationrelatesprimarilytomultipathscatteringoftheradio- ofsmallextentalongtheline-of-sight(Stinebringetal2001),soit
wavesbytheionisedInterstellarMedium(ISM),asitisthisscat- isnotdistributedturbulencebutdiscrete,localisedstructureswhich
teringwhichgivesrisetotheobservedinterferencefringes.Itmay areresponsibleforthesefeatures.However,thephysicalnatureof
bethatthesepathsaredeterminedbyrandomirregularities–e.g. the scattering media remains obscure. Consequently there is now
causedbyturbulence–intheISM.Inthatcaseanygivendynamic someincentivetoexploretheinformationcontentinindividualdy-
spectrum contains information about those random irregularities, namicspectra,inordertobuildadetailedpictureofthescattering
andthereisnostrongmotivationtostudytheindividualspectrain structures.Furthermotivationforinvestigatingindividualdynamic
detailastheyreflectparticularrealisationsofastochasticprocess. spectraisprovidedbystudiesofthepulsarsthemselves:iftheprop-
But some pulsar dynamic spectra exhibit a high level of order in ertiesofthescatteringscreenareaccuratelyknownitispossibleto
their fringe patterns — see Rickett (1991) for an overview. This
(cid:13)c 0000RAS
2 Walker,Koopmans,Stinebring&vanStraten
usethescreenforveryhighresolutionimagingofthepulsarmag- waveinterferencephenomenon,andbysimultaneousoptimisation
netosphere(WolszczanandCordes1987;Gwinnetal1997;Walker ofthehundredsofparameterswhichdescribetheintrinsictemporal
and Stinebring 2005, WS05 henceforth). Precision pulsar timing fluxvariationsofthepulsar.
programs(e.g.Manchester2007)alsoprovideanincentiveforun- Thispaperisorganisedasfollows:inthenextsectionwedetail
derstandingtheparticularinterstellarpropagationpathswhichcon- the improvements we have made to the reconstruction algorithm
tributetoindividualobservations—ifthepropagationdelaysre- describedbyWS05;in§3wepresenttheresultswehaveobtained,
mainuncorrectedinthedatatheyconstituteapotentiallylargesys- using the same test data employed by WS05; and in §4 we con-
tematicerrorinpulsearrivaltimemeasurements. siderwhatourtestdatatellusabouttheinterstellarmedium,with
It has previously been demonstrated that one can iteratively emphasisonthetemporalresponseofthescatteringmedium.
construct a model of the electric field as a function of radio-
frequencyandtime,startingfromtheobserveddynamicspectrum
(WS05).Theprocedureforachievingthisislargelyequivalentto
determining a phase for the electric field in each pixel of the dy- 2 OPTIMISATIONOFTHEE-FIELDMODEL
namic spectrum, because a noisy estimate of the field amplitude
InWS05wedescribedanalgorithmformodellingtheelectricfield
is given directly by forming the square root of the observed in-
structure in the Fourier domain conjugate to the dynamic spec-
tensity. This situation is common to the broad class of problems
trum. The latter is recorded as a function of radio-frequency, ν,
knownas“phase-retrievalproblems”,whicharewellknowninthe
andtime,t,andthecorrespondingconjugatevariablesaredelay,τ,
opticsliterature(e.g.Fienup1982).However,themethodofsolu-
andDoppler-shift,ω.TheWS05algorithmproceedsfromagridof
tion demonstrated for pulsar dynamic spectra appears to be new;
noisymeasurementsoftheelectricfieldenvelope,|U(ν,t)|2,toa
it exploits sparseness of the power distribution in the Fourier do-
main,andasolutionisobtainediterativelybyaddingdiscretenew modelofUe(τ,ω)intermsofdiscretewavecomponents,j:
fieldcomponentsinsuchawayastominimisethedifferencesbe- X
tweenthemodeldynamicspectrumandthedata.Conceptuallythe Ue(τ,ω)= u˜jδ(τ −τj)δ(ω−ωj), (1)
processhasastrongsimilaritytotheCLEANalgorithm(Ho¨gbom j
1974)whichiscommonlyusedinradioastronomicalimagingfor andthecomponents(characterisedbyτ ,ω andthecomplexnum-
j j
deconvolvingthesynthesisedtelescopebeamfroma“dirty”image ber u˜ ) are chosen one-by-one so as to minimise the differences
j
ofthesky;CLEANusuallyworkswelliftheimageisonlysparsely betweenmodelanddata.With∼ 9000componentsthemodelre-
populated with emission. The connection between the two algo- ported by WS05 gives a good visual match to the data, but the
rithmsisreinforcedwhenwerecognisethatdeterminingtheelec- residuals are large in comparison with the noise level, implying
tricfieldfromthedynamicspectrumisalsoequivalenttoadecon- large systematic errors. An inspection of the differences between
volution.Thedynamicspectrumasafunctionofradio-frequency the observed secondary spectrum (i.e. the power spectrum of the
andtimeissimplyI(ν,t)=U∗(ν,t)U(ν,t),whereU istheelec- dynamic spectrum) and the model secondary spectrum – the two
tric field; in the Fourier domain this relationship becomes a con- quantitiesshowninthelowerpaneloffigure1inWS05–reveals
volution Ie = Ue∗ ⊗Ue and so we are deconvolving the Fourier thatthediscrepanciesoccurpredominantlyatthesamelocationsin
transformoftheelectricfieldfromitscomplexconjugate. delay-Dopplerspacewherethereisalreadypowerpresentinboth
Some comments about nomenclature are appropriate at this modelanddata.Thissuggeststhatthesystematicerrorsintroduced
point.Anyprocesswhichallowsustoreconstructtheelectricfield bytheWS05algorithmarenotduetoincorrectlyplacedcompo-
which gave rise to a recorded fringe pattern is sensibly termed nents(i.e.wrongτ ,ω ),oraninsufficientnumberofcomponents
j j
“holography”, and the recorded fringe pattern (i.e. the dynamic inthemodel,butratherduetoerrorsindeterminingthevariousu˜ .
j
spectruminourcase)isthehologram.Inthecaseconsideredhere ItwasnotedinWS05thatglobaloptimisationofthe{u˜ }is
j
itis“digitalholography”becausethereconstructionisdoneinsoft- desirable,inordertoreducesystematicerrorsinthemodel,butdif-
ware,and“in-line”becausetheobject(i.e.thescatteringscreenin ficult to achieve because of the large number of free parameters
ourcase)istransparentandsitsinthebeamwhichformstherefer- whichwouldbeinvolvedinsuchanoptimisation.Inparticularin-
encewave.In-lineholographyissometimescalled“Gaborhologra- versionofa104 ×104 non-sparsematrix–whichisperhapsthe
phy”becauseitisthearrangementwhichwasoriginallyconceived most obvious approach to solution of the linearised least-squares
bytheinventorofholography,DennisGabor. optimisation–iscomputationallychallenging.Furthermoreanyap-
Todatethefidelityofelectricfieldreconstructionsfrompulsar proachbasedonsolutionofthelinearisedproblemwouldrequire
dynamicspectrahasbeenpoor.AlthoughthemodelofWS05suc- iterationinordertosolvethefullnon-linearoptimisationproblem.
cessfullyreproducesthegeneralappearanceoftheirtestdynamic Fortunatelythereareeasiermethods–see,forexample,Nocedal
spectrumtherearesubstantialquantitativedifferencesbetweenthe andWright(1999)–andwehaveusedoneofthesetooptimisethe
model and the data. These errors are seen when the model “sec- WS05electricfieldmodelaswenowdescribe.
ondary spectrum” (i.e the power spectrum of the dynamic spec- Themethodwhichweemployedisaquasi-Newtonmethod,
trum)iscomparedtothedata:thedataexhibitadynamicrangeof in which a demerit function S is minimised by seeking succes-
63dB,butthedifferencebetweentheWS05modelsecondaryspec- sivelycloserapproximationstothesolutionof∂S/∂x =0forall
i
trumandthesecondaryspectrumofthedatahasadynamicrange parametersx overwhichwewishtooptimise.Newton’smethod
i
of 47 dB, implying that the reconstruction has correctly captured requiresknowledge ofthe Hessian(i.e.all thesecond derivatives
onlythetop16dBofthedata.Herewereportmodificationstothe ∂2S/∂x ∂x ), which is computationally expensive when a large
i j
reconstructionprocesswhichhaveyieldeddramaticimprovements number,N,ofparametersisinvolved.Bycontrast,quasi-Newton
totheaccuracyoftheelectricfieldmodel;onthesametestdataset methodsproceedbyapproximatingtheHessian;informationfrom
asusedbyWS05wefindthatourimprovedmodelcapturesthefull current and previous iterations is used to update the approxi-
dynamicrangeofthedata.Thesegainswereachievedbysimulta- mateHessianforsubsequentiterationsyielding,ineffect,afinite-
neousoptimisationofthethousandsofparametersdescribingthe difference representation of the local curvature of S. The update
(cid:13)c 0000RAS,MNRAS000,000–000
InterstellarHolography 3
schemewhichweusedistheBFGS(Broyden-Fletcher-Goldfarb- (ii) Thenumberofnewelectricfieldcomponentspickedateach
Shanno) update. Specifically: we used the L-BFGS-B algorithm iterationisgivenbytheintegerpartof1+N /100,whereN isthe
c c
(Byrdetal1995),whichisa“Limitedmemory”implementationof currentnumberofelectricfieldcomponents.Soinitiallyonlyone
BFGS, using line-search minimisation, in which Bounds are per- new component is picked at each iteration, and when the model
mitted on the parameters x . The term “Limited memory” here containsalargenumberofcomponentsthefractionalincreaseper
i
referstothefactthattheN ×N Hessianisreplacedbytheouter iterationis1%.
productoftwom×N matrices,withmbeingthenumberofprior (iii) Whenthenumberofelectricfieldcomponentsexceeds100,
iterationswhichareemployedinconstructingtheupdate.Because thealgorithmisfreetopickcomponentswhichhaveτ <0.
m∼afew,andmdoesnotgrowwithN,thestoragerequirements (iv) Thealgorithmterminateswhenthereducedχ2 (i.e.χ2 per
of the algorithm are only modest and grow linearly with N. The degreeoffreedom)reachesunity,orwhenthereducedχ2 reaches
L-BFGS-Bcodewaswrittenbyspecialistsinthefieldofnumerical aminimum—whicheveroccursfirst.
optimisation; it is freely available as a set of FORTRAN subrou-
Because of the high dynamic range of the data it is important to
tines1.Inordertomakeuseofthiscodeitisnecessaryfortheuser
maintainhighprecisionintheoptimisation,sowesettheL-BFGS-
tosupplyroutineswhichevaluatethedemeritfunction,S,andits
Bparameternamed“factr”(whichmeasurestheprecisioninunits
partialderivatives,∂S/∂x ,withrespecttoalltheparameters,x ,
i i ofmachineprecision)tothevalue10.Thenumberofprevioussteps
overwhichwewishtooptimise.GiventheseinputstheL-BFGS-B
usedinformingeachBFGSupdateism=10.
code will search for a minimum in S. If a minimum is found by
L-BFGS-B it is not guaranteed to be a global minimum. The L-
BFGS-Bpackageiswelldocumented;thereareclearinstructions
3 RESULTS
onhowtousethecodeandsimpledriversareincludedinthepack-
age so that the user can verify correct performance on their own AdirectcomparisonoftheresultsofWS05withthenewalgorithm
computer. For our application we used a model with two sets of ispossiblebyusingthesametestdataasWS05employed.Those
parameters:thevariousu˜j,eachofwhichisdescribedbytwoun- dataareshowninfigure1,alongwiththemodelgeneratedbythe
boundedrealvariables,representingtherealandimaginaryparts, new algorithm; both model and data are shown in the form of a
andasetofpositive-definiterealnumbersdescribingtheintrinsic dynamic spectrum and its power spectrum (the “secondary spec-
pulsarfluxasafunctionoftime,{fk}.Theinclusionofthevarious trum”).Themodelwasgeneratedusingthealgorithmdescribedin
parametersfk : fk (cid:62)0 ∀ kdemandsthattheoptimisationsoft- §2.Figure1alsoshowstheresidualsbetweendataandmodeldy-
warebeabletohandlevariableswhicharebounded,asisthecase namicspectra,andthepowerspectrumofthoseresiduals;theresid-
fortheL-BFGS-Bpackage.TheWS05algorithmdoesnotattempt ualsappearnoise-likeinbothpanels.Somequantitativemeasures
tosolveforthefk explicitlybutsimplyappliesaFourier-domain ofthesuccessofthenewalgorithmareappropriate.Thenewalgo-
filter to the data in an attempt to remove the intrinsic pulsar flux rithmachievesareducedchi-squaredvalueofχ2 =1.00(thiswas
r
variations.Becausethereisnoclear-cutdistinctionbetweenintrin- thestoppingcriterionwhichwasreachedfirst),versusχ2 = 1.19
r
sicfluxvariationsandthoseduetowaveinterference,theprocedure achieved by WS05, and it does so with only 8,000 electric field
usedbyWS05isquitecrudeandinthepresentworkweuseitonly components versus 8,720 in WS05. Note, however, that the new
asastartingpoint(asperitem(i)in§3.1ofWS05). algorithm does employ an additional 270 real numbers – one for
OurfirstattemptatimprovingontheWS05approachwasto each time sample – to describe the intrinsic pulsar flux variation
take the output of the WS05 algorithm and use it as the starting with time; these numbers were in effect fixed in the WS05 algo-
pointofanoptimisationwithL-BFGS-B.Theresultsweregood, rithm by a Fourier-plane filter acting on the input data. Subtract-
yieldingamodelwhichcapturedmuchmoreofthedynamicrange ingthemodeldynamicspectrumfromthedata,andthenforming
inthedata,andhadalowervalueofthereducedχ2statistic.This the power spectrum of the residuals gives a sensitive test of the
resultwasencouraging,butitwasclearthatwecoulddobetter:the fidelity of the model because it picks out correlated errors in the
componentswhichareidentifiedateachiterationintheWS05al- dynamicspectrummodel.Forthenewalgorithmthelargestvalue
gorithmdependontheelectricfieldmodelatthatpoint,sothesys- oftheresidualpowerisonly11dBabovethemeannoisepowerin
tematicerrorsintheWS05modelarenotcompletelyeliminatedby thedata,whereasthepeakpowerinbothdataandmodelis63dB
postfactooptimisation—spuriouscomponentsremainintheopti- abovethemeannoisepowerinthedata.Thusthealgorithmiscer-
misedmodel,albeitatalowlevel.If,ontheotherhand,theelectric tainlyfreeofsystematicerrorsoverarangeof52dB.
field model at each iteration of the WS05 algorithm is optimised In fact the new algorithm has achieved noise-limited perfor-
(usingL-BFGS-B)thenthesespuriousfeaturescanbeminimised. manceandiscapturingthefulldynamicrangeofthedata;tosee
This approach has the additional benefit of simplicity in the pro- thisweneedonlyexaminethestatisticsofthenoisepower,shown
cessing of the data, requiring only one pass through the data and asahistograminfigure2.Iftheresidualswerepurelynoise-like
onesetofcode.Wethereforeimplementedanewiterativedecom- thentheexpectedprobabilitydistributionfunctionwouldbeanex-
positionalgorithm,inFORTRAN,whichemploystheL-BFGS-B ponential,becausetheresidualpoweristhesumofthesquaresof
packagetooptimisethemodelu˜j andfk.Thenewalgorithmdif- two independent variables each of which has the same Gaussian
fersfromWS05inthefollowingways: distribution. We can see from figure 2 that the residuals conform
closelytothisexpectation;wecanalsoseethatthepeakresidualis
(i) The parameters f (describing the intrinsic pulsar fluxes)
k notintroducedbyanysystematicerrorinthemodellingbutrather
areoptimisedonceforevery100newfieldcomponentswhichare
itissimplythetailofthenoisedistribution.Theresidualsactually
picked. The optimisation over {f } is done separately from that
k exhibitaslightlylowernoiselevelthanthedata(dashedline);this
over{u˜ }.
j canbeunderstoodbyreferencetofigure1.Thenoiseinthedata
isestimatedfromthefloorpowerlevelinthesecondaryspectrum,
specificallyanareainthecornerofthesecondaryspectrum(away
1 http://www.ece.northwestern.edu/∼nocedal fromanyobvioussignalpower)ischosenandthemeanpowerover
(cid:13)c 0000RAS,MNRAS000,000–000
4 Walker,Koopmans,Stinebring&vanStraten
Figure1.Data(top),model(middle),andresidualsbetweendataandmodel(bottom),foranobservationofPSRB0834+06ina1.56MHzbandcentredon
321.00MHz.ThedataweretakenwiththeAreciboradiotelescopeinconjunctionwiththeWAPPbackendsignalprocessingunitsonMJD53009;thereare
1024spectralchannels,and270timesampleseachof10secondsduration.Theleft-handpanelsshowdynamicspectra,whiletheright-handpanelsshowthe
correspondingsecondaryspectra(powerspectraofthedynamicspectra);therangeofthesecondaryspectrais±50mHzinDoppler-shiftand±327.6µsin
delay.Inversegrey-scale(blackispeakintensity)isusedinallcases;thetransferfunctionislinearforthedynamicspectra,andlogarithmicforthesecondary
spectra.Allsecondaryspectrashownherehavethesametransferfunction;thetransferfunctionsfordynamicspectraareidenticalinthecaseofdataand
model,whereastheoutputrangeisstretchedbyafactoroffiveforthedynamicspectrumresidualsinordertorevealtheirstructure.Incomparisonwithfigure
1ofWS05notethatthepresentfigureincludesthemodulatingeffectsoftheintrinsicpulsarfluxvariations.Wehavechosentodisplaythefullsecondary
spectra,includingnegativedelays,eventhoughthesespectraaresymmetricthroughtheorigin,sothatthestructurenearzerodelaycanbebetterseen.
(cid:13)c 0000RAS,MNRAS000,000–000
InterstellarHolography 5
Figure2.Histogramofthenumberofpixelsperbin(ofwidth6×10−8)as
afunctionofpower,forthetwo-dimensionalpowerspectrumoftheresid-
ualbetweenthemodeldynamicspectrumandthedata(solidline).Also
shownisthecorrespondinghistogramforthedata(dashedline).Theresid-
ualpowerdistributioncorrespondsverycloselytotheexpectedexponential
noisedistribution.Themeanpowerintheresidualsisslightlylessthanthat Figure3.Modelelectricfieldamplitudescorrespondingtothemodelspec-
calculatedfromthedatabecausethelatterestimateincludesasmallcontri- trashowninfigure1.Heretheamplitudesareshowningrey-scale,asa
butionfromsignalswhichareincorporatedintothemodel. functionofDoppler-shiftanddelay,withalogarithmictransferfunction;
peakintensityisblack.TheDoppler-shiftrangesoveratotalof100mHz
(270pixels)andthedelayrangesoveratotalof655µs(1024pixels).Ab-
thisareaiscomputed.Themodelsecondaryspectrumalsoexhibits solutevaluesofdelayandDoppler-shiftareinprincipleunknownbutitap-
afloorpowerlevel,eventhoughthemodelisintendedtorepresent pearssensibletoassigntheoriginofcoordinatestothecentreoftheimage
“signal”ratherthan“noise”;thispowerisnotpresentintheresid- (whichisthepeakfieldamplitude)inthisparticularcase.Itisevidentfrom
ualdynamicspectrumsothemeanpowerlevelintheresidualsis thisfigurethattheimagehasbeenwellseparatedfromitscomplexconju-
gate,whichwouldappearasaninvertedparabola.Alsonotableisthelow
slightlylowerthaninthedata.
levelofpoweraroundzerodelay,whichindicatesthattheintrinsicpulsar
Themodelelectricfieldstrengthdeterminedbyournewalgo-
fluxmodulationshavebeenaccuratelydeterminedbythealgorithm.
rithmisshowninfigure3.Visuallythisresultissimilartothatob-
tainedbyWS05(theirfigure2),withthenotableexceptionthatin
WS05allfieldamplitudeswerefixedatzerointheregionτ <0(so ferencesaboutthepropertiesofthescatteringmediumwhichgives
thatregionisnotdisplayedintheirfigure2).Onphysicalgrounds risetothedatashowninfigure1.Exceptingthetemporalanalysis
thisareaisexpectedtobefreeofastronomicalsignals.However, in§4.1ourdiscussionisonlyqualitative;thatisbecausewehave
in-line holography generally involves a certain amount of confu- constructedourholographicimageindelay-Dopplerspace(i.e.the
sionbetweentheimageanditsconjugate,becausetheobservable Fourierspaceconjugatetothefrequency-timespaceinwhichthe
quantityisusuallytheintensityU∗U whichisinvariantunderthe dataarerecorded),whereasprogressinunderstandingthescatter-
replacementU → U∗.Inthepresentcontexttheconjugateimage ingmediumreliesonaknowledgeoftheimageintwo-dimensional
appearsatnegativedelays,becauseformingthecomplexconjugate spatial(angular)coordinatesandthereisnosimplewayofproceed-
of exp[2πi(ντj +ωjt] is equivalent to making the replacements ingfromtheformertothelatter.
τ → −τ ,ω → −ω .Inourfigure3itcanbeseenthatthereis
j j j j
littletraceofaninvertedparabolainthelowerhalfofthefigure–
onlyweakcomponentscanbeseenintheregionτ < 0–indicat-
4 PROPERTIESOFTHESCATTERINGMEDIUM
ingminorconfusionbetweentheimageanditsconjugate(seealso
§4.1).Thiscleanseparationisanotherindicationofthelowlevelof An important qualitative aspect of figure 3 is that power in the
systematicerrorsinherentinthenewalgorithm.Inthesamevein modelelectricfieldistightlyconcentratedaroundaparaboliclo-
wenotethatthereislittlepowerevidentinfigure3neartheτ =0 cus,withdelayproportionaltothesquareoftheDoppler-shift;this
locus(ahorizontallinetangenttotheapexoftheparabola),indicat- muchwasalreadyevidentinWS05.Thisconfirmsthattheunder-
ingthattheintrinsicpulsarfluxvariations,describedbythevarious lyingscatteringishighlyanisotropic,withascatteredimagewhich
f ,havebeenaccuratelymodelled. isverymuchlongerthanitiswide—aconclusionwhichhaspre-
k
Althoughtheoriginofcoordinates(i.e.τ =0,ω =0)infig- viouslybeenarrivedatfromconsiderationofthepropertiesofob-
ure3isinprincipleunknown(WS05),inthecaseofthisparticular served secondary spectra for this and other pulsars (Walker et al
datasetitappearssensibletoassigntheorigintothecentreofthe 2004;Cordesetal2006;TrangandRickett2007).
imagebecausethisiswherethelargestamplitudefieldcomponent Wecanalsoseefromfigure3thatsomepartsoftheparabola
is located and the great majority of the intensity in figure 3 lies showhighintensitylevelswhileothersarealmostcompletelyde-
abovethispoint—consistentwiththebrightestimagecomponent voidofsignal,andinsomeplacesthehigh-lowtransitionsarefairly
beinganunscatteredwave. abrupt.Abruptchangesaresuggestiveofwell-definedboundaries
Having now reached the point where we can form models tothescatteringregions.Fourintensityconcentrationsareseenat
whichareagoodmatchtothedataitisappropriatetodrawsomein- largedelayandpositiveDoppler-shift;thesecorrespondtothefea-
(cid:13)c 0000RAS,MNRAS000,000–000
6 Walker,Koopmans,Stinebring&vanStraten
(cf.equations1and2ofWS05).Thecorrespondingfieldintensity
is
I(τ,t)=U∗U =A(τ) + B(τ,t), (3)
where
B(τ,t):=Xu∗u δ(τ −τ )δ(τ −τ )exp[2πi(ω −ω )t],
ejek j k k j
j(cid:54)=k
(4)
describes the beating between waves of differing Doppler-shifts
(butthesamedelay)and
A(τ):=Xu∗u δ(τ −τ )
ejej j
j
Z ˛ ˛2
= dω ˛˛Ue(τ,ω)˛˛ , (5)
is the mean intensity impulse response function of the scattering
medium.ThefunctionA(τ)convolvedwiththeintrinsicpulsepro-
fileyields(uptoanormalisationfactor)theaveragepulseprofile
forthisobservation.Forourimage,Ue,themeanintensityimpulse
response,A,isshowninthetoppaneloffigure4.Forthosedelays
where Ue includes one wave component which has a much larger
amplitudethantheothercomponentsatthatdelaythebeatterms
willallberelativelysmalland|B|willbesmallcomparedwithA.
Ingeneral,however,Bisnotnegligible;wedeferconsiderationof
Btolaterinthissection.
There are several recognisable features of the A(τ) deter-
minedfromourmodel:thepeakatzerodelayisinaccordwiththe
physicalexpectationthatabrightimageshouldformatthedelay
minimum;theparabolicarcseeninfigure3isresponsibleforthe
extendedscatteringtailstretchingoutbeyondτ =100µs;andthe
peaksnearτ =140,160,230,270µsinfigure4correspondtothe
featureslabelled“A,B,C,D”,respectively,byHilletal(2005)—
thesefeaturesareseenasdiscreteintensityconcentrationsatlarge
Figure4.Themeanintensityimpulseresponseofthescatteringmedium: delayandpositiveDoppler-shiftinfigure3.
toppanel,asdeterminedfromtheholographicimageshowninfigure3; Atlargenegativedelaysinthetoppaneloffigure4weseean
bottompanel,asdeterminedfromtheholographicimageshowninfigure3 intensitylevel∼ 10−5,whereasphysicallyweexpectzerointen-
withlowamplitudecoefficients(|u˜k|<0.004)settozero. sitybecausewavepropagationcanonlyintroducepositivegroup-
delays.Thisissimplyduetonoiseinthereconstruction;notealso
thatasimilarfloorintensitylevelispresentatlargepositivedelays.
turesnamed“A,B,C,D”byHilletal(2005)inasecondaryspectrum
Themodelelectricfieldinevitablyincludesnoisebecausethereis
analysis of data spanning more than three weeks. These features
nocleardistinctionbetweennoiseandsignalcontributionstothe
arepresumablyduetolocalisedplasmaconcentrations;itisnotyet
inputdata.BysettingthelowamplitudecoefficientsofUe tozero
clearwhethertheseconcentrationsarediffractingorrefractingthe
weareabletoeliminatemuchofthisintensityfloor.Specifically,
radio-wavesintothetelescope.Feature“A”liessignificantlyabove
settingthemodelcoefficientsu˜ tozerofor|u˜ |<0.004preserves
j j
thelocusoftheparabola;thisextradelaycouldbethewave-speed
alloftherecognisablesignalfeaturesinA(τ)butlargelyremoves
(dispersive)delayofahighcolumn-densitystructure—aninter-
thenoisefloor;theresultisshowninthelowerpaneloffigure4.
pretation whichcan in futurebe tested by comparingdata at two
At small negative delays the intensity level in figure 4 is up
differentfrequenciesobtainedatthesameepoch.
toanorderofmagnitudehigherthanthenoisefloor.Weinterpret
this as evidence of low-level confusion between the holographic
image, U, and its complex conjugate, U∗ — as noted in §3 we
4.1 Temporalresponseofthemedium expectthisconfusiontobepresentatsomelevel.Forapplications
wheresuppressionoftheconjugateimageiscriticalitispossible
Inferring the spatial structure of the scattering medium from our
to undertake the holographic image reconstruction entirely in the
imageofUe isnotasimpleexerciseanditisbeyondthescopeof
positive-delayhalf-space,butwenotethatthisapproachisexpected
this paper to attempt an in-depth analysis. On the other hand the
tobeproblematicifthebrightestimageisnottheimagewiththe
holographic image shown in figure 3 is well suited to determin-
leastdelay(WS05).
ing the temporal response of the medium which is introduced by
A useful characterisation of the temporal response of the
multi-pathpropagation.Theelectricfieldamplitudeasafunction
mediumistheintensity-weightedaveragedelay,∆,definedby
ofdelay,τ,andobservingtime,t,canbedeterminedbyforming
theinverseFouriertransformofequation1withrespecttoω: R∞dττI(τ,t)
X ∆(t)≡ R0∞dτI(τ,t) . (6)
U(τ,t)= u δ(τ −τ )exp[2πiω t] (2) 0
ej j j
Thetime-averageofthisquantity,(cid:104)∆(cid:105),providesuswithasimple
j
(cid:13)c 0000RAS,MNRAS000,000–000
InterstellarHolography 7
5 DISCUSSION
Over an interval of several months we expect PSR B0834+06 to
exhibitchangesδ(cid:104)∆(cid:105)inthemeandelay(cid:104)∆(cid:105)asthepulsarmoves
behinddifferentregionsofthescatteringscreen.Theelectricfield
image shown in figure 3 clearly demonstrates that this particular
screenhasaverypatchydistributionofscatteringmaterialsothat
observationsmadeatotherepochsareexpectedtoexhibitquitedif-
ferentdistributionsofpoweralongtheparaboliclocus.Asahypo-
thetical example we can imagine that at another epoch the holo-
graphicimagemightlooklikefigure3atpositiveDopplershifts,
but show no scattered power at negative Doppler shifts; in this
casethemeandelaywouldberoughlyhalfofwhatwehavemea-
sured.Wethereforeexpectthatoveranintervalofseveralmonths
PSRB0834+06willexhibitchangesδ(cid:104)∆(cid:105)∼(cid:104)∆(cid:105).Theanticipated
delaychangesδ(cid:104)∆(cid:105)∼15µsareverylargecomparedtothepreci-
sionwithwhichmillisecondpulsarscanbetimed(e.gvanStraten
Figure 5. The interstellar propagation delay, ∆(t) (solid line), as deter- et al 2001). B0834+06 is not a millisecond pulsar, and even if it
minedfromtheholographicimageshowninfigure3withlowamplitude were it would not be used for precision timing experiments pre-
coefficients(|u˜k| < 0.004)settozero.Alsoshownisthe(unweighted) cisely because this line-of-sight is known to exhibit very striking
meandelay,(cid:104)∆(cid:105) = 15.2 µs(dashedline),fortheobservation.Thedot-
effectsfrommulti-pathpropagationintheinterstellarmedium.Itis,
tedlineshowsaweightedmean(14.8 µs),withtheweightforeachtime
however,salutarytoseehowlargetheinfluenceoftheinterstellar
samplebeingequaltotheintrinsicpulsarflux,fk,asdeterminedbythe mediumcanbeonpulsearrivaltimes.Moreovertheline-of-sight
modellingproceduredescribedin§2ofthispaper.
to B0834+06, although unusual, is not unique — B1133+16 and
B1929+10, for example, appear to show similar, striking effects
(PutneyandStinebring2005).Norisitaverydistantpulsar,sothe
gauge of the influence of wave propagation on pulse arrival time
scatteringstructureswhichhavebeenrevealedinthepresentdata
forthisline-of-sightfortheparticulartimeandfrequencyintervals
are probably abundant in the interstellar medium. B0834+06 is,
coveredbyourdata.FortheimageUe showninfigure3wecom-
however,arelativelybrightsourceandother,fainterpulsarsmight
pute(cid:104)∆(cid:105) (cid:39) 17µs.However,thisvalueisclearlyanoverestimate
be viewed through similar scattering media without being recog-
becausenoiseinthereconstruction–thenoisefloorisseeninthe
nisedassuchbecausethescatteredintensityissmall.Thedataof
toppaneloffigure4–biasestheestimateupward.Abetterestimate
PutneyandStinebring(2005)supportthesepointsasmanyofthe
isavailableifweemploytheimageUe withlow-amplitudecoeffi-
pulsarswhichtheystudiedappeartoshowstrongscatteringarising
cients set to zero (as per the lower panel in figure 4); this yields
frommultiple,physicallydistinctregions,andtheynotethatvery
the result (cid:104)∆(cid:105) (cid:39) 15.2µs. Substantial contributions to this mean
sensitiveobservationsarerequiredinordertorevealthepresence
delayarisefromthewholeregion0<τ <120µs,withrelatively
ofthesemedia.Furthermoretheverystructurednatureofthescat-
minorcontributionsfromeachofthediscretefeatures-“A,B,C,
teringmediumseeninfigure3tellsusthatapulsarwhichshows
D”ofHilletal(2005)-seenatlargedelay;togetherthesefeatures
nomeasurableinterstellartimingdelaysatoneepochcouldexhibit
contributeabout10%ofthemeasured(cid:104)∆(cid:105).Thereisnocontribu-
largedelaysatotherepochs.Presumablyinterstellarscatteringme-
tionfromtheartifactsassociatedwiththeconjugateimageasthese
diaexhibitabroadrangeofphysicalproperties,withacorrespond-
occurintheregionτ <0andareexcludedbythedefinitiongiven
inglybroadrangeofinfluenceonpulsearrivaltimemeasurements,
inequation6.
andforanygivenpulsarweshouldnotaskourselves“whether”but
Finallywereturntotheinfluenceofthebeatterms,described
“atwhatlevel”isanextendedscatteringtailpresent?
byB;thesetermscausethepropagationdelay∆(t)tovaryoverthe
courseofour45minutestretchofdata.Asnotedearlier,thebeat Acommonlyusedstrategyformitigatingtheinfluenceofthe
termsarenotnegligibleandforourdatawefindthattherearesub- interstellar medium on pulsar timing experiments is to undertake
stantialvariationsin∆of±40%aroundthemeanvalue(cid:104)∆(cid:105);these the timing observations at high radio frequencies. The rationale
variationsareplottedinfigure5.Largegradientsin∆(t)areseen forthisisbasedontheassumptionthatthescatteringoriginatesin
infigure5sothat,forexample,apulsearrivaltimemeasurement distributed Kolmogorov turbulence, for which the expected prop-
atthestartofourobservationsandonetaken5minuteslaterwould agation delay falls very rapidly with frequency. Distributed Kol-
havedifferedbyapproximately11 µs.Weemphasisethatthebe- mogorov turbulence is not a good model for what we see in fig-
haviourseeninfigure5isspecifictothisline-of-sightandtothe ure 3. Currently we are not able to predict what interstellar tim-
particularsetoffrequenciescoveredinourdata.Weexpectthatthe ingperturbationwouldbemeasured,forthispulsarandthisepoch,
temporalvariationsin∆wouldhavebeensmallerifourdatahad at radio-frequenciesoutside the1.56 MHz observingband ofthe
coveredagreaterbandwidththanthe1.56MHzusedhere.Witha presentdata.Itislikelythattheinterstellardelayswouldbesmaller
broaderobservingbandwidthwewouldhavefinerdelayresolution at higher frequencies – because the cold plasma refractive index
inourimageUe,sowewouldbebetterabletoseparatecomponents declines with frequency and so the scattering angles decrease at
which have similar values of τ and these separated components higherfrequenciesforanygivenplasmastructure–butwecannot
wouldnotbeat(seeequation4),sotheimportanceofBwoulddi- say how much smaller. From data spanning many months, Hem-
minish. It is beyond the scope of this paper to attempt a detailed bergerandStinebring(2008)haveusedasecondaryspectrumanal-
analysisofdelayvariations;herewesimplynotethatthetypeof ysis to estimate multipath propagation delays at frequencies be-
behaviourseeninfigure5ispotentiallyproblematicforprecision tween1,150MHzand1,500MHz,measuringvalueswhichrange
pulsartiming. fromonetotwoordersofmagnitudesmallerthanourestimateof
(cid:13)c 0000RAS,MNRAS000,000–000
8 Walker,Koopmans,Stinebring&vanStraten
(cid:104)∆(cid:105).Theirdatarefertoadifferentline-of-sight(PSRB1737+13) 7 ACKNOWLEDGMENTS
andcannotbedirectlycomparedwithourresult;however,theirfre-
Thanks to Don Backer, Barney Rickett and Bill Coles for help-
quencycoverageisgreatenoughtoallowthemtoestimatethefre-
ful discussions. MAW thanks R. Ramachandran for reminding
quencydependenceofthepropagationdelayfortheirobservations.
him of the correct terminology for the reconstruction undertaken
Theyfindthatthepropagationdelayisconsistentwithapower-law
by WS05. At the University of Sydney this work was supported
frequencydependenceofthemeanpropagationdelay,withpower-
by the Australian Reseach Council, at the Kapteyn Institute, AS-
lawindex−3.6±0.2—slightlylesssteepthanexpectedfordis-
TRONandLeidenbytheNetherlandsOrganisationforScientific
tributedKolmogorovturbulence(power-lawindex−4.4).
Research,andatOberlinbytheNationalScienceFoundation.
Thepossibilityoflarge,transientinterstellarpropagationde-
lays makes it prudent to quantify the interstellar delays at each
epoch where an accurate pulsar timing measurement is desired.
In this paper we have shown how those delays can be quantified
whenarecordingofthepulsardynamicspectrumisavailable.The REFERENCES
techniquewehavepresentedhasthemeritofbeingabletoaccu-
ByrdR.H.,LuP.,NocedalJ.,ZhuC.Y.1995SIAMJ.Sci.Comp.
ratelydeterminetherelativepropagationdelays,evenforcomplex
16,1190
scatteringgeometries,atanyepochofobservation.Thetechnique
CordesJ.M.,RickettB.J,StinebringD.R.,ColesW.A.2006ApJ
doeshaveitslimitationsthough.Foremostamongtheseisthatthe
637,346
delaysareallrelativemeasurementsandtheoriginofcoordinates
FienupJ.R.1982Appl.Opt.21,2758
(i.e.τ = 0)mustbedeterminedbysomeothermeans;removing
GwinnC.R.etal1997ApJL483,L53
thislimitationisahighprioritygoalforfuturework.
HembergerD.A,StinebringD.R.2008ApJL,674,Lxxx
Forcaseswherethegreatmajorityofthescatteringarisesina
HillA.etal2005ApJL619,L171
single,thinscreenweanticipatethatthekeytoprogressliesincon-
Ho¨gbomJ.A.1974A&AS15,417
structingtheholographicimagebymodellingthestructureofthe
ManchesterR.N.2007“40yearsofpulsars”(astro-ph:0710.5026)
phase screen, rather than forming the image Ue in delay-Doppler Nocedal J., Wright S.J. 1999 “Numerical Optimization”
space as we have done in this paper. Modelling the phase screen
(Springer:NewYork)
iscomputationallymuchmoredemandingbutithasseveraladvan-
PutneyM.L.,StinebringD.R.2006ChJA&A6(Sup.2),233
tagesovertheapproachwehavetakenhere.Firstitallowsthein-
RickettB.J.1991ARAA28,561
terstellarpropagationdelaystobedeterminedonanabsolutescale,
StinebringD.R.etal2001ApJL549,L97
becausethepropagationgeometryisfixedinthemodel.Secondly
TrangF.S.,RickettB.J.2007ApJ661,1064
itallowsustocalculatethewavefieldforvariousdifferentobserver
van Straten W., Bailes M., Britton M., Kulkarni S.R., Ander-
locations,thuspermittingdirectcomparisonsbetweendifferentob-
sonS.B.,ManchesterR.N.,SarkissianJ.2001Nat412,158
serving epochs and between different telescopes (e.g. Very Long
WalkerM.A.,MelroseD.B.,StinebringD.R.,ZhangC.M.2004
BaselineInterferometry)atthesameepoch.Thirdlyitpermitsdi-
MNRAS354,43
rectcomparisonbetweendataatdifferentfrequencies.Andbythe
WalkerM.A.,StinebringD.R.2005MNRAS362,1279(WS05)
sametokentherewouldbenodifficultymodellingdatatakenover
WolszczanA.,CordesJ.M.1987ApJL320,L35
a single, wide bandwidth even if the phase screen has large spa-
tialvariationsinwavedispersion.Bycontrastthepresentapproach
is not ideal for wide bandwidth data because of smearing due to
differentialdispersivedelays,acrosstheband,fromarangeofdis-
persionmeasures.Finally,ifwemodelthephasescreenitselfthe
resultsimmediatelyprovidepowerfulconstraintsonmodelsofthe
scatteringmedium,becausethephasestructuretellsustheelectron
column-densitystructureandtheline-of-sightmagneticfieldiffull
polarisationinformationisrecorded.
6 CONCLUSIONS
Interstellar holography is a precise new technique which affords
detailed insights into the interstellar propagation of radio pulsar
signals. The holographic image reconstructed from our test data
reveals a complex scattering structure whose physical nature is
unknown at present. Holographic imaging can be used to deter-
minetheinfluenceofmulti-pathpropagationonpulsearrivaltime
measurementsandthustocorrectforthesepropagationdelays.In-
terstellarpropagationdelaysareunpredictableandcanpotentially
make large contributions to the systematic errors in pulse arrival
time measurements. It is therefore prudent to make provision for
holography in every instance where an accurate measurement of
theunperturbedarrivaltimeisdesired.
(cid:13)c 0000RAS,MNRAS000,000–000