Roche tomography of cataclysmic variables: I. artefacts and techniques

12 Pages·0.54 MB·English

by  C. A. Watson

Mon.Not.R.Astron.Soc.000,000–000 (0000) Printed1February2008 (MNLATEXstylefilev1.4)
Roche tomography of cataclysmic variables: I. artefacts
and techniques
C.A. Watson and V.S. Dhillon
Department of Physicsand Astronomy, University of Sheffield, Sheffield S3 7RH, UK
1
0 Accepted for publication in the Monthly Notices of the Royal Astronomical Society
0
2
23/01/2001

n
a
J
ABSTRACT
9
Roche tomography is a technique used for imaging the Roche-lobe filling secondary
2
starsincataclysmicvariables(CVs).InordertointerpretRochetomogramscorrectly,

onemustdeterminewhetherfeaturesinthereconstructionarereal,orduetostatistical
1
orsystematicerrors.Weexploretheeffectsofsystematicerrorsusingreconstructionsof
v
3 simulateddatasets andshowthat systematicerrorsresultincharacteristicdistortions
1 ofthefinalreconstructionsthatcanbeidentifiedandcorrected.Inaddition,wepresent
5 a new method of estimating statistical errors on tomographic reconstructions using a
1 Monte-Carlobootstrappingalgorithmandshowthismethodtobemuchmorereliable
0 than Monte-Carlo methods which ‘jiggle’ the data points in accordance with the size
1 of their error bars.
0
/ Key words: binaries: close – novae, cataclysmic variables – stars: late-type – stars:
h
imaging – line: profiles.
p
-
o
r
t
s
1 INTRODUCTION CV secondaries can be used to study the solar-stellar con-
a
: nection. It is well known that magnetic activity in isolated
v The secondary, Roche-lobe filling stars in CVs are key to lower-main sequence stars increases with decreasing rota-
Xi our understanding of the origin, evolution and behaviour tion period (e.g. Rutten 1987). The most rapidly rotating
of this class of interacting binary. To best study the sec- isolated stars of this type have rotation periods of ∼ 8
r
a ondary stars in CVs we would ideally like direct images of hours, much slower than the synchronously rotating sec-
thestellarsurface. Thisiscurrentlyimpossible, however,as ondary stars found in most CVs. One would therefore ex-
typical CV secondary stars have radii of 400 000 km and pect CVs to show even higher levels of magnetic activity.
distances of 200 parsecs, which means that to detect a fea- There is a great deal of indirect evidence for magnetic ac-
ture covering 20 per cent of the star’s surface requires a tivity in CVs – magnetic activity cycles have been invoked
resolutionofapproximately1micro-arcsecond,10000times to explain variations in the orbital periods, mean bright-
greaterthanthediffractionlimitedresolution oftheworld’s nesses and mean intervals between outbursts in CVs (see
largesttelescopes.Rutten&Dhillon(1994)describedaway Warner 1995). The magnetic field of the secondary star is
around this problem using an indirect imaging technique also believed to play a crucial role in angular momentum
called Roche tomography which uses phase-resolved spectra loss via magnetic braking in longer-period CVs, enabling
to reconstruct the line intensity distribution on the surface CVs to transfer mass and evolve to shorter periods. Oneof
of the secondary star. the observable consequences of magnetic activity are star-
Obtaining surface images of the secondary star in CVs spots, and their number, size, distribution and variability,
has far-reaching implications. For example, a knowledge of as deduced from surface images of CV secondaries, would
the irradiation pattern on the inner hemisphere of the sec- providecritical testsof stellar dynamomodels inahitherto
ondary star in CVs is essential if one is to calculate stel- untested period regime.
lar masses accurately enough to test binary star evolution However, determining whether features in Roche to-
models (see Smith & Dhillon 1998). Furthermore, the ir- mography reconstructions are real or not is problematic as
radiation pattern provides information on the geometry of theresultingmapsarepronetobothsystematicerrors, due
theaccreting structures around the white dwarf (see Smith to errors in the assumptions underlying the technique, and
1995). Perhaps even more importantly, surface images of statisticalerrors,duetomeasurementerrorsontheobserved
(cid:13)c 0000RAS
2 C.A. Watson and V.S. Dhillon
datapoints.Itisthereforeessentialtoquantifytheeffectsof line profile models (Unruh & Collier Cameron 1995), chro-
theseerrorsonRochetomogramsinordertoproperlyassess mosphericemission (Unruh&CollierCameron 1997; Bruls,
the reality of any features present. In this paper we study Solanki&Schussler1998),gravitydarkeninganddifferential
the artefacts produced by errors in the input parameters rotation (Hatzes et al. 1996). In this section we explore the
using reconstructions of simulated datasets, and present a effectsofsystematicerrorsonRochetomograms ofCVsus-
technique for estimating the statistical errors on the recon- ingsimulated reconstructions inasimilar mannertoMarsh
structions. Subsequent papers will present the application & Horne (1988) in their treatment of systematic errors on
of Roche tomography to real datasets. Dopplertomograms of CV accretion discs.
3.1 The test image
2 ROCHE TOMOGRAPHY Alloursimulations are based on a single test image (Image
A, figure 1). This image mimics the effect of irradiation of
Roche tomography (Rutten & Dhillon 1994; Rutten &
the secondary star by the compact primary, taking into ac-
Dhillon 1996) is analogous to the Doppler imaging tech-
countthedistancetothecompanionandtheincidenceangle
niqueusedtomaprapidlyrotatingsinglestars(e.g.Vogt&
at the surface of the star. The intensity of the side of the
Penrod 1983), detached binary stars (Ramseyer, Hatzes &
starthatisnotirradiatedissettotheintensityatthetermi-
Jablonski1995)andcontactbinaries(Maceronietal.1994).
nator. Also, shielding by (for example) an accretion stream
However,incontrast totheDopplerimaging technique,the
is mimicked by a narrow band stretching along theequator
continuum is ignored in Roche tomography because of the
unknown and variable contribution of the accretion regions fromtheinnerLagrangian (L1)pointtotheterminatorand
covering3.6%ofthetotalsurfacearea.Theintensityofthis
tothespectrum,whichforcesRochetomographytomapab-
band is equal to the intensity of the non-irradiated side of
solutelinefluxes.Apartfrom this,theprinciplesonlydiffer
thesecondary.
in detail.
In addition, an array of ten completely dark spots are
In Roche tomography we assume values for the binary
distributed around the star. Four are located on the trail-
parameters and that the secondary is Roche-lobe filling,
inghemispherecenteredatlatitudesofapproximately+50◦,
locked in synchronousrotation and hasa circularised orbit.
+25◦, 0◦ and –25◦, all with different longitudes. Another
We then model the secondary as a series of quadrilateral
four are located on the leading hemisphere at latitudes of
tilesofapproximatelyequalarealyingonthecriticalRoche
+50◦, +25◦, 0◦ and –36◦, all with the same longitude. The
surface.Eachtileorsurfaceelementisassignedacopyofthe
twofinalspotsarelocatedatthenorthpole(+90◦)andthe
local (intrinsic) specific intensity profile convolved with the
instrumentalresolution.Theseprofilesarescaledtotakeinto L1 point. Each spot covers 0.7% of the total surface area
accounttheprojectedarea,limbdarkeningandobscuration of the secondary, with the exception of the L1 spot which
covers 1%.
and then Doppler shifted according to the radial velocity
Finally, the geometry of the Roche-lobe itself is given
of the surface element at a particular phase. Summing up
the contributions from each element gives the rotationally bythebinary parametersM1 = 1.0M⊙,M2 =0.5 M⊙ and
an orbital period of 2 hours. Unless otherwise stated the
broadened profile at that particular orbital phase.
orbital inclination used was 60◦. In total, the image con-
Inordertoreconstructasurfaceimageofthesecondary,
sists of 3008 elements which,with thesebinary parameters,
the‘inverse’ofthisprocessmustbecarried out.Wedothis
givesamaximumradialvelocitydifferencebetweenanytwo
byiterativelyvaryingthestrengthsoftheprofilecontributed
neighbouring elements of 18.1 kms−1.
byeachelementassumingthattheshapeoftheintrinsicline
profile does not change. This iterative procedure is carried
out (assuming no contribution from the accreting regions) 3.2 Synthetic datasets and fits
until a satisfactory fit to the observed data, defined by the
χ2 statistic (i.e. χ2 ∼ 1), is achieved. As there are a large Unlessotherwisestated,aperfectsynthetictrailedspectrum
was computed assuming an instrumental resolution of 20
number of different intensity distributions that all predict
kms−1 and a negligible intrinsic line width. The whole or-
trailed spectra consistent with the observed one, an addi-
bit was sampled over 50 evenly spaced phased bins, with a
tional constraint is required. This is found by selecting the
velocity bin-width of 10 kms−1 and a negligible exposure
mapofmaximumentropywithrespecttoasuitabledefault
time. No contribution from an accretion disc or other part
map and is performed using the maximum entropy algo-
of the system was considered. A sample trailed spectrum is
rithmdevelopedbySkilling&Bryan(1984).Furtherdetails
showninfigure2andthebestfittothisdatasetisdisplayed
on the Roche tomography process can be found in Dhillon
in Image B, figure1.
& Watson (2000).
In each case the initial guess at the fit was one of uni-
form intensity distribution and a uniform default map was
used except in section 3.14. Due to the nature of the tests,
the fits could not proceed to the same final χ2 on each oc-
3 SYSTEMATIC ERRORS AND ARTEFACTS
casion. Thus the optimal fit was determined by minimising
Therehavebeenanumberofstudiesexploringtheeffectsof χ2 untilthespot features were resolved.
systematic errors on Dopplerimages of single stars. Forex- The reconstructions are plotted using a linear grey-
ample, the significance of so-called polar spots (large, long- scale.Avalueof100correspondstothemaximumintensity
lived, high latitude spots straddling the polar regions) has on the original test image (Image A, figure 1). Any inten-
been tested for robustness against errors such as incorrect sity value ∼ 22% greater than the maximum intensity of
(cid:13)c 0000RAS,MNRAS000,000–000
Roche tomography of CVs: I. artefacts and techniques 3
ImageA.Thetestimage ImageB.Bestfit
ImageC.Systemicvelocity ImageD.Flare
ImageE.Limbdarkening ImageF.Velocitysmearing
ImageG.Eclipsebyanaccretiondisc ImageH.Incorrectsecondarymass
Figure 1.Theeffects ofsystematicerrorsonthetestimage.
(cid:13)c 0000RAS,MNRAS000,000–000
4 C.A. Watson and V.S. Dhillon
ImageI.Incorrectinclination ImageJ.The‘mirroring’effect
ImageK.Phaseunder-sampling ImageL.Ring-likestreaksduetophaseunder-sampling
ImageM.Incomplete phasecoverage ImageN.Resolution
ImageO.Uniformdefaultmap ImageP.Gaussian-blurrdefaultmap
Figure 1 (continued).Theeffects ofsystematicerrorsonthetestimage.
(cid:13)c 0000RAS,MNRAS000,000–000
Roche tomography of CVs: I. artefacts and techniques 5
notallownegativeimagevalues,whichmeansithasgreater
freedom in setting the intensity (and hencethe location) of
thebrightband.Thebrightbandthereforebordersthedark
band,which lies on theequator.
3.4 Time variability
One of the key assumptions of our model is that features
on the surface of the secondary do not vary during the
course of an observation, e.g. a spot does not suddenly ap-
pear/disappear at a particular binary phase. Here we ex-
plore how a variation in the line flux from one region of
the secondary may effect the reconstruction by simulating
the presence of a flare. The flare is situated on the equa-
tor of the leading hemisphere (phase 0.75) and displaced
Figure 2.Synthesized datasetproducedbythetestimage. slightly towards the back of the star. It is only present be-
tween phases 0.8 and 0.92 inclusive, has approximately six
timestheintensityofthesurroundingnon-irradiatedregion
the test image is flagged as white. This is because some
and covers 0.7% of thetotal surface area.
artefactsgenerate‘hot’elementswithhighintensitieswhich
Inthereconstruction(ImageD,figure1)abrightregion
would otherwise dominate thegrey-scale plots.
can be seen at phase 0.75. This region is at roughly the
We are now ready to explore the effects of systematic
samelongitudebutatamuchlowerlatitudethantheactual
errors on the Roche tomography reconstructions. In each
positionoftheflare.Thisshiftinlatitudeiseasilyexplained;
caseafakedatasetwasgeneratedthatwaseithercorrupted
duetotheorbitalinclination,regionsnearthesouthpoleare
insomeformorwasfittedassuminganincorrectparameter.
visible for shorter periods than the equatorial regions. By
placingabrightregion there,thefittingroutineistryingto
3.3 Systemic velocity γ accountforthelimitedphasecoverageoverwhichtheflareis
observed. However, in doing so the model must smear this
Thesystemic velocity istheradial velocity of thecentre-of-
feature in order to match the observed radial velocities of
mass of the binary system and represents the mean shift of
the flare. In addition, bright ‘ringing’ can be seen, which
the line from its rest wavelength. In cases where the wrong
is due to the limited phase coverage over which the flare is
valueforthesystemicvelocityisassumed,thereconstruction
observed;anexplanation ofthiseffect willbegivenin more
process will have difficulty in fitting the line profiles. As
detail in section 3.12.
a result, artefacts will be generated in order to allow the
modeltofittheprofiles aswell aspossible. Duetothehigh
quality of the synthetic datasets used, one might expect a
3.5 Limb darkening
small error in the systemic velocity to produce noticeable
artefacts. Image C (figure1) shows theeffects of ignoring a In section 2 we mentioned that the line profile is scaled to
systemicvelocityof+5kms−1presentinthedata.Although takeintoaccountlimbdarkening.Thefittingprocedurecan-
the spot features and, to some extent, the irradiated inner not,without prior information, account for limb darkening.
face of the secondary are still visible in the reconstruction, This is because limb darkening causes each element in the
dark bandsframed bybright stripes are formed around the maptohaveaphase-dependentintensityvariation.Wemust
equatorial regions. therefore supply the fitting procedure with an appropriate
These equatorial stripes can be explained as follows. limb darkening relation. In this section we show the effects
When fitting the trailed spectrum the data appears to be of assigning incorrect limb darkeningstrengths and laws.
‘skewed’toonesideoftheradialvelocityaxisandtherefore First, data was generated assuming a linear limb dark-
the model cannot possibly fit the edges of the line profiles. ening coefficient of 0.5. This was then fitted assuming no
Atthelower(ormostnegative)radialvelocitiesexpectedby limb darkening (Image E, figure 1). The spot features are
the model there is no data, the data that should have been readily visible, although a dark equatorial band is formed.
present at these velocities having been shifted to a higher This band can be understood because the model calculates
radialvelocity.Thusthemodelattemptstofitzerointensity muchhigherintensities at theedgesof theprofilesthan are
to the regions of the star that contribute to the line profile presentinthefakedataset.Asdescribed insection 3.3,this
at these points, resulting in a dark equatorial band. resultsindarkequatorialbandsinthereconstruction.Inad-
At theother extreme,thehigher (or most positive) ra- dition,theoverallintensitylevelofthemapislessthanthat
dial velocities in the line profiles are shifted beyond theve- ofthetestimageduetothereductioninlightcausedbythe
locities expected by the model. Thus the intensities found limb darkening. Conversely, if one assumes too high a limb
atthehighestradialvelocitiesarefargreaterthanexpected darkening coefficient then the opposite occurs. The model
bythemodelandtheoppositescenarioemerges.Themodel calculates much lower intensities at the edges of the pro-
nowattemptstofithigh intensitiestotheregions that con- filesandinordertocompensatethefittingroutineproduces
tributetothelineprofilesatthesepoints.Asonlytheequa- bright bandsaround theequator.
torial regions ever contribute to the edges of the line pro- Foroursecondtest,weassumedanincorrectlimbdark-
files, a bright equatorial band results. The algorithm does ening law. Data generated using a linear limb darkening
(cid:13)c 0000RAS,MNRAS000,000–000
6 C.A. Watson and V.S. Dhillon
coefficient of 0.5 was then fitted assuming quadratic limb
darkeningco-efficients of a = –0.025, b = 0.635 and
I(µ)=I1(cid:0)1−a(1−µ)−b(1−µ)2(cid:1),
where µ = cosγ, γ is the angle between the line of sight
andtheemergent fluxandI1 isthemonochromatic specific
intensity at µ=1. The reconstruction showed a bright band
around the equatorial regions where the fitting routine has
had to compensate for the increased limb darkening pre-
dicted by thequadratic limb darkeninglaw in this case.
3.6 Velocity Smearing
Velocity smearing is the degradation in spectral resolution
duetothecombinedeffectsofafiniteexposuretimeandthe
Figure 3. Solid line: line profile at phase 0, computed assum-
orbital motion of thesecondary star. Velocity smearing can
ing zero exposure time. Crossed line: line profile around phase
be ignored if it is insignificant compared to the instrumen- 0, computed assuming an exposure time twice that required by
tal resolution. To obtain adequate signal-to-noise, however, equation1.
exposure times, texp, are usually selected which match the
resulting velocity smearing to the instrumental resolution,
vres, via theformula broadened.Thisisduetothealgorithmattemptingtomatch
thewidth of the bumpsin theline profile which, in the ab-
texp ∼Pvres/2πKR, (1) sence of intrinsic broadening, it can only do by increasing
thesize of thespots in thereconstruction.
where P is the orbital period of the binary and K is
R
theradial-velocity semi-amplitude of thesecondary star. In
cases where longer exposure times are unavoidable, signif-
3.8 Eclipse by an accretion disc
icant velocity smearing will occur. It is possible to correct
fortheeffectsofvelocitysmearing inRochetomographyby ThehighestinclinationCVsmayshowsignsofaneclipseof
calculating the line profiles at intermediate phases during the secondary star by an accretion disc around phase 0.5.
an exposure and averaging the results. This slows the iter- Thiscanbeaccountedforduringthereconstructionprocess
ations, however,soin thissection weexplorewhat happens either byremoving the affected spectra (as commonly done
if velocity smearing is neglected during the reconstruction in Doppler tomography during primary eclipse), or by in-
process. cludingtheeclipseinthealgorithm.Ifneitherofthesesteps
Image F (figure 1) shows what happens when an ex- aretaken,artefactswillbeintroducedinthereconstructions,
posure time twice as long as that required by equation 1 as described in this section.
is assumed when generating fake data and then ignored in A model including a cylindrical accretion disc of ra-
thereconstructionprocess.Onceagainthelargesteffectsare dius0.5R andheight0.1R wascreatedandasynthetic
L1 disc
seen in the equatorial regions. Figure 3 shows why this oc- trailed spectrum including the eclipse by the disc was gen-
curs.Theradialvelocitiesoftheedgesofthesmearedprofile erated assuming i = 80◦. Image G (figure 1) shows the re-
extend beyond the velocities expected when the exposure construction if the eclipse by the accretion disc is ignored.
time is ignored. For the same reasons as discussed in sec- Asone might expect,the artefacts generated by theeclipse
tion3.3,thisresultsinbrightanddarkartefactsaroundthe are most apparent around the L1 point. In particular, the
equator. Figure 3 also shows how bumps in the line profile southern regions (where light is blocked by the eclipse) ap-
duetostar-spots(e.g.thediparound0kms−1)aresmeared, pear dark, and there are alternating regions of bright and
resultinginablurringofthestar-spotsinthereconstruction. dark bands radiating from the L1 point. These are due to
the attempts of the fitting algorithm to match the eclipse
phases, resulting in a reduction of intensity on the inner
3.7 Intrinsic line profile width hemisphere which is compensated for at other phases by
the bright bands. Note that the spot features are well re-
Wehaveso far assumed that theintrinsic line-profilewidth
constructed, but exhibit the mirroring effect described in
is negligible in comparison with the rotational broadening.
section 3.11 dueto the high inclination.
This may not always be thecase and the width or shape of
the intrinsic line profile may be important. Here we exam-
ine the effects of assuming an incorrect intrinsic line-profile
3.9 Masses
width.
Data was produced assuming 20 kms−1 instrumental A knowledge of the binary masses are essential in order to
resolution as before but with an intrinsic profile width of define the geometry of the secondary star model. Very few
45kms−1 (FWHM).Thereconstructionwascarriedoutas- accurate mass determinations are available for CVs, how-
sumingnointrinsicline-profilewidthandtheresultwasvery ever,andthoughwecanuseRochetomographytodetermine
similar to Image F (figure 1). Once more the edges of the them (see section 4) we still need to know how an error in
linescannot befittedandbandingishenceseen aroundthe themasses may effect thereconstructions.
equator. The reconstruction also shows that the spots are Image H (figure 1) shows the effect of underestimating
(cid:13)c 0000RAS,MNRAS000,000–000
Roche tomography of CVs: I. artefacts and techniques 7
the secondary mass by 0.1M⊙. The effect of reducing the
secondary’s mass whilst keeping the period and inclination alsoshowsthatthefeaturesbecomelessapparentduetothe
the same is to shift the model to higher radial velocities. factthattheyaresmearedovertwicetheareatheyoriginally
Hencethespotfeaturesonthereconstructionappearshifted covered.
towards the L1 point and lower radial velocities in order to
compensate for this. In addition, as the radial velocities of
thespotfeaturesintheoriginaldatasetnolongercorrespond 3.12 Phase sampling
toawelldefinedlocationonthesecondary,thespotsappear
Upuntilnowoursimulations havebeencarried out over50
smearedoveralargerarea.Atphases0.25and0.75theouter
evenlysampledphasebinscoveringawholebinaryorbit.In
hemisphereofthemodelisshiftedtohigherradialvelocities
practice, however, it may not always be possible to obtain
thanpresentinthedatasetandhenceadarkbandisformed
thisnumberofphasesduetosignal-to-noiserequirements,or
neartherear ofthestar. Atotherphases themodelcannot
theremaybegapsintheorbitalcoveragebecauseofweather
matchthewidthofthelineprofilesduetothereducedradius
conditions. Here we model such scenarios.
of the star and this results in a bright equatorial band.
First we explore how under-sampling in phase can af-
fectthereconstruction.ImageK(figure1)wasreconstructed
from 6evenly spaced phases (5independent)usingzero ex-
3.10 Inclination ⋆
posure lengths and it is dominated by ring-like streaks.
Up until now we have always assumed the correct orbital The effect is analogous to the streaks observed in Doppler
inclination(60◦)ineachofourtests.Inreality,however,we tomography (Marsh & Horne 1988) and can be understood
may not knowthe inclination accurately. Image I (figure1) by considering that lines of constant radial velocity on the
showstheeffectofassuming anincorrect inclination of55◦. secondarystarataparticularphasecanbeintegratedalong
Once again, the most obvious artefacts appear at the to construct a line profile. These lines of constant radial
equatorial regions. The lower orbital inclination shifts the velocity are ring-like in shape and if there are only a few
modeltolowerradialvelocities.Asaconsequence,thegreat- phases, or if the profile is particularly bright at a certain
est mismatch between the observed and computed trailed phase(asinthecase ofaflare;section 3.4), thestreakswill
spectra occur around phases 0.25 and 0.75. At the lowest notdestructivelyinterfere,leavingring-likeartefacts on the
radialvelocitiesexpectedbythemodelthereisnodataand Roche tomogram. To demonstrate this effect we produced
a similar scenario occurs to that discussed in section 3.3; mapswherethoseelementswiththesameradial velocityas
darkbandsarefittedaroundtheequatorontheinnerhemi- thespot features were represented as agrey-scale whilst all
sphereofthestartocompensateforthelackofdataatthese otherelementswereleftblank.Thiswasdoneforeachofthe
velocities.Brightregionsthenappeararoundthedarkband 6 phases in the dataset and the resulting images were then
inordertofittheintensitiesseenatotherphases.Athigher superimposed on top of one another so that regions where
radial velocities around phases 0.25 and 0.75, the observed streaksoverlap areshowndarker(ImageL,figure1).Itcan
intensityisfargreaterthanexpectedbythemodel,andthe be seen that the features in this image closely resemble the
opposite occurs on the outerhemisphere of thestar. artefacts in Image K (figure1).
Thespotfeaturesremainvisible,butaresystematically ImageM(figure1)showstheeffectofincompletephase
shiftedtowardstheouterhemisphereofthestar.Thisshifts coverage, with phases between 0.4 and 0.6 missing (which
them to regions of the star with higher radial velocities, may occur if certain phases are discarded due to, for ex-
which is necessary to counteract the shift of the model to ample, an eclipse by the disc; see section 3.8). There is no
lower radial velocities caused by assuming an incorrect in- problem recovering the features around the L1 point where
clination.Theoppositeoccursiftoohighavaluefortheor- the majority of the information is missing. Although the
bital inclination is assumed; the dark bands become bright other spot features are also readily recovered they appear
andthespotfeaturesaresystematically shiftedtowardsthe elongated or streaked due to the same process described in
innerhemisphere and lower radial velocities. theabove paragraph.
3.13 Instrumental resolution
3.11 The mirroring effect
The spectral resolution of the data governs the size of the
Thefollowingsectionsdescribeartefactsthatarenotstrictly
features that can beimaged, as shown by theblurred spots
dueto systematic errors but are a result of factors that are
in Image N (figure 1) where we have reconstructed the test
eitherwithinourcontrol(e.g.instrumentalresolution,phase
image using data of lower resolution (60 kms−1) than be-
sampling)orareinherenttothereconstructionprocessitself.
fore. The number of surface elements in the model should
We begin with the problem that, although the radial
bechosentomatchthespectralresolution.Ifthemodelhas
velocities contain latitudinal information, they cannot con-
toofewelements,theseparation betweenadjacent elements
strainwhetherafeatureislocatedinthenorthernorsouth-
invelocityspace(particularlyaroundphases0.25and0.75)
ern hemisphere. This information can only be supplied by
is greater than the resolution of the data. This means that
self-obscuration, but this becomes ambiguous at higher in-
clinations until, at an inclination of 90◦, a feature in the
northernhemisphereisself-obscured inanidenticalmanner ⋆
Although such a scenario would probably only occur if long
to a feature located at the same latitude in the southern exposuretimeswerebeingused, wehave assumedzeroexposure
hemisphere.Thisresultsin a‘mirroring’ of featuresin both length in order to separate the effects of velocity smearing (de-
hemispheres as demonstrated in Image J (figure 1), which scribedinsection3.6)andthoseduetophaseunder-sampling.
(cid:13)c 0000RAS,MNRAS000,000–000
8 C.A. Watson and V.S. Dhillon
therewillbepartsofthelineprofiletowhichnosurfaceele- Bymodellingtheshapeandmappingthesurfaceinten-
mentsinthemodelcontribute(ataparticularphase)andno sity distribution of the secondary star, Roche tomography
acceptablefitwillbefound.If,ontheother-hand,toomany can provideconstraints on the CV masses. This is because,
elementsinthemodelareused,therewill benoproblemin as discussed in section 3.9, incorrect values of M1 and M2
fittingthe data butit will slow thereconstruction process. will introduce artefacts in the reconstructions which will
lower the entropy of the final solution. The correct values
of the component masses are therefore those that produce
3.14 Default maps
themap of highest entropy.
This is demonstrated in figure 4, which displays the
Inthereconstructionprocessthefinalmapistheoneofmax-
‘entropylandscape’ ofthetestdatashownin figure2.Each
imum entropy (relative to an assumed default map) which
square in the entropy landscape corresponds to the maxi-
isconsistent withthedata.Ifthedataaregood thenthefi-
mum entropy obtained in a reconstruction assuming a par-
nalmapisconstrained bythedataandwillnot bestrongly
ticular pair of component mass values, a uniform default
influenced by the default. In this case the choice of default
map and iterating to the same χ2 value on each occasion.
makeslittledifferencetothefinalmap.Ifthedataarenoisy,
The reconstruction with the highest entropy was obtained
however,thedataconstraintswillbeweakandthemapwill
bestronglyinfluencedbythedefault.Thedefaultcanthere- for component masses of M1 = 1.0 M⊙ and M2 = 0.5 M⊙,
which areidentical tothemasses usedto construct thetest
fore be regarded as containing prior information about the
data. This shows the entropy-landscape method to be an
map,andregionsofthemapthatarenotconstrainedbythe
extremelyeffectivemethod ofdeterminingaccurate compo-
data will converge to thedefault value.
nentmasses which accounts for both geometrical distortion
Toinvestigatetheeffectofdifferentdefaultmapsonthe
and complicated intensity distributions on the secondary
reconstruction, it is necessary to degrade the quality of the
star.Notethatthecorrectinclinationwasusedforeveryre-
testdataset.Wethereforegenerateddataover15phasebins
with a resolution of only 60 kms−1 anda S/N ratio of only construction infigure4,butinpractice theinclination may
not be accurately known. A similar method to the entropy
50.Thedatasetwasthenfittedusingauniformdefaultmap
landscapecan thenbeusedtodeterminetheinclination or,
and a Gaussian-blurr default map. The former was created
in cases where the secondary eclipses the white dwarf, the
bysettingeveryelement in thedefault totheaverage value
inclination can be varied for each combination of compo-
of the map, and hence constrains large-scale surface struc-
nentmassesinordertoremainconsistentwiththeobserved
ture. The latter was created bysmoothing the map using a
eclipse width.
Gaussian-blurringfunctionandhenceconstrainssmall-scale
Figure 4 also shows a ridge of high entropy running
surface structure.
alongalineof(almost)constantq.Componentmasseslying
Image O (figure 1) shows the reconstruction using the
alongthisdiagonalarethereforemoredifficulttodistinguish
uniform default map and, as expected, it can be seen that
between than component masses lying perpendicular to it.
thespotsare moreblurredthanin thereconstruction using
The gradient of the diagonal can be derived by considering
the Gaussian-blurr default map (Image P, figure 1). Fea-
tures below –60◦ latitude are never visible and hence in that the best fit is obtained when the centre of the model
line profile matches the centre of the observed profile. The
the fitting routine they are assigned a value determined by
centre of the line profile is governed by the radial velocity
the default map. In the case of the uniform default map,
of the secondary, which relates directly to the distance of
this is the average element value of the map. In the case of
thecentre-of-massofthesecondaryfrom thecentre-of-mass
theGaussian-blurrdefaultmap,however,thevalueassigned
depends upon the intensity values of the neighbouring ele- of the binary, a2. Thus only those combinations of M1 and
mentsand,for elements that are not visible, this is close to M2 that maintain constant a2 can provide satisfactory fits
to the data. Using Kepler’s law and assuming a constant
zero intensity. This explains why the reconstruction using
the Gaussian-blurr default has a dark southern limb and, period and a2 we can show that the ridges in the entropy
landscape correspond to component masses which obey the
to compensate for this, a bright region above the L1 point.
relation
Hence,althoughtheuseofaGaussian-blurrdefault ismore
likely to resolve small features, it is important to be aware M1 =constant. (2)
of the effects that visibility may have on thefinal outcome. (M1+M2)2/3
This simple argument is complicated in the case of irra-
diation. Rutten & Dhillon (1994) found that the ridge of
4 DETERMINING PHYSICAL PARAMETERS:
high entropy was shifted towards slightly higher values of
THE ENTROPY LANDSCAPE
M2.Theirexplanationforthiswasthattheintensitydistri-
If the centre-of-light and centre-of-mass of the secondary bution was strongly peaked towards the L1 point, causing
are not coincident, the star’s radial-velocity curve will be a large difference with the uniform default map. This dif-
distortedinsomewayfromthepuresinewavewhichrepre- ference was reduced by increasing M2, allowing the bright
sentsthemotion of its centre-of-mass. The observed radial- region to shift away from the L1 point and to spread over
velocitycurvesofCVsecondariessufferfromthisdistortion, moreelements,reducingthedeviationbetweenthemapand
due to both geometrical effects caused by the Roche-lobe thedefaultwhilststillfittingthetrailedspectrumtowithin
shapeandnon-uniformitiesinthesurfacedistributionofthe its uncertainties. Despite this, the mass determination was
line strength due to, for example, irradiation. If these fac- still correct to within ∼2per cent.
tors are ignored, they will lead to systematic errors in the Anotherfeatureofnoteintheentropylandscapeoffig-
determination of thebinary star masses. ure4isthat,astheprimarymassincreases,thereisalarger
(cid:13)c 0000RAS,MNRAS000,000–000
Roche tomography of CVs: I. artefacts and techniques 9
5 STATISTICAL ERRORS
Althoughitispossibletopropagatethestatisticalerrorson
thedatapointsthroughtheiterativeprocessinordertocal-
culatethestatistical errorsoneachelementinthemap,the
resultingerror-barsareunreliableduetothefact thatnoise
in Roche tomograms is correlated. This correlation arises
from the projection of bumps in the profile along arcs of
constant radial velocity on thesecondary starand from the
effects of blurring by the default map. Doppler imaging of
single stars has largely relied on consistency tests, such as
comparingimagereconstructionsusingeithersubsetsofob-
servations(e.g.Barnesetal.1998)ordifferentspectrallines
(e.g. Strassmeier & Bartus 2000), in order to address this
issue. So far, the only attempt at a formal statistical error
analysis in Doppler imaging has been undertakenusing the
Occamian approach (Berdyugina 1998).
Figure 4. Entropy landscape of the test data presented in fig-
InthecaseofRochetomography,onewillundoubtedly
ure2. The contours arelines of constant entropy and have been
plottedtohighlighttheridgeofhighentropy. be working with lower signal-to-noise datasets and fewer
spectral lines than conventional Doppler imaging. For in-
stance, the well-studied K-star AB Dor has a mean visual
range of secondary masses that can be fitted satisfactorily. magnitudeof∼7(Cameron&Foing1997),whereasthesec-
This can be explained by generalising equation 2 to ondary stars in CVs never attain magnitudes of below ∼10
andareusuallyseveral magnitudesfainter thanthis.Roche
q+1∝M1/2,
1 tomograms of CV secondaries are therefore not as strongly
which shows that as the primary mass increases, the mass constrained as Doppler images of single stars by the input
ratio, q = M2/M1, must also increase in order to hold a2 data. This makes the determination of statistical errors of
constant. It can be shown that the secondary star radius is prime importance in Roche tomography, because only then
related to q and a2 via theequation can one test whether a surface feature is real or just a spu-
rious artefact dueto noise.
R2/a2 ∝q1/3(1+q)2/3, We have addressed the problem of statistical error de-
which shows that the radius of the secondary star must in- terminationinRochetomographyusingaMonte-Carlostyle
crease if q increases. This in turn means that the model approach. Monte-Carlo techniques rely on the construction
line-width is greater than the observed line-width and this ofalarge(typicallyhundreds,inthecaseofRochetomogra-
overlap allows a greater range of secondary star geometries phy) sample of synthesized datasets which have been effec-
for which a satisfactory fit can be obtained. In contrast, at tivelydrawnfromthesameparentpopulationastheoriginal
low primary masses the allowed value of q is decreased, re- dataset,i.e.asiftheobservationshavebeenrepeated many
sulting in a low secondary star radius and hence a small hundredsoftimes.Thislargesampleofsynthesizeddatasets
model line-width. Atthelowest primary masses in figure 4, is then used to create a large sample of Roche tomograms,
themodelline-widthsbecometoonarrowtofittheobserved resultinginaprobabilitydistributionforeachelementinthe
profile to the desired χ2, no matter what intensity pattern map.Themaindifficultywiththistechnique,asidefromthe
is mapped onto thesecondary star’s surface, explaining the demands on computer time, lies in the construction of the
abruptendof thehigh entropyridgeon thelower left-hand sampleofsynthesizeddatasets.Oneapproach(usedbyRut-
portion of figure 4. tenetal.1994inspectraleclipsemapping)isto‘jiggle’each
It is important to note that the maximum entropy al- datapointaboutitsobservedvalue,byanamountgivenby
gorithmweuserequiresthatallinputdataispositive.Asa its error bar multiplied by a number output by a Gaussian
result, it is necessary to invert absorption line profiles and, random-number generator with zero mean and unit vari-
in the case of low signal-to-noise data, add a positive con- ance. This process adds noise to the data, however, which
stant in order to ensure that there are no negative values means that the synthesized datasets are not being drawn
in the continuum-subtracted spectra. If this positive con- from the same parent population as the observed dataset
stant is not added, and any negative points are either set – noise is being added to a dataset which has already had
to zero or ignored during the reconstruction, the fit will be noise added to it during themeasurement process. In prac-
positively biased in the line wings, resulting in a broader tice,thismeansthatfittingthesampledatasetstothesame
computed profile and hence larger stellar masses in the en- level of χ2 as the original dataset is either impossible (i.e.
tropy landscape. It is a simple matter to account for this the iteration does not converge) or results in maps domi-
positive constant during the reconstruction process by em- nated by noise, which grossly overestimates the true error
ploying a virtual image element which contributes a single oneachelement.Wehaveverifiedthatthisisthecaseusing
value to all data points, effectively cancelling out the con- simulations similar to those described in section 5.1.
stant.Thisvirtualimageelementcanalsobeusedtocorrect Instead we have opted for thebootstrapping technique
for small errors in the continuum subtraction, which might ofEfron(1979); seeEfron&Tibshirani(1993) forageneral
otherwise systematically affect the masses derived from the introductiontothemethod.Weimplementedthebootstrap
entropy landscape. asfollows.Fromourobservedtrailedspectrumcontainingn
(cid:13)c 0000RAS,MNRAS000,000–000
10 C.A. Watson and V.S. Dhillon
datapointsweformoursynthesizedtrailedspectrumbyse-
lecting,atrandomandwithreplacement,ndatavaluesand
placing these at their original positions in the new synthe-
sized trailed spectrum. For points that are not selected we
settheassociatederrorbartoinfinityandtherebyeffectively
omitthesepointsfromthefit.Forpointsthathavebeense-
lected onceormorethan once(typically37% ofthepoints)
wedividetheirerror barsbythesquareroot of thenumber
of times they were picked. The advantage of bootstrap re-
sampling over jiggling is that the data is not made noisier
bytheprocessasonlytheerrorsbarsonthedatapointsare
manipulated (i.e. the data values remain unchanged). In so
doing, theamount of noise present in the data is conserved
anditisthereforepossibletofitthesynthesizedtrailedspec-
Figure 5. The test image used to construct the noisy dataset
tra to the same level of χ2 as the observed data, giving a
describedinsection5.1
much more reliable estimate of the statistical errors in the
maps.
5.1 Testing the bootstrap
In order to assess the viability of the bootstrap technique
we ran a simple test. A model with parameters M1 = 1.0
M⊙, M2 = 0.5 M⊙, i = 60◦ and P = 2.78 hours was cre- A B
atedofuniformintensitydistribution,withtheexceptionof
a single spot covering 0.8% of the total surface area on the
trailing hemisphere (see figure 5). From this model a syn-
thetictrailedspectrumwasgeneratedwith30evenlyspaced
phasebinsand60kms−1instrumentalresolution.Eachdata
pointinthespectrumwasthen‘jiggled’ by10%oftheorig-
inaldatavaluein ordertosimulate noise.Aclosefittothis
dataset was then carried out using thecorrect parameters. Figure6.TheRochetomographyreconstructionofthetestim-
As can be seen, the reconstruction (figure 6) shows ageinfigure5usingnoisydata.Althoughthespotcanbeseenat
many spurious artefacts due to noise. Fitting to a higher phase0.25,manyspuriousartefactsduetonoisearealsopresent,
χ2, whilst removing some of the lower-level noise artefacts, includingfeaturesAandB.
does not remove the spurious spot features marked A and
B in figure 6. The triangular symbols in figure 7 show the
intensities in each vertical slice through A, B and the real First, it can be seen that in regions where the data do not
spot. It can be seen that, from this information alone, it is constrain the intensity distribution (such as the southern
impossible to deducewhich spot features are real. hemisphere, which is never visible), the reconstruction al-
The curves in figure 7 show the confidence intervals ways converges to the default map value, which is why the
foundafterreconstructing200bootstrap samplesoftheim- confidence intervals converge around ‘SP’ in figure 7. Sec-
age.Asthedistributionofelementintensitiesfoundisoften ond, the slices through features A and B demonstrate how
non-normal and not centred on the intensity found in the noise is correlated in the Roche tomography reconstruction
original fit (see figure 8), it is incorrect to calculate a sum- process,becauseelementssurroundingthepeakintensityare
mary error statistic like the standard deviation. Instead we all shifted in thesame direction as thepeak itself.
take the mode of the distribution to represent the ‘most
probable’ value,anddefinethe95% confidenceinterval(for
example)astheregionwhichencloses95% ofthebootstrap
6 CONCLUSIONS
values above and below the mode. It can be seen that the
confidence intervals in figure 7 widen significantly around We have shown that any feature on a Roche tomogram
the noise features A and B, and encompass the true inten- mustbesubjectedtotwotestsbeforeitsrealitycanbecon-
sity level (indicated by the horizontal line). The confidence firmed. The first test is to determine whether the feature
intervals around the real spot feature, however, do not ex- is statistically significant and is performed using a Monte-
hibitsuchawideningandclearlydonotencompassthetrue Carlotechnique.Theconventionalmethod,wherehypothet-
intensity level of the non-spotted regions. This test there- ical datasetsareconstructed by‘jiggling’ thedatain accor-
fore confirms that the bootstrapping method works; it has dance with their error bars, was found to add noise and
correctly identifiedthereal spot featureand correctly iden- hencegrosslyover-estimatetheerrorsinthereconstruction.
tifiedfeaturesAandBasartefactsduetonoise.Inasimilar Thesolutiontothisproblemwasfoundbyusingabootstrap
test,statisticalerrorsdeterminedbythemethodof‘jiggling’ re-sampling algorithm and, through the use of simulations,
failed to distinguish between the real spot and features A we have shown this technique to correctly distinguish be-
and B. tween real features and artefacts due to noise. The second
There are two additional features of note in figure 7. testistocomparethefeaturewiththeappearanceofknown
(cid:13)c 0000RAS,MNRAS000,000–000