在得到X波段图像之后,我们可以通过拟合光谱来测量天体温度、金属丰度等信息,就要用到XSPEC。作为X射线光谱分析的通用软件,已经被包含在NASA的高能物理套件HEASOFT之中。是X射线分析套件XANADU的主要组件。这个名字本来是元上都的英译,马可波罗曾在他的游记中极力描绘那里的奢华繁荣,后又经英国诗人柯勒律治题诗歌颂,而演变为东方仙境的代名词。
言归正传,下面我们采用上文得到的 MACS J0257.6-2209 定标文件来提取光谱。如果不想重复前面的步骤,只需提取压缩包primary目录下的evt2.fits 文件,执行下面的命令即可。
1 |
dmcopy "evt2.fits[events][energy=500:7000][bin x=2000:5500:2,y=2000:5500:2]" img_tot.fits |
因为星系团是展源,为了消除背景AGN的影响,我们先用wavdetect做小波探测找出图中的所有点源。
1 2 3 4 5 6 7 8 9 10 11 12 13 |
punlearn wavdetect pset wavdetect infile=img_tot.fits pset wavdetect outfile=sources_tot.fits pset wavdetect scellfile=cell_tot.fits pset wavdetect imagefile=image_tot.fits pset wavdetect defnbkgfile=background_tot.fits pset wavdetect expfile=none pset wavdetect scales="2 4 8 16" pset wavdetect clobber=yes pset wavdetect sigthresh = 1e-06 pset wavdetect regfile=sources_tot.reg pset wavdetect mode=h wavdetect |
可以直接在事件文件中查看找到了多少点源
1 |
dmlist sources_tot.fits blocks |
也可以直接在ds9中显示出来。
1 |
ds9 img_tot.fits -region sources_tot.reg& |

检查是否有误判的区域,确认无误后,就可以剔除点源来提取光谱了。
1 2 3 |
dmcopy "evt2.fits[exclude sky=region(sources_tot.reg)]" field_nopt.fits dmcopy "field_nopt.fits[events][energy=500:7000][bin x=2000:5500:2,y=2000:5500:2]" img_cluster.fits ds9 img_cluster.fits& |
在ds9中选出整个星系团。保存位置文件(cl.reg )。然后就可以提取光谱了。
1 2 |
dmcopy "field_nopt.fits[events][sky=region(cl.reg)]" evt2_cl.fits dmextract "evt2_cl.fits[events][bin pi]" cl.pha clobber=yes |
但是,我们选取的星系团区域中,并不全是来自星系团的光子。空间高能粒子,银河系X射线背景,无法分辨的遥远AGN都混杂其中。同可见光观测要做平场一样,我们还需要背景文件。获取背景文件有两种途径,一种是望远镜对空白天区的长时间曝光blank-sky background,但由于望远镜本身的状态变化,不同时期的背景文件并不通用,需要根据状态参数找到对应的修正文件。这个方法有些复杂,我们这里采用另外一种:利用CCD上的空白区域来提取背景。既然图像上的点源和星系团都已经挑出来了,剩下的就都是背景了。但选取时仍需要注意,尽量选取靠近星系团的区域,尽量使用目标所在的CCD,同时不要和图像边缘靠得太近。但如果星系团很大,占据了整个画面,就只有前一个方法可用了。
1 2 |
dmcopy "field_nopt.fits[events][sky=region(bck.reg)]" evt2_bck.fits dmextract "evt2_bck.fits[events][bin pi]" bck_cl.pha |
由于X射线望远镜对不同能段光子的有效探测面积和能量分辨率都不同,而且还随时间和位置变化。因此每次观测,我们都需要根据星系团所在位置计算响应文件。包括修正光子分布的Redistribution Matrix Function (RMF)和记录量子转换效率和有效面积的Ancillary Response File (ARF)。更具体的介绍可以参考MIT工作组的文档。
1 2 3 4 5 6 7 |
dmcopy infile="evt2_cl.fits[energy=300:2000][bin det=8]" outfile=cl_sou.wmap clobber=yes punlearn mkwarf mkwarf infile=cl_sou.wmap outfile=arf_cl.fits weightfile=cl_sou.wfef spectrumfile="none" egridspec=0.3:11.0:0.01 pbkfile=pbk0.fits clobber=yes punlearn mkacisrmf mkacisrmf infile=CALDB outfile=rmf_cl.fits wmap=cl_sou.wmap energy=0.3:11.0:0.01 channel=1:1024:1 chantype=PI gain=CALDB mode=h |
生成的响应文件可以通过rmfimg转换为图像来检查。
1 2 3 4 5 6 7 8 9 10 11 12 |
pset rmfimg infile = rmf_cl.fits pset rmfimg outfile = imgrmf_cl.fits pset rmfimg arf = arf_cl.fits pset rmfimg arfout = imgarf_cl.fits pset rmfimg product = no pset rmfimg verbose = 3 pset rmfimg clobber = yes pset rmfimg mode = h rmfimg ds9 imgrmf_cl.fits& ds9 imgarf_cl.fits& |

左侧为ARF右侧为RMF,横轴均为能道编号,纵轴为光子原始能量。对于ARF,由于硬X射线掠射角增加,望远镜有效面积减小,探测器量子效率也同时降低,因此呈现自下而上的递减趋势。而RMF在理想状态下应为狭窄的对角直线,每个光子都进相应的能道,但实现起来非常困难,只能根据所在能道反推原始能量分布。
现在我们已经准备好了全部文件,将所有这些文件同光谱关联到一起。
1 2 |
pset grppha clobber=yes grppha infile=cl.pha outfile=cl_grp.pha chatter=0 comm=' group min 1 & chkey RESPFILE rmf_cl.fits & chkey ANCRFILE arf_cl.fits & chkey BACKFILE bck_cl.pha & exit' |
在拟合光谱之前,还要考虑银河系晕中的中性氢消光,这与天体位置有关。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
$ dmlist evt2_cl.fits header | grep NOM 0055 RA_NOM 44.4640975717 Real8 Nominal RA 0056 DEC_NOM -22.1583582964 Real8 Nominal Dec 0057 ROLL_NOM 355.2881562242 Real8 Nominal Roll $ nh 2000 44.4640975717 -22.1583582964 >> Leiden/Argentine/Bonn (LAB) Survey of Galactic HI LII , BII 210.101242 -60.925470 Requested position at X and Y pixel 346.67 21.64 Search nH in 4 X 4 box Each pixel is 0.675 deg 0.675 deg nH calculated using all points within 1.0000 deg from input position RA DEC Dist nH 43.8807 -22.2531 0.5485 2.12E+20 43.8991 -21.7153 0.6856 2.41E+20 44.7581 -22.9792 0.8649 1.82E+20 44.7833 -22.4343 0.4044 1.97E+20 44.7998 -21.8885 0.4117 1.97E+20 LAB >> Average nH (cm**-2) 2.06E+20 LAB >> Weighted average nH (cm**-2) 2.06E+20 |
我们采用最后一个加权平均值。终于可以运行XSPEC了。
cpd /xw
#设置输出到屏幕,如果保存为文件则是 cpd file.ps /cps
data cl_grp.pha
# 载入光谱
show rate
# 显示计数率,曝光时间等信息
query yes
setplot energy
# 設定显示模式为能谱
abund angr
# 設定太阳金属丰度参考值,这里取Anders和Grevesse于1989发表的数据。详细列表。
notice **-**
ignore bad
ignore 0.0-0.5
ignore 8.0-**
# 显示光谱 0.5-8.0 keV 的部分
setplot rebin 3 20
# 显示时合并部分通道,方便查看。后面两个值分别是最小信噪比和最大bin size。
plot ldata
# 对数显示,就可以看到文章开头的光谱了。statistic cstat
# 采用 C-statistic 统计
model tbabs*mekal
# 定义模型,下面按提示输入参数。注意这里的中性氢密度nH 单位是10^22 atoms/cm^–2。详见Xspec用户手册。温度(2)、金属丰度(4) 以及归一化参数(7)是我们要拟合的,随便给个初值。红移在BAX星系团数据库中可以查到: 0.322(Simbad 和 NED中还没有)。
123456789101112131415 Input parameter value, delta, min, bot, top, and max values for ...1 0.001 0 0 100000 1e+061:TBabs:nH>0.02061 0.01 0.0808 0.0808 79.9 79.92:mekal:kT>51 -0.01 1e-06 1e-05 1e+19 1e+203:mekal:nH>0.0011 -0.01 0 0 1000 10004:mekal:Abundanc>0.30 -0.01 0 0 10 105:mekal:redshift>0.32216:mekal:switch>11 0.01 0 0 1e+24 1e+247:mekal:norm>0.0001freeze 1 3 5 6
thaw 2 4 7
#設定自由参数为2 4 7
renorm
fit 20 1e-2
#默认采用 Levenberg-Marquardt 算法拟合,需指定迭代次数,精度可选。
12345678910 Model TBabs <1> * mekal <2> Source No.: 1 Active/OnModel Model Component Parameter Unit Valuepar comp1 1 TBabs nH 10^22 2.06000E-02 frozen2 2 mekal kT keV 10.7532 +/- 0.9207163 2 mekal nH cm-3 1.00000E-03 frozen4 2 mekal Abundanc 0.438082 +/- 0.1290115 2 mekal redshift 0.322000 frozen6 2 mekal switch 1 frozen7 2 mekal norm 2.24587E-03 +/- 9.69028E-05steppar 2 0.1 10 20 4 0.1 0.6 20
# 如果为单参则直接给出一维分布
plot contour
# 显示结果如右图所示。可以看出温度在10keV左右误差较小,金属丰度不确定性较大。下面计算这两个参数的 $$1\sigma$$误差,虽然前面的拟合结果中也有,但并不可靠。如果没有steppar一步,算误差时还可能发现新的最优值,所以还是再算一下比较放心。
1234 error 1.0 2 4Parameter Confidence Range (1)2 9.84033 11.6813 (-0.912888,0.928126)4 0.312238 0.573214 (-0.125844,0.135133)setplot device spec.eps /cps
# 输出到指定文件
setplot rebin 3 20
#setplot 只改变显示效果,不影响实际数据,因此每次绘图都要重设
plot ldata res
# 增加残差部分
quit
这样我们就完成了X射线光谱拟合,得到了星系团MACS J0257.6-2209 的温度 T= 10.75 +/- 0.92keV 以及金属丰度 Z=0.44 +/- 0.13 。

参考资料:《Spectral fitting: XSPEC》2009