Xspec11 for beginners

From WikiVirgo

Jump to: navigation, search

by Dima Iakubovskyi

После получения спектров в рентгеновском диапазоне возникает желание вытянуть оттуда физические параметры, такие как поток излучения (в таком-то диапазоне), параметры моделей и т.д. Наиболее распространенным приложением, которое позволяет это сделать, является пакет xspec11, входящий в пакет программ по астрофизике высоких энергий HEASOFT.

В этой статье я попытаюсь последовательно изложить основные возможности пакета xspec11, не претендуя при этом ни на полноту изложения, ни на ее глубину.

Contents

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\cdot 10^{20}cm − 2, что является оценкой сверху для галактического объекта. С другой стороны, несмотря на то, что SN1006 - объект галактический, он находится довольно высоко над галактическим диском (на высоте около 550 пк), поэтому величина столбовой плотности газа перед SN1006 не должна сильно отличаться от галактической. Мы ограничим ее в области 6.0-7.5\cdot 10^{20}cm − 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, а числом фотонов от фонаN_{back} \gg N_{source}. Поэтому ошибка на больших энергиях большая, и χ2 занижается. Поэтому при доминировании фона нет ничего страшного в том, что χ2 будет меньше 1, и, вообще говоря, ничего почетного для модели от того, что χ2 будет в районе 1.

Другая причина несоответствия модели состоит в том, что не был правильно вычтен космический фон. Данное обстоятельство является важным, т.к, при наличии количества N фотонов в бине (N \gg 1 для низких энергий) фон должен быть вычтен с относительной точностью не менее 1/\sqrt{N}, чтобы не испортить χ2. Поскольку данное исследование здесь не было проведено, приведенные результаты имеют исключительно иллюстративный характер.

Полезные дополнительные процедуры

Полезными процедурами, кроме уже перечисленных, являются:

  • flux - подсчитывает поток излучения в заданном диапазоне. Может посчитать и его ошибку.
  • fakeit Если данных нет, приходится их симулировать. Полезно для написания пропозалов или (что, в принципе, то же самое) для теоретиков, которые хотят предсказать некий эффект.
  • steppar Если процедура error занимает много времени, единственный выход - оценить ошибку таким образом.
  • goodness Оценивает качество фита с помощью Монте Карло симуляций.
  • identify Идентифицирует линии излучения в некоторых моделях.
  • ftest Сравнивая две модели на одном и том же наборе данных, указывает, какая из них лучше (это не всегда тривиальная процедура).
  • improve По отзывам, штука полезная, однако требует MINUIT. В данный момент на ВИРГО не установлена.
  • iplot Изменяет изображение спектра (перемасштабирование, подписывание осей и прочее).
  • hardcopy Печатает изображение спектра (в .ps формате)
  • systematic Задает систематическую ошибку в спектре. Часто бывает необходимой.
  • time Выдает CPU время, потраченное на рассчеты.

Переменные xspec11

У xspec11 есть несколько важных переменных. Ниже я их перечислю:

  • chatter. Самый важный параметр, от него зависит полнота выдаваемой на экран информации. По умолчанию равен 10.
  • abund. Модель абсолютного содержания элементов в Солнце. По умолчанию это angr.
  • cosmo. Космологические параметры - H0, q0, ΩΛ,0. Необходима для корректного вычисления потока от далеких объектов.
  • xsect. Сечения реакций (необходимы для правильного вычисления поглощения и теплового излучения). По умолчанию это bcmc.
  • method. Метод минимизации. По умолчанию это leven.
  • 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

На этой оптимистичной ноте я заканчиваю свое повествование :-)

Software

Personal tools