Cas A Supernova Remnant:
data reduction and analysis

Supernova remnants are extended sources that are very bright in the x-rays. This manual takes you through the data reduction and spectral analysis of a supernova remnant using data from two different satellites.

The data reduction is done in xselect. You will look at the observation, define regions and extract spectra.

Ftools are used to create files which can convert the spectra you extracted in terms of the pulse height of an event into the implied energy of the incoming photon.

Finally you will look at the spectra in XSPEC, and model it with various classic models.

Cas A, Rosat analysis

• Go to the HEASARC w3-browse page. Do a Basic search for Cas A in Rosat.
• Take the PSPCB observation not the PSPCC observation. You will need just the basic science events.
• Extract the files using tar $-$xvf. Cd into the subdirectory where it put the data. Uncompress the data files using gunzip and begin xselect, as you have done before.

• To define regions and look at the spectra enter xselect read in events and extract an image in sky coordinates:
  xsel:ASCA $>$ read events rp500125n00_bas.fits
  Got new mission: ROSAT
  $>$ Reset the mission ? $>$[yes] yes
  xsel:ROSAT$-$PSPC $>$ set image sky
  xsel:ROSAT$-$PSPC $>$ extract image
  xsel:ROSAT$-$PSPC $>$ plot image

• NOTE: the other choice of coordinate systems is detector coordinates. You want to specify sky coordinates for ROSAT because 1) the position of objects on the sky is what is relevant to your analysis, and 2) ROSAT wobbles, which distorts the image in detector coordinates. The transformation of each photon event to sky coordinates accounts for this wobble, and the spatially dependent response functions which I will talk about later are correctly made from the sky coordinates for ROSAT. (Other detectors, in particular those on board ASCA, need regions to be defined in detector coordinates, so this is an important issue to keep in mind.)

• When you plot an image an saoimage window will appear. There are buttons across the top of the image: Scale, Color, Cursor, Pan, etc. Clicking on these will change the lower line of buttons to specific related tasks.

• Click on the ``Color'' button, then click on the ``cmap'' button to get various colormaps that will make the remnant easier to see than the gray-scale does. Pick a map by clicking on its button, (write and read obviously are not maps, everything else is.) I think the best one to visualize the details in an image is ``A'', but use whatever you like.

• Click on the ``Scale'' button. Try all the scaling factors on saoimage, to see what the entire shape of the remnant is like. Is it just the small compact source in the center and the rest is scattered emission?

• Use the ``Pan'' menu to zoom in on the remnant.

• Click on the ``Cursor'' button to define regions. Click on the button with a circle on it to define circular regions.
• Use the left mouse button to define the center of the circle, use the ``middle'' mouse button to tell saoimage where the edge of the circle should be. The equivalent to a middle mouse button is holding both buttons down at the same time. Save a region by typing ``s'' in the window. Exclude a region by typing ``e'' in the window. To make an annulus save the outer circle then exclude the inner circle.

• Define three annuli, one inner one, one of the bright shell, and one just outside the bright shell. Note: You need to exclude the inner regions after saving the outer regions.

• To save a region file click on the ``region'' button in the ``cursor'' menu, and then click the ``write'' button. You can also read in previously made region files using the read button. Toggle back to the region defining cursor buttons by hitting ``region'' again.

• I usually look at my regions in emacs to make sure they match up like I want them to. You also want the regions defined so that the program can make reasonable decisions on which pixels to include and which to exclude. For this you want the center position and the radii of the regions to be half-integer values. You also want the annuli to be concentric, and the outer circle of the innermost region to line up with the inner circle of the shell region, and the outer circle of the shell region to line up with the innermost circle of the outer region.

• So, look at all the regions, with ``more'' or ``emacs''. For circular regions the first two coordinates are the X,Y position of the center, and the third coordinate is the radius of the circle.

• Define the center of all the circles by the closest half integer to the position you picked when defining the shell region. Use emacs to make them all the same.

• Define the radii of the circles from the inside out. So change the innermost region's radius to the closest half-integer then use this radius for the excluded region in the shell region. Likewise change the outer radius of the shell region to the closest half-integer, and use this radius for the innermost circle in the outer region. Lastly change the outermost radius in the outer region to the nearest half-integer.

• Note: this prescription is just so that all the annuli line up, are concentric, and cover the entire surface of the remnant. Any prescription that does these three things is fine; you usually want to go back and define slightly different regions later anyway.

• To extract spectra from these regions, make sure you are in the right coordinate system (sky or detector), put in a filter region, and extract a spectrum in pi channels, like this:
  
  xsel:ROSAT$-$PSPC $>$ filter region inner.reg
  xsel:ROSAT$-$PSPC $>$ set phaname pi
  xsel:ROSAT$-$PSPC $>$ extract spectrum
  xsel:ROSAT$-$PSPC $>$ save spectrum cas_inner.pi
  $>$ Group ( or rebin ) the spectra before outputting? $>$[yes] no
  xsel:ROSAT$-$PSPC $>$ clear events
  xsel:ROSAT$-$PSPC $>$ clear region
  xsel:ROSAT$-$PSPC $>$ filter region shell.reg
  xsel:ROSAT$-$PSPC $>$ extract spectrum
  xsel:ROSAT$-$PSPC $>$ save spectrum cas_shell.pi
  $>$ Group ( or rebin ) the spectra before outputting? $>$[no] no
  xsel:ROSAT$-$PSPC $>$ clear events
  xsel:ROSAT$-$PSPC $>$ clear region
  xsel:ROSAT$-$PSPC $>$ filter region outer.reg
  xsel:ROSAT$-$PSPC $>$ extract spectrum
  xsel:ROSAT$-$PSPC $>$ save spectrum cas_outer.pi
  $>$ Group ( or rebin ) the spectra before outputting? $>$[no] no
• Note about pi channels: Pha stands for pulse height, pi stands for pulse-height invariant channels. The pha column is what is filled by the detector on-board ship, the pi column is a conversion of the pha data either on board or at the flight center that is normalized in some way. The response functions are generally made for converting the pi counts into energies, so this is the mode we want to extract the spectrum in.

• Now you must make files that enable you to convert the counts in the channels to incoming photons of different energies and allow for the efficiency of the detector at that energy, the scatter of the photons which is related to the point-spread function, and the efficiency of the detector as a function of position on the detector which also varies with energy.
• There are two files related to this, the response matrix function, RMF, and the ancillary response function, arf. The RMF gives the efficiency of the detector as a function of energy and is provided by HEASARC at:
  ftp://legacy.gsfc.nasa.gov/caldb/data/rosat/pspc/cpf/matrices/
  you want: pspcb_gain2_256.rmf
  $-$ valid for PSPCB data taken AFTER 1991 Oct 14
• You must make an ARF file for each of the regions you defined. The ftool that does this is pcarf. It reuires two calibration files to run, the RMF file and the crffile, spec_rsp from
  ftp://legacy.gsfc.nasa.gov/caldb/data/rosat/pspc/cpf/pspcb_v2.spec_resp

• After obtaining these two things there is an important change you need to make to your spectral files. This is necessary because of an unfixed glitch in the processsing. You need to change the CDELT1 and CDELT2 keywords to 2.595021e$-$4 times whatever the value of your WMREBIN keyword is (usually 15). You can use FKEYPRINT to find the value of WMREBIN and FPARKEY to put the new value into CDELT1 and CDELT2.
  
  fkeyprint cas_inner.pi,cas_shell.pi,cas_outer.pi WMREBIN fparkey 3.8925315E$-$03 cas_inner.pi+0 CDELT1
  
  where 3.8925315E$-$03 equals 15 times 2.595021e$-$4, cas_inner.pi is one of the spectral files and 0 is the extension number in which CDELT1 and CDELT2 reside.

  Do this for both CDELT1 and CDELT2 for all of your spectral files.

• Now you are ready to make arf's for your regions using pcarf. At your prompt type a command of the form:
  
  pcarf phafile rmffile outfile crffile
  example:
  pcarf cas_inner.pi pspcb_gain2_256.rmf cas_inner.arf pspcb_v2.spec_resp
• Now you are ready to put these files in XSPEC and start analyzing them. Please refer to the brief XSPEC guide at the end of this document for assistance.

• Enter xspec by typing xspec at your prompt

• Enter the data and response files and ignore the channels that are flagged as bad like so:
  XSPEC$>$ data 1:1 cas_inner.pi
  XSPEC$>$ resp 1:1 pspcb_gain2_256.rmf
  XSPEC$>$ arf 1:1 cas_inner.arf
  XSPEC$>$ data 2:2 cas_shell.pi
  XSPEC$>$ resp 2:2 pspcb_gain2_256.rmf
  XSPEC$>$ arf 2:2 cas_shell.arf
  XSPEC$>$ ig bad
• Note: the 1:1 and 2:2 in these commands are as follows: the first number is the data-group number, and the second number is the datafile number. Suppose you wanted to enter three data-files that were all going to have exactly the same model and two more that would have a different model. You would then designate them like this:
  XSPEC$>$ data 1:1 filename1
  XSPEC$>$ data 1:2 filename2
  XSPEC$>$ data 1:3 filename3
  XSPEC$>$ data 2:4 filename4
  XSPEC$>$ data 2:5 filename5
  The response, arf and background files (if there is a background file) only absolutely need the datafile number, but it doesn't harm anything to include both the group number and file number.

• Set the plot device to an xwindow: cpd /xw
• Plot by energy rather than channels: setplot energy
• Plot the data on a logarithmic scale: plot ldata

• Try to fit various basic spectral models to the continuum emission in this spectrum; bremsstrahlung, powerlaw, blackbody. You will also want an absorbing column of hydrogen between you and the source, which modifies these basic models by a multiplicative factor. Enter in a model like so:
  
  mo wabs(bbody)
  or
  mo wabs(brems)
  or
  mo wabs(power)
  Just press return through the section where it asks you for starting values of the parameters.

• Try fitting both datasets with one model to start with, but allow the normalizations to be separately free. Do this with a newpar command:
  newpar 6 1.0
  where 6 is the parameter number of the normalization and 1.0 is a nominal starting point. (Note that 6 may not be the parameter number of your normalization).

• Fit these models to the data using few iterations and not too much precision to get the fit started:
  fit 10 10.0
  where 10 is the number of iterations and 10.0 is the reduction in chi-squared required from one iteration to the next. The reasons you want to start with a low level of precision and work up is that you have a large number of data-points so you want the fit to go quickly, and you don't want the fitting program to get stuck in some local minimum of chi-squared. You can increase the level of precision by entering a smaller number down to about 1E-6, but it is unwise to raise the number of iterations, at least for the purposes of preliminary analysis.

• Try all of these models first. See where in the spectrum they fit well to the data and where they don't.

• Try modeling both datasets with one model normalizing them with respect to their areas. Just use the radii of the regions you defined to determine the areas, for example:
  inner: radius r=7.25 area proportional to r$^{2}$ = 52.5625
  shell: outer radius R= 24.00 area proportional to (R$^{2}$$-$r$^{2}$) = 523.4375
  normalization of shell = normalization of inner times (R$^{2}$$-$r$^{2}$) divided by r$^{2}$ = 9.9583829

• To tie the normalizations back together again with a factor use:
  newpar 6 = 3 9.9583829
  where 3 is the normalization of the first data-group and 6 is the normalization of the second data-group and 9.9583829 is the appropriate area scaled factor. You may or may not find that this is the appropriate normalization factor. If not then it means that one of the regions is emitting more photons per square cm in our direction than the other region is. Think about why this might be, why is the canonical supernova remnant shell-like?

• You can save the data files you have used or the model you have fitted using the save command in xspec. Refer to the appendix of this document for a description of save, or use the help command within xspec.

ASCA Cas A analysis

• For the sake of brevity I have extracted the ASCA spectrum, background files, RMF and ARF's for you. The process is similar to that for ROSAT but not the same. For extensive detail on how to do ASCA extraction consult the ABC guide at: http://heasarc.gsfc.nasa.gov/docs/asca/abc/abc.html

• The imaging capabilities of ASCA are not as good as that of ROSAT, so I have just extracted one region over the whole remnant. The real value of ASCA is its spectral range and resolution, as you will see.

• Enter xspec by typing xspec at your prompt

• Enter the data, background and response files using the data file script that I made in xspec:
  XSPEC$>$ @snr_all_data.xcm

• Set the plot device to an xwindow: cpd /xw
• Plot by energy rather than channels: setplot energy
• Plot the data on a logarithmic scale: plot ldata

• Notice the bumpy structure of the spectrum. From left to right, starting with the large bump just below 2keV we have the emission lines of Si, S, Ar, Ca, and the large bump on the right is Fe. Just seeing these emission lines tells you that there is hot thermal, highly ionized gas producing this emission.

• Try to fit various basic spectral models to the continuum emission in this spectrum; bremsstrahlung, powerlaw, blackbody. You will also want an absorbing column of hydrogen between you and the source, which modifies these basic models by a multiplicative factor. Enter in a model like so:
  
  mo wabs(bbody)
  or
  mo wabs(brems)
  or
  mo wabs(power)
  Just press return through the section where it asks you for starting values of the parameters.
• Do a show par. Given in the show par along the top are the model components with component numbers in parenthesis. The parameters of each component are then listed. The leftmost number is the parameter number, use this for newpar commands. The second number is not useful, it describes, out of the subset of parameters that are separate degrees of freedom, what number parameter is this, just disregard it. The third number tells you which model number is this parameter a part of. Then there is a description of what the parameter is, then its current value, then the estimated error on that value or frozen if the parameter is frozen.

• Plot the model against the data using plot ldata. When in doubt always look at the data.

• Fit these models to the data using few iterations and not too much precision to get the fit started:
  fit 10 100.0
  where 10 is the number of iterations and 100.0 is the reduction in chi-squared required from one iteration to the next. Later you will want to increase the precision of the fit. You could also have just tried to renormalize your model with a renorm command, and this is useful if the model is simply incredibly way off to begin with but once the fit is in the same ballpark as the data, it's better to do a fit.

• Try all of these models first. See where in the spectrum they fit well to the data and where they don't.

• With a wabs(bbody) model in place try adding in gaussian lines for the ion emissions.

• Add in a newcomponent model between the 1st and 2nd models like this:
  XSPEC$>$ addcomp 2 gaus
  Input parameter value, delta, min, bot, top, and max values for ...
  Mod parameter 2 of component 2 gaussian LineE keV
  6.500 5.0000E$-$02 0. 0. 1.0000E+06 1.0000E+06
  you enter 1.9
  Mod parameter 3 of component 2 gaussian Sigma keV
  0.1000 5.0000E$-$02 0. 0. 10.00 20.00
  you enter 0.01
  Mod parameter 4 of component 2 gaussian norm
  1.000 1.0000E$-$02 0. 0. 1.0000E+24 1.0000E+24
  you enter 0.1
  Here you want to choose a Line Energy (the first parameter) that corresponds approximately to one that you see on the graph of the spectrum. The second parameter, the line width or Sigma, should be set to be 0.01 because this is the resolution appropriate for ASCA. The third parameter, the norm, should be set to 0.1 or less to start off with so that it doesn't totally dominate the spectrum, but it doesn't disappear either.

• Freeze the line width parameter, Sigma, because that is a mission dependent parameter for our case and we don't want it to vary.

• Note also, you should be careful where you add a new component to the model. If there is a multiplying model already present like wabs, you will probably want this to modify the new component you are adding so placing the new component between the multiplying model and the additive model (brems, pow, or bbody) is the safest way to do this. Addcomp 2 implies that the new model should be added between models 1 and 2.

• When you add a new component always plot the data again.
• You may want to tweak the parameters yourself before letting XSPEC do a fit. To do this use the newpar command. For example if you can see that the line energy you inputted is to the left of the actual emission line in the spectra, and a show par command gives you:
  11 4 5 gaussian LineE keV 1.750 +/$-$ 0.4215E$-$03
  for the gaussian in question, you can change the Line Energy manually by typing:
  newpar 11 1.8

• After entering a gaussian, freezing the Sigma at 0.01, plotting the data, and tweaking the parameters yourself a little, go ahead and fit it, say fit 10 100.0.

• Add in another gaussian line for the Fe line, at around 6.5 keV. Do all the same steps as for the Si line.

• If you are using a bbody model you'll see that even with the gaussian line, the fit at the higher energies is not that good. First try increasing the precision of the fitting mode, say fit 10 1.0.

• It doesn't help, right? So try a different continuum model. You can do this without taking away the gaussian lines by using a addcomp to put in the brems model followed by a delcomp to remove the bbody model. Make sure you've adding in the new component and taken the old one away before running a fit or renorm.

• After fitting it you should see that the bremss model does much better for the higher energy part of the spectrum.

• Play around with modelling the rest of the emission lines.

Appendix: Xspec command basics

• To enter XSPEC simply type xspec at your prompt. You will now be in xspec and get an XSPEC$>$ prompt.
• Basic commands in XSPEC
  • help When in doubt within XSPEC simply type help. This will give you a list of available commands. Typing the name of any of these commands will give an explanation of their usage. You can get more complete descriptions of all the commands I am about to list by using the help command. To get out of help simply keep pressing return.

  • exit exit xspec, quit also works.

  • data filename Use this command to enter the spectral file. Type data followed by the file name (or full location and name of the file if it is in a different directory than your current working directory)

  • backgrnd filename Specifies the background file for the data file.

  • response filename Either the full multiplied together rsp file or just the rmf file.

  • arf filename The arf file.

  • ig bad ignore channels that are flagged as bad.

  • ig 1:0.0$-$0.2 10.0$-$** ignore channels for dataset 1, in the given energy ranges, where the decimal place indicates that it is an energy range not a channel number and ** indicates to go to the edge of the existing data. This command can be used with channel numbers also which would be indicated by the lack of a decimal point in the range specification.

  • not notice channels or energy ranges. Same syntax as ignore.

  • cpd $\backslash$ xw Set plotting device to be an xwindow.

  • setplot energy Sets the coordinate system to the energy of the channel rather than just the number of the channel.

  • plot (ldata resid) Plot the data in the window opened by cpd. You can qualify how you would like the plot done by specifying ldata, which makes the y-axis logarithmic, or resid, which plots the residuals if there is already a spectral model in use.

  • model Enter in a model or combination of models, combined using (,), +, *, or /. Some models such as brem, gauss, powerlaw, (bremsstrahlung, gaussian, and powerlaw) are additive models, they are directly associated with the emission of photons. Some models, specifically wabs (an absorbing hydrogen column model), are multiplicative models, they modify existing models by an energy dependent factor. An example of entering in a spectral model would be: mo wabs(brems + gaus + gauss +gau) + wabs(pow)
    which is equivalent to
    model wabs*(bremss + gaussian + gaussian + gaussian) + wabs(powerlaw)
    because xspec has command completion and will just fill in the rest of the model name automatically without typing tab or anything. After you have typed this XSPEC will prompt you for values of the parameters to start fitting with. You can enter a value, make one parameter equal to a factor times a previous parameter by typing `` = 13 2.2'', if you want the current parameter to equal 2.2 times parameter 13, or just press return for the default value.

  • delcomp 3 Delete the third component in the model. When the parameters of a model are given, the model is described at the top of the display with numbers in brackets next to each model component.

  • addcomp 2 gauss Adds a gaussian component to the model between the current 1st and 2nd components. If there is a current multiplicative model modifying the additive models you should be careful to place the new model component such that it is within the parentheses if you want it to also be modified by this same multiplicative model.

  • show (par, allfile) Show par lists the model parameters, show allfile shows information about the files such as what the file names are, the count rates and the noticed channels, show all gives all this information.

  • freeze 1$-$4,7,39 Freeze model parameters 1,2,3,4,7 and 39. The number of a model parameter is listed at the far left-hand side after a show par or a fit has been run.
  • thaw thaw model parameters, same syntax as freeze. It is best to freeze all but the most important parameters first and then slowly thaw the rest of the parameters if the model is a complicated one.

  • fit 10 0.01 Fit the model to the data, using 10 iterations and requiring a 0.01 reduction in chi-squared between iterations or else a fit is found. Both the number of iterations and the level of precision are optional parameters, although if you want to change the level of precision you must also type in the number of iterations. The default of the number of iterations and precision level is their most recent value.

  • newpar 6 11.0 change the value of parameter 6 to the value of 11.0. If you want to set one parameter equal to a factor times another parameter use newpar 15 = 6 0.55, which would make parameter 15 equal to 0.55 times the value of parameter 6.

  • save model filename Saves the current model and the values of all its parameters as filename.xcm. You can re-enter a previously saved model by typing
    @filename.xcm
    at the xspec prompt.
  • save files 'filename' Saves the filenames of the data, background, response files, and the ignored channels as filename.xcm. You can re-enter the dataset you are working on after you've saved it by typing
    @filename.xcm
    at the xspec prompt.

• Word to the wise, or rather the speedy: Both xselect and xspec have tab completion, so if you start typing a filename you can just press tab and it will finish it for you.