Xspec11 for beginners
From WikiVirgo
После получения спектров в рентгеновском диапазоне возникает желание вытянуть оттуда физические параметры, такие как поток излучения (в таком-то диапазоне), параметры моделей и т.д. Наиболее распространенным приложением, которое позволяет это сделать, является пакет xspec11, входящий в пакет программ по астрофизике высоких энергий HEASOFT.
В этой статье я попытаюсь последовательно изложить основные возможности пакета xspec11, не претендуя при этом ни на полноту изложения, ни на ее глубину.
xspec11 vs xspec
Хотя на данный момент существует xspec 12-той версии (он же просто xspec), он довольно существенно отличается от xspec 11-той версии (в частности, не все команды под 11-тый xspec адекватно воспринимаются 12-тым). По моим ощущениям, недостатков в 12-том xspec'е больше, а достоинств меньше, чем в 11-том. Поэтому я и использую 11-тый.
Начало работы
Запуск этого приложения происходит с помощью команды
$ xspec11
Перед запуском xspec11 надо инициализировать соответствующие переменные. В ВИРГО классах это обычно происходит путем набора в рабочем терминале
$ source /virgo/scripts/login.sh
в случае bash
$ source /virgo/scripts/login
если вы используете tcsh
Обычно также установен алиас - "virgo"
Подключение спектров
Поскольку xspec11 работает со спектрами, ему надо эти спектры показать. Наиболее простой метод это сделать - показать PHA спектры, полученные ранее командой grppha, к примеру
XSPEC> data 1:1 spectrum1.grp
Часто имеет смысл одновременно моделировать несколько спектров. В таком случае в xspec11 имеются две возможности. Если параметры модели спектра будут одни и те же (например, анализ совместных наблюдений двумя камерами), то второй спектр добавляется с помощью команды
XSPEC> data 1:2 spectrum2.grp
Как увидим далее, в плане моделирования этот случай является более простым. Если же параметры модели предполагаются разными, второй спектр добавляется с помощью команды
XSPEC> data 2:2 spectrum2.grp
Остальные спектры добавляются по такому же правилу.
Визуализация подключенных спектров
Для визуализации спектров необходимо, во-первых, перенаправить вывод на графический терминал,
XSPEC> cpd /xw
Во-вторых, надо построить спектр:
XSPEC> plot dataПолученный таким образом спектр будет не очень наглядным.
Во-первых, на оси абсцисс будет отложен номер канала, а не энергия фотона. Связь между первым и вторым, к сожалению, неоднозначна. Поэтому ее описывают с помощью responce matrix, которая для различных инструментов называется по-разному (например, для XMM-Newton/EPIC и Swift/XRT это RMF, Swift/BAT - DRM, RXTE/PCA - RSP). Для перехода в единицы энергии необходимо, чтобы соотв. responce matrix присутствовала в подключенном файле. Это условие является необходимым для дальнейшего анализа! Переключение оси абсцисс с каналов в кэВы производится с помощью команды
XSPEC> setplot energyПри этом энергия отображается в логарифмическом масштабе, поэтому каналы, начинающиеся с нуля энергий, будут отображены некорректно. К счастью, responce matrix достоверно известна только начиная с некоторой минимальной энергии. Для XMM-Newton/EPIC и Swift/XRT это, например, 0.3 кэВ. Поэтому для корректного отображения спектров (и для корректного дальнейшего анализа) необходимо выбросить все события с энергией меньше 0.3 кэВ:
XSPEC> ignore *:**-0.3
К слову, у responce matrix есть не только нижний, но и верхний предел. Для каждого конкретного инструмента он свой. Например, для MOS камер в XMM-Newton/EPIC или Swift/XRT он составляет 10 кэВ. Соотв. исключение выглядит так:
XSPEC> ignore *:**-0.3 10.0-**
Во-вторых, поскольку количество фотонов спадает с энергией, удобнее представлять ось ординат в логарифмическом масштабе:
XSPEC> plot ldataРазмерность по оси ординат - нормированные события в единицу энергии и единицу времени. Для перевода в единицы физического (т.е., независимого от внутренних свойств инструмента) потока (например, в keV cm − 1sec − 1keV − 1) необходимо знать эффективную прощадь, которая хранится в ARF файле. Важно, чтобы ARF файл был подключен. В частности, это делается внутри grppha.
Важно заметить, что физический поток не может быть изображен без процедуры моделирования спектра. К этому обстоятельству я вернусь позже.
В третьих, изображается не спектр источника, а результат вычитания фона (снова-таки, добавленного в grppha). Напомню, вычитание фона - отнюдь не безобидная процедура. В частности, она может быть неточной. При этом она тем менее точная, чем больше относительный вклад фона (особенно на высоких энергиях) и чем больше площадь регионов источника и фона (по сравнению с полем зрения инструмента). Кроме того, доминирование фона на высоких энергиях портит χ2 статистику, о чем я расскажу немного позже.
Моделирование спектров - рекомендации по выбору моделей
Для определения физических параметров необходимо спектры смоделировать. Моделирование в xspec11 производится с помощью счетного множества встроенных аддитивных, мультипликативных и конволютивных моделей и их комбинаций (а также еще одного счетного множества моделей, определяемых пользователем с помощью команды mdefine). Сама процедура выбора модели при этом часто находится на стыке науки и искусства. Ниже я перечислю несколько советов, используемых обычно при выборе модели.
- Во-первых, на низких энергиях спектр источника модифицируется поглощением, в основном в межзвездной среде. Есть несколько моделей межзвездного поглотителя. Я рекомендую использовать для этого модель phabs (предыдущая модель, wabs, несколько устарела).
Эта модель состоит из одного параметра - эквивалентной столбовой плотности (column density) атомов водорода. Данная величина может быть получена двумя путями. Для внегалактических источников подходят галактические обзоры. Я обычно использую данные обзора LAB, доступные здесь он-лайн (стоит отметить, что для ярких источников данные по столбовой плотности часто известны более точно; для этого надо просмотреть соотв. статьи). Для галактических источников данная величина является лишь ограничением сверху, а сама column density должна определяться из моделирования. Как вариант, возможна ситуация, когда кто-то сделал его до Вас, и Вы ему по каким-то причинам верите.
Кроме того, для особо нетривиальных объектов возможно наличие существенного поглощения внутри самой области излучения. В этом случае приходится либо использовать имеющиеся конволютивные модели, либо самим их изобретать.
- Во-вторых, необходимо проверить наличие тепловой компоненты в спектре. На нее непосредственно указывают эмиссионные линии в спектре объекта. В большинстве объектов видна, впрочем, и нетепловая компонента, поэтому излучение в общем случае может быть описано как сумма тепловой и нетепловой компонент.
- В-третьих, надо начинать с самых простых моделей, уточняя их по мере необходимости. Например, простейшая нетепловая модель - это powerlaw, описывающая степенной спектр (с фиксированным спектральным индексом). Следующие по сложности модели - это cutoffpl или bknpower - учитывают экспоненциальное обрезание и резкую перемену спектрального индекса, соотвественно. Простейшая тепловая модель - bbody - описывает излучение абсолютно черного тела. Эмиссионные линии могут быть описаны феноменологически, например, с помощью модели gaussian.
- В-четвертых, здесь имеется (пока незавершенная) классификация моделей, которые обычно используются для объектов определенного класса. Если Ваш объект принадлежит к известному классу, имеет смысл поискать модели среди перечисленных там. Если это не помогает, смело обращайтесь к первоисточнику
Моделирование спектров - пример источника нетеплового излучения
Теперь можно приступить к моделированию спектра. В качестве пробного объекта рассмотрим блазар Mrk 421 - яркий объект с отсутствием эмиссионных линий в рентгеновском диапазоне. Излучение данного объекта в интервале 0.3-10.0 кэВ достаточно точно описывается степенным законом.
Допустим, мы хотим смоделировать спектр с помощью модели, состоящей из нетеплового излучения. Если параметры модели ввести "по умолчанию", получится что-то вроде
XSPEC> model phabs*powerlaw Model: phabs<1>( powerlaw<2> ) Input parameter value, delta, min, bot, top, and max values for ... 1 0.001 0 0 1E+05 1E+06 1:phabs:nH> 1 0.01 -3 -2 9 10 2:powerlaw:PhoIndex> 1 0.01 0 0 1E+24 1E+24 3:powerlaw:norm> --------------------------------------------------------------------------- --------------------------------------------------------------------------- Model: phabs<1>( powerlaw<2> ) Model Fit Model Component Parameter Unit Value par par comp 1 1 1 phabs nH 10^22 1.00000 +/- 0.00000 2 2 2 powerlaw PhoIndex 1.00000 +/- 0.00000 3 3 2 powerlaw norm 1.00000 +/- 0.00000 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Chi-Squared = 6.0385754E+08 using 1697 PHA bins. Reduced chi-squared = 356468.4 for 1694 degrees of freedom Null hypothesis probability = 0.00
Поскольку объект внегалактический, мы можем определить column density (см. выше). Два ближайших измерения LAB (указанных здесь) дают значения column density: 1.88е20 и 1.94е20 cm − 2. С другой стороны, Mrk 421 - источник достаточно яркий, и столбовая плотность, определенная здесь по параметрам линии 21 см, оказалась равной 1.61е20 cm − 2. Для фиксации столбовой плотности наберем
XSPEC>newpar 1 1.61e-2 3 variable fit parameters Chi-Squared = 2.9676317E+08 using 1196 PHA bins. Reduced chi-squared = 248753.7 for 1193 degrees of freedom Null hypothesis probability = 0.00 XSPEC>freeze 1 Number of variable fit parameters = 2
Остальные параметры оставлены свободными. Следующим шагом будет моделирование:
XSPEC>query yes Querying disabled - assuming answer is yes
(этот параметр позволяет моделировать столько раз, сколько нужно для "схождения" величины χ2)
XSPEC>fit Chi-Squared Lvl Fit param # 2 3 427813. -3 2.796 0.1259 331095. -1 1.873 0.1976 15885.2 -2 2.135 0.1910 4464.42 -3 2.214 0.1901 4446.15 -4 2.215 0.1896 4446.14 -5 2.215 0.1896 --------------------------------------------------------------------------- Variances and Principal axes : 2 3 9.90E-07 | -1.00 -0.09 2.15E-08 | 0.09 -1.00 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Model: phabs<1>( powerlaw<2> ) Model Fit Model Component Parameter Unit Value par par comp 1 1 1 phabs nH 10^22 1.610000E-02 frozen 2 2 2 powerlaw PhoIndex 2.21541 +/- 0.990957E-03 3 3 2 powerlaw norm 0.189637 +/- 0.171366E-03 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Chi-Squared = 4446.144 using 1697 PHA bins. Reduced chi-squared = 2.623094 for 1695 degrees of freedom Null hypothesis probability = 0.00
Результат получился не очень достоверный ("плохой фит"). Для того, чтобы понять причину плохого фита, необходимо узнать, на каких интервалах энергий набирается χ2. Это можно сделать, набрав
XSPEC>plot ldata delchiОказалось, что основной вклад в χ2 дают бины с энергией менее 1 кэВ. Исключим временно фотоны с энергией меньше 0.7 кэВ из анализа и зафитируем результат:
XSPEC>ignore 0.0-0.7 ignoring channels 1 - 95 in dataset 1 Chi-Squared = 2811.878 using 1619 PHA bins. Reduced chi-squared = 1.738947 for 1617 degrees of freedom Null hypothesis probability = 0.00 XSPEC>fit Chi-Squared Lvl Fit param # 2 3 1625.81 -3 2.258 0.1977 1625.58 -4 2.258 0.1977 1625.58 2 2.258 0.1977 --------------------------------------------------------------------------- Variances and Principal axes : 2 3 is [[belongstoCategory::Software]] 1.87E-06 | -0.99 -0.14 2.78E-08 | 0.14 -0.99 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Model: phabs<1>( powerlaw<2> ) Model Fit Model Component Parameter Unit Value par par comp 1 1 1 phabs nH 10^22 1.610000E-02 frozen 2 2 2 powerlaw PhoIndex 2.25766 +/- 0.135639E-02 3 3 2 powerlaw norm 0.197721 +/- 0.249466E-03 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Chi-Squared = 1625.585 using 1619 PHA bins. Reduced chi-squared = 1.005309 for 1617 degrees of freedom Null hypothesis probability = 0.436
Как видно, результат значительно улучшился. Поэтому можно говорить с достоверностью, что спектр Mrk 421 на интервале 0.7-10.0 кэВ представляет "чистый" степенной закон со спектральным индексом 2.258. Замечу, что величину, стоящую за знаком "+/-", в общем случае нельзя считать 1σ ошибкой. Для определения более реалистичной 1σ ошибки необходимо набрать
XSPEC>error 1.0 2 Parameter Confidence Range ( 1.000) 2 2.25633 2.25900 ( -1.335382E-03, 1.342535E-03)
Замечу, что точность описания 1σ ошибки с помощью метода "+/-" сильно падает с ростом числа параметров.
Следующим шагом является анализ спектра на энергиях ниже 0.7 кэВ. Для этого подключим события с энергией 0.3-0.7 кэВ к уже смоделированному спектру и посмотрим на его относительное поведение относительно модели на низких энергиях:
XSPEC>notice 0.3-0.7 A total of 79 more channels will be noticed Net count rate (cts/s) for file 1 313.3 +/- 0.2860 ( 99.1% total) using response (RMF) file... ../epn_ti40_sdY9_v6.8.rmf using auxiliary (ARF) file... pn-timing.arf using background file... pn-timing_back.pi Chi-Squared = 7474.199 using 1698 PHA bins. Reduced chi-squared = 4.406957 for 1696 degrees of freedom Null hypothesis probability = 0.00 XSPEC>plot ratio
Относительное отношение спектра к степенной модели уменьшается практически монотонно на энергиях меньше 0.7 кэВ. Это является свидетельством "обрезания" спектра на низких энергиях. Поскольку причиной такого обрезания не может быть поглощение в Галактике, существует три возможных механизма объяснения поведения спектра на энергиях меньше 0.7 кэВ:
- поглощение в рентгене больше поглощения, полученного по данным радионаблюдений, поскольку рентген производится на более близких расстояниях к активному ядру галактики, чем радио;
- спектр излучения данного источника в рентгене начинает спадать на энергиях менее 0.7 кэВ, отклоняясь от степенного закона. Это является свидетельством наличия максимума излучения на еще более низких энергиях;
- совокупность первых двух механизмов.
Для проверки первой гипотезы столбовую плотность необходимо "отпустить", и зафитировать с учетом изменяющейся столбовой плотности (поскольку Mrk 421 находится довольно далеко, z=0.031, есть смысл модифицировать поглощение с учетом его красного смещения с помощью модели zphabs):
XSPEC>addcomp 2 zphabs Input parameter value, delta, min, bot, top, and max values for ... 1 0.001 0 0 1E+05 1E+06 2:zphabs:nH>0.01 0 -0.01 0 0 10 10 3:zphabs:redshift>0.031 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Model: phabs<1>*zphabs<2>( powerlaw<3> ) Model Fit Model Component Parameter Unit Value par par comp 1 1 1 phabs nH 10^22 1.610000E-02 frozen 2 4 2 zphabs nH 10^22 1.000000E-02 +/- 0.00000 3 5 2 zphabs redshift 3.100000E-02 frozen 4 2 3 powerlaw PhoIndex 2.21541 +/- 0.990956E-03 5 3 3 powerlaw norm 0.189637 +/- 0.171367E-03 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Chi-Squared = 4721.297 using 1697 PHA bins. Reduced chi-squared = 2.787070 for 1694 degrees of freedom Null hypothesis probability = 0.00 XSPEC>fit Chi-Squared Lvl Fit param # 2 3 4 1844.14 -3 2.290 0.2065 2.1802E-02 1838.13 -4 2.292 0.2071 2.2110E-02 1838.13 -5 2.292 0.2071 2.2104E-02 --------------------------------------------------------------------------- Variances and Principal axes : 2 3 4 1.97E-08 | -0.05 0.83 -0.55 3.71E-06 | 0.96 0.19 0.20 7.90E-08 | -0.27 0.52 0.81 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Model: phabs<1>*zphabs<2>( powerlaw<3> ) Model Fit Model Component Parameter Unit Value par par comp 1 1 1 phabs nH 10^22 1.610000E-02 frozen 2 4 2 zphabs nH 10^22 2.210404E-02 +/- 0.457016E-03 3 5 2 zphabs redshift 3.100000E-02 frozen 4 2 3 powerlaw PhoIndex 2.29153 +/- 0.185286E-02 5 3 3 powerlaw norm 0.207122 +/- 0.407236E-03 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Chi-Squared = 1838.131 using 1697 PHA bins. Reduced chi-squared = 1.085083 for 1694 degrees of freedom Null hypothesis probability = 7.745E-03
Для проверки второй гипотезы необходимо, наоборот, оставить столбовую плотность на уровне радиоданных, а изменить модель континуума. Одна из возможностей - т.наз. логарифмически-параболическая модель (примененная к Mrk 421, например, здесь):
XSPEC>mdefine logpar E**(-(a+b*log(E))) XSPEC>mo phabs*logpar Model: phabs<1>( logpar<2> ) Input parameter value, delta, min, bot, top, and max values for ... 1 0.001 0 0 1E+05 1E+06 1:phabs:nH>1.61e-02 1 0.1 1E-22 1E-22 1E+22 1E+22 2:logpar:a> 1 0.1 1E-22 1E-22 1E+22 1E+22 3:logpar:b> 1 0.01 0 0 1E+24 1E+24 4:logpar:norm> --------------------------------------------------------------------------- --------------------------------------------------------------------------- Model: phabs<1>( logpar<2> ) Model Fit Model Component Parameter Unit Value par par comp 1 1 1 phabs nH 10^22 1.610000E-02 +/- 0.00000 2 2 2 logpar a 1.00000 +/- 0.00000 3 3 2 logpar b 1.00000 +/- 0.00000 4 4 2 logpar norm 1.00000 +/- 0.00000 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Chi-Squared = 9.4665640E+07 using 1697 PHA bins. Reduced chi-squared = 55915.91 for 1693 degrees of freedom Null hypothesis probability = 0.00 XSPEC>freeze 1 Number of variable fit parameters = 3 XSPEC>renorm Chi-Squared = 300504.9 using 1697 PHA bins. Reduced chi-squared = 177.3937 for 1694 degrees of freedom Null hypothesis probability = 0.00 XSPEC>fit Chi-Squared Lvl Fit param # 2 3 4 33192.5 -2 2.327 8.1001E-02 0.1861 2434.59 -3 2.150 0.1211 0.1939 2197.20 -4 2.158 0.1265 0.1934 2196.89 -5 2.158 0.1275 0.1934 2196.89 1 2.158 0.1275 0.1934 --------------------------------------------------------------------------- Variances and Principal axes : 2 3 4 2.09E-08 | -0.09 -0.07 0.99 9.13E-07 | 0.89 0.44 0.11 9.56E-06 | 0.44 -0.90 -0.02 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Model: phabs<1>( logpar<2> ) Model Fit Model Component Parameter Unit Value par par comp 1 1 1 phabs nH 10^22 1.610000E-02 frozen 2 2 2 logpar a 2.15805 +/- 0.161322E-02 3 3 2 logpar b 0.127539 +/- 0.280236E-02 4 4 2 logpar norm 0.193424 +/- 0.191494E-03 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Chi-Squared = 2196.889 using 1697 PHA bins. Reduced chi-squared = 1.296865 for 1694 degrees of freedom Null hypothesis probability = 1.196E-15
Причина неудачи данной модели состоит в неправильном описании поведения спектра на высоких энергиях - в реальном спектре не видно смягчения спектра на высоких энергиях, предсказываемой данной моделью. Более подходящей оказывается модель Банда, grbm, используемая для анализа спектров гамма-всплесков:
Model: phabs<1>( grbm<2> ) Model Fit Model Component Parameter Unit Value par par comp 1 1 1 phabs nH 10^22 1.610000E-02 frozen 2 2 2 grbm alpha -1.54412 +/- 0.446323E-01 3 3 2 grbm beta -2.26433 +/- 0.160980E-02 4 4 2 grbm tem keV 1.40702 +/- 0.144323 5 5 2 grbm norm 3.308631E-04 +/- 0.175369E-03 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Chi-Squared = 1626.638 using 1697 PHA bins. Reduced chi-squared = 0.9608024 for 1693 degrees of freedom Null hypothesis probability = 0.874
Эта модель представляет собой плавный переход между двумя степенными законами с показателями степени − α и − β.
XSPEC>addcomp 2 zphabs Input parameter value, delta, min, bot, top, and max values for ... 1 0.001 0 0 1E+05 1E+06 2:zphabs:nH>1e-02 0 -0.01 0 0 10 10 3:zphabs:redshift>0.031 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Model: phabs<1>*zphabs<2>( grbm<3> ) Model Fit Model Component Parameter Unit Value par par comp 1 1 1 phabs nH 10^22 1.610000E-02 frozen 2 6 2 zphabs nH 10^22 1.000000E-02 +/- 0.00000 3 7 2 zphabs redshift 3.100000E-02 frozen 4 2 3 grbm alpha -1.54412 +/- 0.446323E-01 5 3 3 grbm beta -2.26433 +/- 0.160980E-02 6 4 3 grbm tem keV 1.40702 +/- 0.144323 7 5 3 grbm norm 3.308631E-04 +/- 0.175369E-03 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Chi-Squared = 3794.486 using 1697 PHA bins. Reduced chi-squared = 2.242604 for 1692 degrees of freedom Null hypothesis probability = 0.00 XSPEC>fit Chi-Squared Lvl Fit param # 2 3 4 5 6 1617.88 -1 -1.666 -2.270 1.654 1.7203E-04 4.8128E-03 1617.87 -1 -1.666 -2.270 1.654 1.7201E-04 4.7902E-03 --------------------------------------------------------------------------- Variances and Principal axes : 2 3 4 5 6 6.98E-14 | 0.00 0.00 0.00 1.00 0.00 5.17E-08 | 0.37 -0.18 0.22 0.00 0.89 1.03E-06 | -0.50 0.67 -0.35 0.00 0.43 3.47E-03 | 0.58 0.03 -0.82 0.00 -0.03 5.70E-06 | 0.53 0.72 0.41 0.00 -0.17 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Model: phabs<1>*zphabs<2>( grbm<3> ) Model Fit Model Component Parameter Unit Value par par comp 1 1 1 phabs nH 10^22 1.610000E-02 frozen 2 6 2 zphabs nH 10^22 4.790249E-03 +/- 0.209154E-02 3 7 2 zphabs redshift 3.100000E-02 frozen 4 2 3 grbm alpha -1.66562 +/- 0.339496E-01 5 3 3 grbm beta -2.26951 +/- 0.252971E-02 6 4 3 grbm tem keV 1.65363 +/- 0.480576E-01 7 5 3 grbm norm 1.720141E-04 +/- 0.576008E-04 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Chi-Squared = 1617.874 using 1697 PHA bins. Reduced chi-squared = 0.9561903 for 1692 degrees of freedom Null hypothesis probability = 0.900
Таким образом, указанная модель довольно неплохо описывает поведение на энергиях 0.3-10.0 кэВ. Впрочем, полученный результат не является окончательным. Дальнейшее исследование вопроса требует дополнительных данных (например, со спектрометров RGS), а также более точного понимания характеристик инструмента (RMF и ARF) на низких энергиях. Поэтому вопрос о статистической достоверности отклонения спектра от степенного на низких энергиях требует дополнительного исследования. Данная же задача носит учебный характер.
Моделирование спектров - пример источника теплового излучения, ч.1
Рассмотрим теперь особенности моделирования источника теплового излучения. Примером является одна из областей остатка сверхновой SN1006. Этот остаток сверхновых является "оболочечным" - основной вклад в его рентгеновское излучение приходит из области, близкой к границе расширяющегося остатка сверхновых. В данной области была выбрана часть, отвечающая за рентгеновское излучение, из которой получен MOS1/MOS2 спектр с помощью метода ESAS. Здесь я приведу пример моделирования данного спектра.
Во-первых, загрузим файлы спектра:
XSPEC>da 1:1 mos1S001-obj-reg00-grp.pi Net count rate (cts/s) for file 1 0.3195 +/- 5.5078E-03( 97.6% total) using response (RMF) file... mos1S001-reg00.rmf using auxiliary (ARF) file... mos1S001-reg00.arf using background file... mos1S001-back-reg00.pi 1 data set is in use XSPEC>da 1:2 mos2S002-obj-reg00-grp.pi Net count rate (cts/s) for file 2 0.3290 +/- 5.5811E-03( 97.1% total) using response (RMF) file... mos2S002-reg00.rmf using auxiliary (ARF) file... mos2S002-reg00.arf using background file... mos2S002-back-reg00.pi 2 data sets are in use XSPEC>setpl en XSPEC>ig bad ignoring channels 64 - 87 in dataset 1 ignoring channels 66 - 222 in dataset 2 XSPEC>ig 0.0-0.3 10.0-** ignoring channels 1 - 14 in dataset 1 ignoring channels 1 - 14 in dataset 2 ignoring channels 63 - 87 in dataset 1 ignoring channels 89 - 222 in dataset 2
Далее необходимо выбрать подходящую спектральную модель. Начнем, как обычно, с поглощения. Результаты обзора LAB дают значение галактической столбовой плотности возле SN1006 около 7.5
cm − 2, что является оценкой сверху для галактического объекта. С другой стороны, несмотря на то, что SN1006 - объект галактический, он находится довольно высоко над галактическим диском (на высоте около 550 пк), поэтому величина столбовой плотности газа перед SN1006 не должна сильно отличаться от галактической. Мы ограничим ее в области 6.0-7.5cm − 2.
Поскольку мы выбираем область, где тепловое излучение доминирует, необходимо выбрать модель, адекватно описывающую тепловое излучение. Такой моделью для узких "прифронтовых" участков вблизи оболочкоподобных остатков сверхновых является, в частности, модель vpshock. Зададим параметры модели:
XSPEC>mo phabs*vpshock Model: phabs<1>( vpshock<2> ) Input parameter value, delta, min, bot, top, and max values for ... 1 0.001 0 0 1E+05 1E+06 1:phabs:nH>0.065 0.001 0.06 0.06 0.075 0.075 1 0.01 0.0808 0.0808 79.9 79.9 2:vpshock:kT>1.0 0.01 0.1 0.1 3.0 3.0 1 -0.01 0 0 1 1 3:vpshock:H> 1 -0.01 0 0 1E+03 1E+04 4:vpshock:He> 1 -0.01 0 0 1E+03 1E+04 5:vpshock:C> 1 -0.01 0 0 1E+03 1E+04 6:vpshock:N> 1 -0.01 0 0 1E+03 1E+04 7:vpshock:O> 1 -0.01 0 0 1E+03 1E+04 8:vpshock:Ne> 1 -0.01 0 0 1E+03 1E+04 9:vpshock:Mg> 1 -0.01 0 0 1E+03 1E+04 10:vpshock:Si> 1 -0.01 0 0 1E+03 1E+04 11:vpshock:S> 1 -0.01 0 0 1E+03 1E+04 12:vpshock:Ar> 1 -0.01 0 0 1E+03 1E+04 13:vpshock:Ca> 1 -0.01 0 0 1E+03 1E+04 14:vpshock:Fe> 1 -0.01 0 0 1E+03 1E+04 15:vpshock:Ni> 0 -1E+08 0 0 5E+13 5E+13 16:vpshock:Tau_l> 1E+11 1E+08 1E+08 1E+08 5E+13 5E+13 17:vpshock:Tau_u> 0 -0.01 0 0 5 5 18:vpshock:Redshift> 1 0.01 0 0 1E+24 1E+24 19:vpshock:norm> --------------------------------------------------------------------------- --------------------------------------------------------------------------- Model: phabs<1>( vpshock<2> ) Model Fit Model Component Parameter Unit Value par par comp 1 1 1 phabs nH 10^22 6.500000E-02 +/- 0.00000 2 2 2 vpshock kT keV 1.00000 +/- 0.00000 3 3 2 vpshock H 1.00000 frozen 4 4 2 vpshock He 1.00000 frozen 5 5 2 vpshock C 1.00000 frozen 6 6 2 vpshock N 1.00000 frozen 7 7 2 vpshock O 1.00000 frozen 8 8 2 vpshock Ne 1.00000 frozen 9 9 2 vpshock Mg 1.00000 frozen 10 10 2 vpshock Si 1.00000 frozen 11 11 2 vpshock S 1.00000 frozen 12 12 2 vpshock Ar 1.00000 frozen 13 13 2 vpshock Ca 1.00000 frozen 14 14 2 vpshock Fe 1.00000 frozen 15 15 2 vpshock Ni 1.00000 frozen 16 16 2 vpshock Tau_l s/cm^3 0.00000 frozen 17 17 2 vpshock Tau_u s/cm^3 1.000000E+11 +/- 0.00000 18 18 2 vpshock Redshift 0.00000 frozen 19 19 2 vpshock norm 1.00000 +/- 0.00000 --------------------------------------------------------------------------- --------------------------------------------------------------------------- *WARNING*: ZNIBIN: some low energy bins before stored model Response starts at 3.00E-02 and stored model at 5.44E-02 Chi-Squared = 2.2663709E+11 using 99 PHA bins. Reduced chi-squared = 2.3856535E+09 for 95 degrees of freedom Null hypothesis probability = 0.00 XSPEC>reno Chi-Squared = 3778.220 using 99 PHA bins. Reduced chi-squared = 39.77073 for 95 degrees of freedom Null hypothesis probability = 0.00 XSPEC>query yes Querying disabled - assuming answer is yes XSPEC>fit Chi-Squared Lvl Fit param # 1 2 17 19 3199.96 -3 6.1076E-02 2.814 1.2205E+10 4.4084E-05 1314.10 -4 6.0072E-02 2.948 3.9159E+09 1.3663E-04 896.455 -5 6.0027E-02 2.655 4.1366E+09 2.9431E-04 687.758 -2 6.3286E-02 2.490 3.9427E+09 1.8163E-04 496.636 -3 6.0951E-02 2.627 4.0390E+09 2.6211E-04 424.141 -4 7.3372E-02 1.584 3.7351E+09 2.8832E-04 412.064 -1 7.3872E-02 2.023 3.5545E+09 2.6232E-04 400.075 -2 6.7806E-02 2.279 3.7303E+09 2.4611E-04 395.137 -2 6.2757E-02 2.410 3.8440E+09 2.3532E-04 392.294 -1 6.1153E-02 2.454 3.8890E+09 2.3673E-04 392.261 0 6.1161E-02 2.454 3.8883E+09 2.3611E-04 392.260 0 6.1149E-02 2.454 3.8888E+09 2.3610E-04 --------------------------------------------------------------------------- Variances and Principal axes : 1 2 17 19 4.08E-11 | 0.00 0.00 0.00 1.00 2.61E-05 | 1.00 0.02 0.00 0.00 7.95E-02 | -0.02 1.00 0.00 0.00 3.46E+16 | 0.00 0.00 -1.00 0.00 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Model: phabs<1>( vpshock<2> ) Model Fit Model Component Parameter Unit Value par par comp 1 1 1 phabs nH 10^22 6.114933E-02 +/- 0.931170E-02 2 2 2 vpshock kT keV 2.45427 +/- 0.325043 3 3 2 vpshock H 1.00000 frozen 4 4 2 vpshock He 1.00000 frozen 5 5 2 vpshock C 1.00000 frozen 6 6 2 vpshock N 1.00000 frozen 7 7 2 vpshock O 1.00000 frozen 8 8 2 vpshock Ne 1.00000 frozen 9 9 2 vpshock Mg 1.00000 frozen 10 10 2 vpshock Si 1.00000 frozen 11 11 2 vpshock S 1.00000 frozen 12 12 2 vpshock Ar 1.00000 frozen 13 13 2 vpshock Ca 1.00000 frozen 14 14 2 vpshock Fe 1.00000 frozen 15 15 2 vpshock Ni 1.00000 frozen 16 16 2 vpshock Tau_l s/cm^3 0.00000 frozen 17 17 2 vpshock Tau_u s/cm^3 3.888751E+09 +/- 0.185891E+09 18 18 2 vpshock Redshift 0.00000 frozen 19 19 2 vpshock norm 2.361034E-04 +/- 0.433997E-04 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Chi-Squared = 392.2598 using 99 PHA bins. Reduced chi-squared = 4.129050 for 95 degrees of freedom Null hypothesis probability = 9.254E-38
Качество фита, как видно, желает лучшего. В частности, имеются ярко выраженные излишки (residuals) на энергиях 0.3-0.4 кэВ, 0.7-0.8 кэВ, 1.4-2.1 кэВ, 8.0-10.0 кэВ. Причиной первого и последнего излишков является неполное вычитание фона (напомню, в методе ESAS вычитается только большая часть инструментального континуума, инструментальные линии, излишек мягких фотонов и космический фон должны моделироваться отдельно).
Поскольку данный пример я показываю скорее в образовательных целях, чем в научных (к тому же, для научных целей надо использовать несколько наблюдений, для получения достаточно достоверных утверждений), я не буду детализировать модель фона, а вместо этого выброшу данные на энергиях 0.3-0.4 кэВ и 8.0-10.0 кэВ из дальнейшего анализа. Следующим шагом является понимание того факта, что в остатке сверхновой содержания (abundances) некоторых элементов могут являться отличными от солнечного, и этот факт необходимо учитывать. В частности, это достоверно известно для магния и кремния (замечу, что в xspec11 по умолчанию используется таблица содержания солнечных элементов ):
XSPEC>err max 3. 9-10 Parameter Confidence Range ( 2.706) 9 3.92239 6.07017 ( -1.07831 , 1.06946 ) 10 5.63675 8.49970 ( -1.44834 , 1.41462 )
Тот факт, что содержание магния и кремния достоверно превышает солнечное в несколько раз, говорит в пользу того, что в тепловом излучении от области ударной волны доминирует излучение от выбросов сверхновой (ejecta). Поэтому для более реалистичного моделирования мы должны бы были "отпустить" большинство наблюдаемых содержаний элементов. Проблема в том, что таким образом мы добавляем к модели 11 (!) свободных параметров. Делать это в данном случае сильно рискованно, потому что, во-первых, у нас не очень большое количество спектральных бинов (после вычитания 0.3-0.4 кэВ и 8.0-10.0 кэВ областей их осталось 86), во-вторых, что самое главное, с каждым дополнительным независимым параметром xspec работает все труднее (в частности, сильно повышается вероятность χ2 попасть в локальный минимум). Поэтому количество свободных параметров необходимо пытаться уменьшать всеми возможными методами. Одним из таких методов является использование результатов теоретических предсказаний и численных моделирований содержаний тяжелых элементов. Не вдаваясь в подробности, приведу результат:
XSPEC>newpar 5-15 1 -0.01 0.1 0.1 10 10 5:vpshock:C>4.3 1 -0.01 0.1 0.1 10 10 6:vpshock:N>2.6 1 -0.01 0.1 0.1 10 10 7:vpshock:O>4.4 1 -0.01 0 0 1E+03 1E+04 8:vpshock:Ne>1.5 4.997 0.1 0.5 0.5 50 50 9:vpshock:Mg>15 7.063 0.1 0.5 0.5 50 50 10:vpshock:Si>49.9 1 -0.01 0.5 0.5 50 50 11:vpshock:S>= par10 Equating parameter :S to parameter :Si * 1 1 -0.01 0.1 0.1 10 10 12:vpshock:Ar>= par7 * 2 Equating parameter :Ar to parameter :O * 2 1 -0.01 0.1 0.1 10 10 13:vpshock:Ca>= par7 * 3 Equating parameter :Ca to parameter :O * 3 1 -0.01 0.1 0.1 10 10 14:vpshock:Fe>1 1 -0.01 0.1 0.1 10 10 15:vpshock:Ni>1 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Model: phabs<1>( vpshock<2> ) Model Fit Model Component Parameter Unit Value par par comp 1 1 1 phabs nH 10^22 6.000000E-02 +/- 0.283890E-01 2 2 2 vpshock kT keV 2.08885 +/- 0.314160 3 3 2 vpshock H 1.00000 frozen 4 4 2 vpshock He 1.00000 frozen 5 5 2 vpshock C 4.30000 frozen 6 6 2 vpshock N 2.60000 frozen 7 7 2 vpshock O 4.40000 frozen 8 16 2 vpshock Ne 1.50000 frozen 9 8 2 vpshock Mg 15.0000 +/- 0.00000 10 9 2 vpshock Si 49.9000 +/- 0.00000 11 9 2 vpshock S 49.9000 = par 10 12 7 2 vpshock Ar 8.80000 = par 7 * +2.00000 13 7 2 vpshock Ca 13.2000 = par 7 * +3.00000 14 10 2 vpshock Fe 1.00000 frozen 15 11 2 vpshock Ni 1.00000 frozen 16 12 2 vpshock Tau_l s/cm^3 0.00000 frozen 17 13 2 vpshock Tau_u s/cm^3 4.611996E+09 +/- 0.218822E+09 18 14 2 vpshock Redshift 0.00000 frozen 19 15 2 vpshock norm 2.002326E-04 +/- 0.528642E-04 --------------------------------------------------------------------------- --------------------------------------------------------------------------- 6 variable fit parameters Chi-Squared = 39012.12 using 86 PHA bins. Reduced chi-squared = 487.6516 for 80 degrees of freedom Null hypothesis probability = 0.00 XSPEC>freeze 5-15 Model parameter 5 is already frozen Model parameter 6 is already frozen Model parameter 7 is already frozen Model parameter 8 is already frozen Model parameter 11 is already frozen Model parameter 12 is already frozen Model parameter 13 is already frozen Model parameter 14 is already frozen Model parameter 15 is already frozen Number of variable fit parameters = 4
Моделирование спектров - пример источника теплового излучения, ч.2
Следующим шагом является понимание излишка на 0.7-0.75 кэВ. Согласно одной из гипотез, этот излишек является следствием т.наз. large "shoulder" of the He-like oxygen line complex. Не вдаваясь в подробности, отмечу, что его можно смоделировать тремя узкими гауссианами со связанными нормировками, как показано ниже:
XSPEC>addcomp 2 ga Input parameter value, delta, min, bot, top, and max values for ... 6.5 0.05 0 0 1E+06 1E+06 2:gaussian:LineE>0.73 0.1 0.05 0 0 10 20 3:gaussian:Sigma>0.0 1 0.01 0 0 1E+24 1E+24 4:gaussian:norm>1e-05 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Model: phabs<1>( gaussian<2> + vpshock<3> ) Model Fit Model Component Parameter Unit Value par par comp 1 1 1 phabs nH 10^22 7.499971E-02 +/- 0.267359E-01 2 17 2 gaussian LineE keV 0.730000 +/- 0.00000 3 18 2 gaussian Sigma keV 0.00000 +/- 0.00000 4 19 2 gaussian norm 1.000000E-05 +/- 0.00000 5 2 3 vpshock kT keV 0.983738 +/- 0.132246 6 3 3 vpshock H 1.00000 frozen 7 4 3 vpshock He 1.00000 frozen 8 5 3 vpshock C 4.30000 frozen 9 6 3 vpshock N 2.60000 frozen 10 7 3 vpshock O 4.40000 frozen 11 16 3 vpshock Ne 1.50000 frozen 12 8 3 vpshock Mg 15.0000 frozen 13 9 3 vpshock Si 49.9000 frozen 14 9 3 vpshock S 49.9000 = par 13 15 7 3 vpshock Ar 8.80000 = par 10 * +2.00000 16 7 3 vpshock Ca 13.2000 = par 10 * +3.00000 17 10 3 vpshock Fe 1.00000 frozen 18 11 3 vpshock Ni 1.00000 frozen 19 12 3 vpshock Tau_l s/cm^3 0.00000 frozen 20 13 3 vpshock Tau_u s/cm^3 1.132493E+10 +/- 0.147843E+10 21 14 3 vpshock Redshift 0.00000 frozen 22 15 3 vpshock norm 5.250225E-05 +/- 0.138422E-04 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Chi-Squared = 239.4636 using 86 PHA bins. Reduced chi-squared = 3.031185 for 79 degrees of freedom Null hypothesis probability = 4.630E-18 XSPEC>addcomp 3 ga Input parameter value, delta, min, bot, top, and max values for ... 6.5 0.05 0 0 1E+06 1E+06 5:gaussian:LineE>0.714 0.1 0.05 0 0 10 20 6:gaussian:Sigma>0.0 1 0.01 0 0 1E+24 1E+24 7:gaussian:norm>= par4 * 4 Equating parameter :norm to parameter vpshock:norm * 4 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Model: phabs<1>( gaussian<2> + gaussian<3> + vpshock<4> ) Model Fit Model Component Parameter Unit Value par par comp 1 1 1 phabs nH 10^22 7.499971E-02 +/- 0.267359E-01 2 17 2 gaussian LineE keV 0.730000 +/- 0.00000 3 18 2 gaussian Sigma keV 0.00000 +/- 0.00000 4 19 2 gaussian norm 1.000000E-05 +/- 0.00000 5 20 3 gaussian LineE keV 0.714000 +/- 0.00000 6 21 3 gaussian Sigma keV 0.00000 +/- 0.00000 7 19 3 gaussian norm 4.000000E-05 = par 4 * +4.00000 8 2 4 vpshock kT keV 0.983738 +/- 0.132246 9 3 4 vpshock H 1.00000 frozen 10 4 4 vpshock He 1.00000 frozen 11 5 4 vpshock C 4.30000 frozen 12 6 4 vpshock N 2.60000 frozen 13 7 4 vpshock O 4.40000 frozen 14 16 4 vpshock Ne 1.50000 frozen 15 8 4 vpshock Mg 15.0000 frozen 16 9 4 vpshock Si 49.9000 frozen 17 9 4 vpshock S 49.9000 = par 16 18 7 4 vpshock Ar 8.80000 = par 13 * +2.00000 19 7 4 vpshock Ca 13.2000 = par 13 * +3.00000 20 10 4 vpshock Fe 1.00000 frozen 21 11 4 vpshock Ni 1.00000 frozen 22 12 4 vpshock Tau_l s/cm^3 0.00000 frozen 23 13 4 vpshock Tau_u s/cm^3 1.132493E+10 +/- 0.147843E+10 24 14 4 vpshock Redshift 0.00000 frozen 25 15 4 vpshock norm 5.250225E-05 +/- 0.138422E-04 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Chi-Squared = 218.2895 using 86 PHA bins. Reduced chi-squared = 2.834928 for 77 degrees of freedom Null hypothesis probability = 1.891E-15 XSPEC>addcomp 4 ga Input parameter value, delta, min, bot, top, and max values for ... 6.5 0.05 0 0 1E+06 1E+06 8:gaussian:LineE>0.723 0.1 0.05 0 0 10 20 9:gaussian:Sigma>0.0 1 0.01 0 0 1E+24 1E+24 10:gaussian:norm>= par4 * 2 Equating parameter :norm to parameter gaussian:norm * 2 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Model: phabs<1>( gaussian<2> + gaussian<3> + gaussian<4> + vpshock<5> ) Model Fit Model Component Parameter Unit Value par par comp 1 1 1 phabs nH 10^22 7.499971E-02 +/- 0.267359E-01 2 17 2 gaussian LineE keV 0.730000 +/- 0.00000 3 18 2 gaussian Sigma keV 0.00000 +/- 0.00000 4 19 2 gaussian norm 1.000000E-05 +/- 0.00000 5 20 3 gaussian LineE keV 0.714000 +/- 0.00000 6 21 3 gaussian Sigma keV 0.00000 +/- 0.00000 7 19 3 gaussian norm 4.000000E-05 = par 4 * +4.00000 8 22 4 gaussian LineE keV 0.723000 +/- 0.00000 9 23 4 gaussian Sigma keV 0.00000 +/- 0.00000 10 19 4 gaussian norm 2.000000E-05 = par 4 * +2.00000 11 2 5 vpshock kT keV 0.983738 +/- 0.132246 12 3 5 vpshock H 1.00000 frozen 13 4 5 vpshock He 1.00000 frozen 14 5 5 vpshock C 4.30000 frozen 15 6 5 vpshock N 2.60000 frozen 16 7 5 vpshock O 4.40000 frozen 17 16 5 vpshock Ne 1.50000 frozen 18 8 5 vpshock Mg 15.0000 frozen 19 9 5 vpshock Si 49.9000 frozen 20 9 5 vpshock S 49.9000 = par 19 21 7 5 vpshock Ar 8.80000 = par 16 * +2.00000 22 7 5 vpshock Ca 13.2000 = par 16 * +3.00000 23 10 5 vpshock Fe 1.00000 frozen 24 11 5 vpshock Ni 1.00000 frozen 25 12 5 vpshock Tau_l s/cm^3 0.00000 frozen 26 13 5 vpshock Tau_u s/cm^3 1.132493E+10 +/- 0.147843E+10 27 14 5 vpshock Redshift 0.00000 frozen 28 15 5 vpshock norm 5.250225E-05 +/- 0.138422E-04 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Chi-Squared = 218.1070 using 86 PHA bins. Reduced chi-squared = 2.908093 for 75 degrees of freedom Null hypothesis probability = 6.815E-16 XSPEC>freeze 2-3,5-6,8-9 Number of variable fit parameters = 5
При этом единственным дополнительным параметром, который мы ввели, является нормировка гауссиана с энергией 0.73 кэВ. Проверим его на статистическую значимость:
XSPEC>err 4 Parameter Confidence Range ( 2.706) 4 6.732174E-06 1.129377E-05 ( -2.298774E-06, 2.262822E-06)
Излишек на 1.4-1.6 кэВ объясняется невычтенной узкой инструментальной линией Si Kα (1.49 кэВ). Добавим ее в модель:
XSPEC>addcomp 5 ga Input parameter value, delta, min, bot, top, and max values for ... 6.5 0.05 0 0 1E+06 1E+06 11:gaussian:LineE>1.49 0.1 0.05 0 0 10 20 12:gaussian:Sigma>0.0 1 0.01 0 0 1E+24 1E+24 13:gaussian:norm>1e-05 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Model: phabs<1>( gaussian<2> + gaussian<3> + gaussian<4> + gaussian<5> + vpshock<6> ) Model Fit Model Component Parameter Unit Value par par comp 1 1 1 phabs nH 10^22 7.500000E-02 +/- 0.334087E-01 2 17 2 gaussian LineE keV 0.730000 frozen 3 18 2 gaussian Sigma keV 0.00000 frozen 4 19 2 gaussian norm 9.106036E-06 +/- 0.300654E-05 5 20 3 gaussian LineE keV 0.714000 frozen 6 21 3 gaussian Sigma keV 0.00000 frozen 7 19 3 gaussian norm 3.642414E-05 = par 4 * +4.00000 8 22 4 gaussian LineE keV 0.723000 frozen 9 23 4 gaussian Sigma keV 0.00000 frozen 10 19 4 gaussian norm 1.821207E-05 = par 4 * +2.00000 11 24 5 gaussian LineE keV 1.49000 +/- 0.00000 12 25 5 gaussian Sigma keV 0.00000 +/- 0.00000 13 26 5 gaussian norm 1.000000E-05 +/- 0.00000 14 2 6 vpshock kT keV 0.935248 +/- 0.106334 15 3 6 vpshock H 1.00000 frozen 16 4 6 vpshock He 1.00000 frozen 17 5 6 vpshock C 4.30000 frozen 18 6 6 vpshock N 2.60000 frozen 19 7 6 vpshock O 4.40000 frozen 20 16 6 vpshock Ne 1.50000 frozen 21 8 6 vpshock Mg 15.0000 frozen 22 9 6 vpshock Si 49.9000 frozen 23 9 6 vpshock S 49.9000 = par 22 24 7 6 vpshock Ar 8.80000 = par 19 * +2.00000 25 7 6 vpshock Ca 13.2000 = par 19 * +3.00000 26 10 6 vpshock Fe 1.00000 frozen 27 11 6 vpshock Ni 1.00000 frozen 28 12 6 vpshock Tau_l s/cm^3 0.00000 frozen 29 13 6 vpshock Tau_u s/cm^3 1.158069E+10 +/- -1.00000 30 14 6 vpshock Redshift 0.00000 frozen 31 15 6 vpshock norm 5.187108E-05 +/- 0.177663E-04 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Chi-Squared = 167.1128 using 82 PHA bins. Reduced chi-squared = 2.258281 for 74 degrees of freedom Null hypothesis probability = 3.724E-09 XSPEC>freeze 11-12 Number of variable fit parameters = 6 XSPEC>fit Chi-Squared Lvl Fit param # 1 2 13 15 19 26 160.564 -1 7.4718E-02 0.8971 1.1708E+10 5.3510E-05 9.7623E-06 3.1238E-06 158.485 -1 7.4971E-02 0.8627 1.2127E+10 5.4678E-05 9.2726E-06 8.6493E-06 156.998 -1 7.4996E-02 0.8465 1.2207E+10 5.5355E-05 9.7185E-06 4.2181E-06 156.500 -1 7.5000E-02 0.8366 1.2329E+10 5.5780E-05 9.4058E-06 7.8205E-06 150.478 0 7.5000E-02 0.8309 1.2345E+10 5.4941E-05 9.3253E-06 6.2483E-06 149.976 0 7.5000E-02 0.8309 1.2436E+10 5.4644E-05 9.4187E-06 6.2489E-06 149.768 0 7.5000E-02 0.8315 1.2503E+10 5.4482E-05 9.4391E-06 6.2390E-06 149.663 0 7.5000E-02 0.8320 1.2550E+10 5.4382E-05 9.4428E-06 6.2301E-06 149.603 0 7.5000E-02 0.8325 1.2581E+10 5.4315E-05 9.4425E-06 6.2233E-06 149.589 0 7.5000E-02 0.8319 1.2613E+10 5.4294E-05 9.4512E-06 6.2228E-06 149.566 0 7.5000E-02 0.8323 1.2624E+10 5.4265E-05 9.4431E-06 6.2165E-06 149.542 0 7.5000E-02 0.8320 1.2624E+10 5.4229E-05 9.4275E-06 6.2164E-06 149.526 0 7.5000E-02 0.8319 1.2625E+10 5.4201E-05 9.4175E-06 6.2168E-06 149.516 0 7.5000E-02 0.8318 1.2626E+10 5.4177E-05 9.4114E-06 6.2173E-06 149.508 0 7.5000E-02 0.8318 1.2628E+10 5.4158E-05 9.4080E-06 6.2179E-06 --------------------------------------------------------------------------- Variances and Principal axes : 1 2 13 15 19 26 2.51E-12 | 0.00 0.00 0.00 0.96 0.27 0.09 4.35E-12 | 0.00 0.00 0.00 -0.07 -0.06 1.00 7.92E-12 | 0.00 0.00 0.00 0.28 -0.96 -0.04 1.51E-04 | 0.94 0.34 0.00 0.00 0.00 0.00 1.14E-02 | -0.34 0.94 0.00 0.00 0.00 0.00 1.09E+05 | 0.00 0.00 -1.00 0.00 0.00 0.00 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Model: phabs<1>( gaussian<2> + gaussian<3> + gaussian<4> + gaussian<5> + vpshock<6> ) Model Fit Model Component Parameter Unit Value par par comp 1 1 1 phabs nH 10^22 7.500000E-02 +/- 0.378901E-01 2 17 2 gaussian LineE keV 0.730000 frozen 3 18 2 gaussian Sigma keV 0.00000 frozen 4 19 2 gaussian norm 9.408034E-06 +/- 0.310952E-05 5 20 3 gaussian LineE keV 0.714000 frozen 6 21 3 gaussian Sigma keV 0.00000 frozen 7 19 3 gaussian norm 3.763213E-05 = par 4 * +4.00000 8 22 4 gaussian LineE keV 0.723000 frozen 9 23 4 gaussian Sigma keV 0.00000 frozen 10 19 4 gaussian norm 1.881607E-05 = par 4 * +2.00000 11 24 5 gaussian LineE keV 1.49000 frozen 12 25 5 gaussian Sigma keV 0.00000 frozen 13 26 5 gaussian norm 6.217873E-06 +/- 0.213054E-05 14 2 6 vpshock kT keV 0.831823 +/- 0.100795 15 3 6 vpshock H 1.00000 frozen 16 4 6 vpshock He 1.00000 frozen 17 5 6 vpshock C 4.30000 frozen 18 6 6 vpshock N 2.60000 frozen 19 7 6 vpshock O 4.40000 frozen 20 16 6 vpshock Ne 1.50000 frozen 21 8 6 vpshock Mg 15.0000 frozen 22 9 6 vpshock Si 49.9000 frozen 23 9 6 vpshock S 49.9000 = par 22 24 7 6 vpshock Ar 8.80000 = par 19 * +2.00000 25 7 6 vpshock Ca 13.2000 = par 19 * +3.00000 26 10 6 vpshock Fe 1.00000 frozen 27 11 6 vpshock Ni 1.00000 frozen 28 12 6 vpshock Tau_l s/cm^3 0.00000 frozen 29 13 6 vpshock Tau_u s/cm^3 1.262839E+10 +/- 329.576 30 14 6 vpshock Redshift 0.00000 frozen 31 15 6 vpshock norm 5.415771E-05 +/- 0.217733E-04 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Chi-Squared = 149.5081 using 82 PHA bins. Reduced chi-squared = 1.967212 for 76 degrees of freedom Null hypothesis probability = 1.016E-06
С другой стороны, повышение содержания элементов понизило температуру плазмы - ведь для описания тепловых линий с бОльшим содержанием тяжелых элементов требуется более низкая температура. В результате модель перестала хорошо описывать спектр на больших энергиях. Поэтому в дальнейшем анализе мы ограничимся исследованием области 0.4-3.0 кэВ, в которой доминирует тепловое излучение:
XSPEC>ig 3.0-** ignoring channels 61 - 87 in dataset 1 ignoring channels 63 - 222 in dataset 2 Chi-Squared = 149.4546 using 82 PHA bins. Reduced chi-squared = 1.966508 for 76 degrees of freedom Null hypothesis probability = 1.030E-06 XSPEC>fit Chi-Squared Lvl Fit param # 1 2 13 15 19 26 Fit parameter 1 has pegged 149.454 3 7.5000E-02 0.8415 1.2584E+10 5.3711E-05 9.4070E-06 6.1832E-06 149.454 3 7.5000E-02 0.8414 1.2584E+10 5.3711E-05 9.4070E-06 6.1832E-06 --------------------------------------------------------------------------- Variances and Principal axes : 1 2 13 15 19 26 2.48E-12 | 0.00 0.00 0.00 0.96 0.27 0.08 4.35E-12 | 0.00 0.00 0.00 -0.07 -0.06 1.00 7.91E-12 | 0.00 0.00 0.00 0.28 -0.96 -0.04 1.44E-04 | 0.92 0.40 0.00 0.00 0.00 0.00 5.53E-03 | -0.40 0.92 0.00 0.00 0.00 0.00 1.39E+06 | 0.00 0.00 -1.00 0.00 0.00 0.00 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Model: phabs<1>( gaussian<2> + gaussian<3> + gaussian<4> + gaussian<5> + vpshock<6> ) Model Fit Model Component Parameter Unit Value par par comp 1 1 1 phabs nH 10^22 7.500000E-02 +/- 0.316499E-01 2 17 2 gaussian LineE keV 0.730000 frozen 3 18 2 gaussian Sigma keV 0.00000 frozen 4 19 2 gaussian norm 9.407020E-06 +/- 0.298745E-05 5 20 3 gaussian LineE keV 0.714000 frozen 6 21 3 gaussian Sigma keV 0.00000 frozen 7 19 3 gaussian norm 3.762808E-05 = par 4 * +4.00000 8 22 4 gaussian LineE keV 0.723000 frozen 9 23 4 gaussian Sigma keV 0.00000 frozen 10 19 4 gaussian norm 1.881404E-05 = par 4 * +2.00000 11 24 5 gaussian LineE keV 1.49000 frozen 12 25 5 gaussian Sigma keV 0.00000 frozen 13 26 5 gaussian norm 6.183179E-06 +/- 0.211465E-05 14 2 6 vpshock kT keV 0.841450 +/- 0.683292E-01 15 3 6 vpshock H 1.00000 frozen 16 4 6 vpshock He 1.00000 frozen 17 5 6 vpshock C 4.30000 frozen 18 6 6 vpshock N 2.60000 frozen 19 7 6 vpshock O 4.40000 frozen 20 16 6 vpshock Ne 1.50000 frozen 21 8 6 vpshock Mg 15.0000 frozen 22 9 6 vpshock Si 49.9000 frozen 23 9 6 vpshock S 49.9000 = par 22 24 7 6 vpshock Ar 8.80000 = par 19 * +2.00000 25 7 6 vpshock Ca 13.2000 = par 19 * +3.00000 26 10 6 vpshock Fe 1.00000 frozen 27 11 6 vpshock Ni 1.00000 frozen 28 12 6 vpshock Tau_l s/cm^3 0.00000 frozen 29 13 6 vpshock Tau_u s/cm^3 1.258448E+10 +/- 1179.21 30 14 6 vpshock Redshift 0.00000 frozen 31 15 6 vpshock norm 5.371057E-05 +/- 0.166385E-04 --------------------------------------------------------------------------- --------------------------------------------------------------------------- Chi-Squared = 149.4543 using 82 PHA bins. Reduced chi-squared = 1.966504 for 76 degrees of freedom Null hypothesis probability = 1.030E-06 XSPEC>err 1,4,13,14,29,31 Parameter Confidence Range ( 2.706) Apparent non-monotonicity in Chi-Squared space detected Current bracket values : 6.717281E-02 6.253970E-02 and delta Chi-Squared : 2.66106 4.31946 but latest trial 6.701802E-02 gives Chi-Squared = 2.63644 Suggest that you check this result using the steppar command Chi-Squared when model parameter 1= 7.500000E-02 is 149.4328, which is < previous minimum 149.4543 (critical delta = 0.0100) Chi-Squared when model parameter 1= 7.500000E-02 is 149.4160, which is < previous minimum 149.4328 (critical delta = 0.0100) Apparent non-monotonicity in Chi-Squared space detected Current bracket values : 6.723657E-02 6.267269E-02 and delta Chi-Squared : 2.67221 4.30505 but latest trial 6.495463E-02 gives Chi-Squared = 4.36517 Suggest that you check this result using the steppar command Parameter pegged at hard limit 7.500000E-02 with delta ftstat= 0.00000 1 6.495463E-02 7.500000E-02 ( -1.004537E-02, 0.00000 ) 4 7.154419E-06 1.168864E-05 ( -2.262671E-06, 2.271549E-06) 13 4.473682E-06 7.924309E-06 ( -1.739065E-06, 1.711563E-06) 14 0.730027 0.943957 ( -0.103749 , 0.110181 ) Number of trials exceeded before convergence. Current trial values : 1.269185E+10 1.065784E+10 and delta Chi-Squared : 0.00000 3.91721 Minimization may have run into a problem, check your result Minimization may have run into a problem, check your result Number of trials exceeded before convergence. Current trial values : 1.269185E+10 1.675987E+10 and delta Chi-Squared : 0.00000 7.64384 Apparent non-monotonicity in Chi-Squared space detected Current bracket values : 1.512598E+10 1.512601E+10 and delta Chi-Squared : 2.67252 2.75343 but latest trial 1.512599E+10 gives Chi-Squared = 2.67247 Suggest that you check this result using the steppar command 29 1.096681E+10 1.512599E+10 ( -1.725041E+09, 2.434145E+09)
Стоит отметить, что величина χ2 заметно больше 1.0. Данное обстоятельство вызвано, прежде всего, несовершенством тепловой модели, описывающей систему. Однако обратное, в общем случае, не верно! Дело в том, что, используя статистику χ2, подразумевается гауссово распределение ошибок. Это верно в случае доминирования сигнала от источника. Однако на высоких энергиях, в особенности для источников с тепловым спектром, доминирует фон, и ошибка задается не числом фотонов от источника Nsource, а числом фотонов от фона
. Поэтому ошибка на больших энергиях большая, и χ2 занижается. Поэтому при доминировании фона нет ничего страшного в том, что χ2 будет меньше 1, и, вообще говоря, ничего почетного для модели от того, что χ2 будет в районе 1.
Другая причина несоответствия модели состоит в том, что не был правильно вычтен космический фон. Данное обстоятельство является важным, т.к, при наличии количества N фотонов в бине (
для низких энергий) фон должен быть вычтен с относительной точностью не менее 
, чтобы не испортить χ2. Поскольку данное исследование здесь не было проведено, приведенные результаты имеют исключительно иллюстративный характер.
Полезные дополнительные процедуры
Полезными процедурами, кроме уже перечисленных, являются:
- flux - подсчитывает поток излучения в заданном диапазоне. Может посчитать и его ошибку.
- fakeit Если данных нет, приходится их симулировать. Полезно для написания пропозалов или (что, в принципе, то же самое) для теоретиков, которые хотят предсказать некий эффект.
- steppar Если процедура error занимает много времени, единственный выход - оценить ошибку таким образом.
- goodness Оценивает качество фита с помощью Монте Карло симуляций.
- identify Идентифицирует линии излучения в некоторых моделях.
- ftest Сравнивая две модели на одном и том же наборе данных, указывает, какая из них лучше (это не всегда тривиальная процедура).
- iplot Изменяет изображение спектра (перемасштабирование, подписывание осей и прочее).
- hardcopy Печатает изображение спектра (в .ps формате)
- systematic Задает систематическую ошибку в спектре. Часто бывает необходимой.
- time Выдает CPU время, потраченное на рассчеты.
Переменные xspec11
У xspec11 есть несколько важных переменных. Ниже я их перечислю:
- chatter. Самый важный параметр, от него зависит полнота выдаваемой на экран информации. По умолчанию равен 10.
- cosmo. Космологические параметры - H0, q0, ΩΛ,0. Необходима для корректного вычисления потока от далеких объектов.
- xsect. Сечения реакций (необходимы для правильного вычисления поглощения и теплового излучения). По умолчанию это bcmc.
- statistic. Используемая статистика. Для гауссовых распределений подходит chi, используемая по умолчанию.
- weight. Разные бины могут входить с разным весом. По умолчанию вес каждого бина одинаков (standard).
Сохранение сессии
Сохранить последние результаты фитирования:
XSPEC> save all work.xcmНедостатком является то, что xspec11 не записывает туда результат команд setplot и query. Приходится редактировать work.xcm. Еще один недостаток - не записываются ответы команд flux, error и т.д. Поэтому важно сохранять логи, например, запуская xspec11 так:
$ xspec11 |& tee xspec.logВыход из xspec11
Стандартный выход:
XSPEC> exit Do you really want to exit? (y)y XSPEC: quit
При этом результаты последнего фита автоматически сохраняются в файле xautosav.xcm.
Запуск сохраненной сессии
XSPEC>@work.xcmЗапуск xspec11 в рамках скрипта
$ xspec11 < work.xcm |& tee xspec.log
При всей удобности процедуры отмечу, что в конце файла work.xcm, запускаемого таким образом, обязано стоять следующее:
exit
y
В противном случае xspec11 ругнется и напишет огромный лог (который в худшем случае заполнит все доступное пространство на диске). Впрочем, это не единственное место, где он может ругнуться (нет на месте требуемого файла, например), поэтому применять его надо с осторожностью! И на всякий случай, если несчастье все же произошло, "держите под рукой" команду
$ killall xspec11На этой оптимистичной ноте я заканчиваю свое повествование :-)
