Happy New Year!

No, I'm not at work today but realized I never put in the final results from the exposure map redo.  So, the  new exposure map is now done.  However, it now appears that the best-fit column density for scattering is much bigger than I was expecting.  It's also not clear if the shape of the halo is perhaps indicating a non-smooth-distribution of dust, which would be exciting if true.  So, it looks like more work is needed, hardly a surprise.  One thing I can do is combine results with the GX5-1 ACIS data that's already been published and see how everything hangs together.  Also, I need to start work on the GX9+1 data as well, and redo all the steps I just did for GX5-1.  GX9+1 has the advantage that I also got some ACIS data as well as HRC-I, so I need to do all of the standard ACIS steps as well to get an energy-dependent halo.

So, more work to be done, although it will likely not happen until I get back from the AAS meeting in DC next week. 

Rethinking those first results...

Well, it always pays to check and recheck.  It looks like my original 'weights' on the exposure map file were done incorrectly; the correct method is shown in this thread. Redoing the calculation reduces the total effective area on axis from ~62 cm^2 to more like 46 cm^2, increasing the required flux to get the observed number of counts. 
As a side note, this is a major benefit for keeping a log file (mine are all called 'LOG') that shows all the steps necessary to recreate each major step in an analysis.  This way, when you find a mistake, it's reasonably easy to go back and redo everything without having to re-derive each step. 
So, now I'm re-running mkinstmap (fast) and mkexpmap (very, very slow...)

First results...

OK, the first surface brightness fit is ready to go.  There will probably be issues with it down the road, but this was my first shot.  Note that the fit uses the BARE-GR-B model from Zubko, Dwek & Arendt (2004) as the dust model.  There is clearly some problem at small angles - perhaps due to a problem with the PSF, perhaps it's saying something about where the dust is.  This is the point where things get interesting.

Technical: Putting it all together: More on the flux issue

As mentioned previously, my best fit to the HETG flux Fx(1.5-9.5 keV) = 3.7 ph/cm^2/s

The HRC response was calculated using weights from 0.5-9.8 keV, although it drops off to < 5% for E > 5.5 keV.  The flux below 1.5 keV is negligible.

Of course, this value includes the scattered flux since the HETG data are in CC mode; in reality, only 80-85% of this is truly from the source, the other 15-20% is scattered and not seen directly. So we add a correction to make it just be the true direct source rate.  I get this 15-20% value from a halo calculation which gives 5% scattering for 1e22 cm^-2 worth of MRN77-type dust, so for 3-4 x10^22 (exact value unknown) it's 15-20%.  I get slightly lower values for some of the reasonable ZDA04 values.  Taking 15% as 'typical', though, this 3.7 ph/cm^2/s becomes 3.1 ph/cm^2/s of direct flux.

Now, within 2" of the GX5-1, there are 500795 cts.  The 'ONTIME' is 4785 s, for a total count rate is 104.6 cts/s, and this should correspond to 90% of the total direct counts, or to 116.3 cts/s total.  Note that I have not taken deadtime into account yet. The effective area here is 62.2 cm^2 (taken from the exposure map), so that means the flux required to get the total observed count rate is 116.3/62.2 = 1.87 ph/cm^2/s.   This is an OK match with the HETG results.  If the deadtime is 30.5%, that can be treated as either reducing the observation time (which increases the measured count rate). Alternatively, if we wish to compare to the 1.87 ph/cm^2/s, we can simply decreasing the estimated flux we got from other sources.  This then converts 3.1 ph/cm^2/s to 2.1 ph/cm^2/s, which is 12%  higher than the 1.87 value.  That's within the range I'm expecting between the  HETG flux and the HRC flux based on the RXTE ASM data.


So, I will use 1.87 ph/cm^2/s as the 'flux' of my HRC-I source, since I only care about the post-deadtime-corrected value, not the actual flux.

Pulling it all together - the flux

Next up is to normalize the surface brightness by the source flux. This has already been done for the 3C273 observation I'm using for the PSF, but remains to be done for the  GX5-1 observation.  The tricky bit is that the 'source flux' depends both on the time the observation was done (as discussed earlier) and on how it's measured.  Chandra can measure the flux from a point source quite easily and accurately, but RXTE is a non-imaging detector and so a flux measured with RXTE will include both the point source and a significant amount of the 'scattered' flux.  Bringing these two measurements into alignment requires a determination of how much flux is scattered - which is, in fact, what our ultimate goal is.

Comparing previous results, Smith et al. (2006) found Fx(1-10 keV) = (5.2±0.1)x10^-8 erg/cm^2/s for their Chandra/ACIS observation of GX5-1 (ObsID 109), while Ueda et al (2005) found Fx(1-10 keV)=4.3x10^-8 erg/cm^2/s from their HETG data (ObsID 716).   From the fit I did of the new HETG spectrum I got Fx(1-10 keV)=3.7x10^-8 erg/cm^2/s, about 40% less than the ACIS observations and 16% less than the HETG flux.  This variation can be compared to the RXTE all-sky-monitor count rates:

ObsID   Description    1.3-3 keV     5-12.1 keV
                                     cts/s              cts/s
109      ACIS-S            13.4±1.0      35.1±1.3
716      1st HETG        14.2±1.1      32.3±1.3
5888    2nd HETG       12.3±1.3      40.8±2.2
7029    HRC-I              11.2±3.4      31.5±6.2

Of course, this emission includes the direct and scattered flux, so it's not possible to convert these into fluxes directly (using, say, WebPIMMS) and get a sensible comparison.  But the range of count rates here is 25% or less, which for an X-ray binary is practically constant.  It also suggests that while the flux from the 2nd HETG observation might be slighter larger than during the HRC-I observation, the difference is less than 25%.

In fact, plugging the total ASM count rate for the HRC-I observation (63.38 cts/s) into WebPIMMS, along with kT=10.5 keV, NH=4e22 cm^-2 (from Smith et al. 2006), gives an HRC-I count rate (including all scattered photons) of 131 cts/s; adding 55 cts/s for the background gives 186 cts/s, just the telemetry limit.  This is reasonably close to the estimate I got from the HETG fit, although keep in mind that was just the source itself and this is the source + scattered flux so I would have expected it to be larger.  Plus the deadtime measurement from the data itself of 33% suggests this is an underestimate since a true count rate of 186 cts/s should have a lower deadtime than 33%.  Hmm.  Have to think more about this.

Pulling it all together - the Point Spread Function

At this point I have a radial profile for the source, in physical units - ie, units that can be compared to theoretical predictions.  However, there are three components to the surface brightness at any point: (1) the point spread function (PSF) of the mirror, (2) the diffuse background, and (3) the actual scattered X-rays that I care about.  I need to come up with some way to estimate the value of (1) and (2).  As it turns out, the diffuse background (term #2) is basically flat, so I can fit that pretty easily by just adding in a constant term.  However, the PSF (term #1) is energy-sensitive and position-dependent.

There are a number of ways to get Chandra's PSF.  The calibration database contains images of point sources that show the PSF, but these only have information near the source, out to ~10'' or so.  There is also the Chandra Ray Trace code (ChaRT) that can model the Chandra PSF for any set of energies.  The problem there is that ChaRT (also called SAOsac) is tuned up using the same calibration database images, and so it tends to give the wrong answers when looking more than 10'' away from the source.  The problem here is that I need the PSF from about 2'' to about 100'' or even 1000'' away from the source, and the models just won't do that.  I've included Figure 4 from Smith, Edgar & Shafer (2002) to demonstrate the problem.  This figure compares observations of Her X-1, a nearly-perfect point source compared to the best-possible ChaRT/SAOsac model.  Although near the source all is well, even by 30'' - 40'' away there are problems.

The solution is to use long observations of bright point sources, which inherently measure the PSF.  In this case, the quasar 3C273 observed with the HRC-I (ObsID 461) is about as good as it gets.  I used this source in Smith (2008) when measuring the halo of GX13+1, and it worked reasonably well.  The biggest problem is that the PSF is energy dependent, and the spectrum of 3C273 is not the same as GX13+1 or GX5-1 for that matter.

However, the energy-dependence isn't that great, and there isn't a better choice, so one works with what one can. Of course, trusting statements like that is the road to disaster. So I've attached two figures that prove it.  The first is the PSF measured from Her X-1 (and modelled with SAOsac) fit to a power-law for a range of energies.  The units in the y-axis are 'arcmin^-2', which are really the surface brightness in units of photons/cm^2/s/arcmin^2 divided by the source flux, in units of photons/cm^2/s. This gives the rather strange, but handy, units of arcmin^-2.  This figure shows that the PSF really is pretty well fit with a power law - and that the ChaRT/SAOsac fits are awful far off-axis.  The second figure simply plots the actual fit parameters for the power-law fits. This plot shows that the power-law slope is roughly constant around alpha=1.9, while the amplitude of the PSF increases almost linearly with energy.  The 'bobble' around 2 keV is almost certainly due to the Silicon K edge at that energy, since both the mirrors and the detectors include lots of silicon.  The upshot of all this is that I feel confident that the biggest problem with using 3C273 as a surrogate for the PSF of GX5-1 (and GX9+1 for that matter) will be a slight offset up or down in the total PSF power.  For example, if 3C273 has more emission at high energies than GX5-1, then the predicted PSF from 3C273 will be slightly larger than the 'real' GX5-1 PSF, since the power-law fits show that the PSF gets more intense at higher energies.

This opens up the question: If the power-law fits work so well, why not just use them?  Well, I do, for ACIS halo measurements.  But these energy-dependent fits can only be made using an energy-sensitive detector like the ACIS.  However, ACIS suffers from pileup when observing bright sources, making measurements of the near-source PSF impossible.  Measurements in the 10''-30'' range can only be done well with the HRC, which has no energy resolution.  Fits to the PSF of 3C273 show that this simple power-law fit don't work so well, as this figure shows.  This fit uses a Gaussian model for the core of the PSF, and then two power-laws, one with a slope of 4 (primarily for the 1-10'' region) and a second with a slope of 2.4 for the 10-100''.  There is also a constant term to handle the background.  And, even with all of these terms, the fit isn't all that good at 10''.  So, I think it's best to use real observations rather than fits at this point.

Getting the radial profile

This step is simple in concept, and a pain in execution unless you've got a script that does it all.  Having worked on these halos for a while (check out my CV), I now have a number of these scripts.  Effectively, what one does is determine how many annuli one wants, and then determines the number of counts in each annulus (from the event file) and the effective area as well (from the exposure map).  Dividing the two gives the surface brightness in units of photons/cm^2/s/pixel, which can be converted into any other convenient units later.  Here are the commands I use to make the annulus file:

sherpa
rad = 10^[-0.3:3.0:0.03]
fp = fopen("annuli.dat","w");
() = fprintf(fp,"%f %f\n",0.0,rad[0]);
for (i=0;i
()=fclose(fp);

And then I just use a script I wrote, make_radprof_hrc.sl, to calculate the actual surface brightness:

make_radprof_hrc.sl hrc_evt2.fits srcfree.reg ../expmap 270.28383 -25.07775 annuli.dat

This script assumes that the exposure map is named hrc_expmap.fits.  The central Right Ascension and Declination are given in the command line as 270.28383 and -25.07775; these values were obtained from the event file itself, using ds9.  Note that this position is close to, but not exactly at, the simbad position for GX5-1 (270.284167, -25.079167).  The separation is 5.2'' - which is rather large for Chandra.

This deserves a closer look to see if there is some reason the positioning of the source would be incorrect.  Normally, Chandra positions are better than 1'', so a 5+'' offset is odd.  Now, it's possible simbad is wrong - remember that Chandra has the highest angular resolution of any X-ray satellite, and that GX5-1 is first and foremost an X-ray source.  Fortunately, simbad lists sources, and in particular lists this Ebisawa et al (2003) paper as the source.  Checking this paper, we find that they reference Liu et al (2001) as their source for the position, with the source position accuracy given as 3'' with a caveat that this is only an estimate.  However, Liu et al. reference the paper "Infrared observations of Galactic bulge X-ray sources" (Hertz & Grindlay 1984), which sounds as though it should have quite accurate positions.  Nonetheless, a quick scan of Hertz & Grindlay shows that they got their positions from Einstein observations, and the infrared observations in fact did not detect GX5-1 at all. 

Upshot: The Chandra position is almost certainly superior to the position listed in Simbad.  You can't trust everything on the web, after all...