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