Fitting XMM-Newton data with the APEC model

The Astrophysical Plasma Emission Code (APEC, Smith et al. 2001) is a state-of-the-art code to model the X-ray emission of optically thin astrophysical plasma in collisional equilibrium. APEC has been widely used in the literature to study the spectra of a wide variety of sources, such as galaxy clusters and groups, supernova remnants, or stellar coronae. The code self-consistently models the continuum (bremsstrahlung) emission as well as the pseudo-continuum and over 10,000 individual emission lines.

Astromodels includes a native Python implementation of APEC based on pyatomdb <https://atomdb.readthedocs.io/en/master/>__. As the model is only available when an installation of pyatomdb can be found, the user must make sure pyatomdb is installed and the required atomic data files have been downloaded.

Setting up pyatomdb

pyatomdb is easily installable using pip:

$] pip install pyatomdb  When running the code for the first time, the user will be requested to download the ATOMDB database locally $] python -c "import pyatomdb"


This command will prompt the user to choose a directory to store the atomic data (usually $HOME/atomdb). Once the data have been downloaded, the ATOMDB environment variable must be set to indicate the location of the atomic data $] export ATOMDB=$HOME/atomdb  OK, we’re all set! Loading the APEC model Once pyatomdb is properly installed, it is available as a “Function1D” object [1]:  %%capture from threeML import * modapec = APEC()  The intensity of the various lines in the model is set relative to the Solar abundance. Therefore, one must set the Solar abundance table to predict the line intensity, using the init_session method of the APEC class. By default, i.e. if the init_session method is ran with no argument, the code defaults to Anders & Grevesse (1989). [2]:  modapec.init_session(abund_table='AG89')  Will not thermally broaden lines  If, for instance, we wish to initialize the APEC model to use the Lodders & Palme (2009) table, we pass modapec.init_session(abund_table='Lodd09')  [3]:  modapec.display()  • description: The Astrophysical Plasma Emission Code (APEC, Smith et al. 2001) contributed by Dominique Eckert • formula:$n.a.$• parameters: • K: • value: 1.0 • desc: Normalization in units of 1e-14/(4*pi*(1+z)^2*dA*2)*EM • min_value: 1e-30 • max_value: 1000.0 • unit: • is_normalization: True • delta: 0.1 • free: True • kT: • value: 1.0 • desc: Plasma temperature • min_value: 0.08 • max_value: 64.0 • unit: • is_normalization: False • delta: 0.1 • free: True • abund: • value: 1.0 • desc: Metal abundance • min_value: 0.0 • max_value: 5.0 • unit: • is_normalization: False • delta: 0.01 • free: False • redshift: • value: 0.1 • desc: Source redshift • min_value: 0.0 • max_value: 10.0 • unit: • is_normalization: False • delta: 0.001 • free: False The parameters of the model are set in the following way [4]:  modapec.kT.value = 3.0 # 3 keV temperature modapec.K.value = 1e-3 # Normalization, proportional to emission measure modapec.redshift.value = 0. # Source redshift modapec.abund.value = 0.3 # The metal abundance of each element is set to 0.3 times the Solar abundance  Now let’s see how the model depends on the input temperature… [5]:  import numpy as np energies = np.logspace(-1., 1.5, 1000) # Set up the energy grid ktgrid = [0.2,0.5,1.0,2.0,3.0,5.0,7.0,9.0,12.0,15.0] # Temperature grid  [6]:  import matplotlib.pyplot as plt import matplotlib.colors as colors import matplotlib.cm as cmx %matplotlib inline plt.clf() fig=plt.figure(figsize=(13,10)) ax = fig.add_axes([0.12, 0.12, 0.85, 0.85]) for item in (ax.get_xticklabels() + ax.get_yticklabels()): item.set_fontsize(18) nspec = len(ktgrid) values = range(nspec) cm = plt.get_cmap('gist_rainbow') cNorm = colors.Normalize(vmin=0, vmax=values[-1]) scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=cm) ccc=[] for i in range(nspec): ccc.append(scalarMap.to_rgba(i)) for i in range(nspec): modapec.kT.value = ktgrid[i] plt.plot(energies,modapec(energies),color=ccc[i],label='kT=%g keV'%(ktgrid[i])) plt.xscale('log') plt.yscale('log') plt.xlabel('Energy [keV]',fontsize=28) plt.ylabel('Photon Flux',fontsize=28) plt.axis([0.1,15.,1e-8,1.0]) plt.title('Z=0.3$Z_\odot$and varying temperature',fontsize=28) plt.legend(fontsize=22, ncol=2)  [6]:  <matplotlib.legend.Legend at 0x7f4927f66940>  <Figure size 432x288 with 0 Axes>  Now let’s see how the model depends on metallicity for a temperature of 1 keV… [7]:  Zgrid = [0., 0.1, 0.3, 0.5, 1., 2.] # Metallicities wrt Solar modapec.kT.value = 1.0  [8]:  plt.clf() fig=plt.figure(figsize=(13,10)) ax = fig.add_axes([0.12, 0.12, 0.85, 0.85]) for item in (ax.get_xticklabels() + ax.get_yticklabels()): item.set_fontsize(18) nspec = len(Zgrid) values = range(nspec) cm = plt.get_cmap('gist_rainbow') cNorm = colors.Normalize(vmin=0, vmax=values[-1]) scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=cm) ccc=[] for i in range(nspec): ccc.append(scalarMap.to_rgba(i)) for i in range(nspec): modapec.abund.value = Zgrid[i] plt.plot(energies,modapec(energies),color=ccc[i],label='$Z/Z_\odot$=%g'%(Zgrid[i])) plt.xscale('log') plt.yscale('log') plt.xlabel('Energy [keV]',fontsize=28) plt.ylabel('Photon Flux',fontsize=28) plt.axis([0.1,15.,1e-7,0.1]) plt.title('kT=1 keV and varying metallicity',fontsize=28) plt.legend(fontsize=22)  [8]:  <matplotlib.legend.Legend at 0x7f494f7edee0>  <Figure size 432x288 with 0 Axes>  Comparison with XSPEC To test the implementation, let’s simulate a spectrum with XSPEC and fit it within 3ML using the APEC model. We simulate an XMM-Newton/EPIC-pn spectrum with a temperature of 5 keV, a metallicity of 0.3 Solar, a normalization of unity and a redshift of 0.1. The simulated model is absorbed by photo-electric absorption (PhAbs model) with a column density of $$10^{21}$$ cm$$^{-2}$$. We start by declaring the model, [9]:  phabs = PhAbs() phabs.NH.value = 0.1 # A value of 1 corresponds to 1e22 cm-2 phabs.NH.fix = True # NH is fixed phabs.init_xsect(abund_table='AG89') modapec.kT = 3.0 # Initial values modapec.K = 0.1 modapec.redshift = 0.1 modapec.abund = 0.3 mod_comb = phabs * modapec  Now we load the simulated spectrum in 3ML… [10]:  xmm_pha = get_path_of_data_file("datasets/xmm/pnS004-A2443_reg2.fak") xmm_rmf = get_path_of_data_file("datasets/xmm/pnS004-A2443_reg2.rmf") xmm_arf = get_path_of_data_file("datasets/xmm/pnS004-A2443_reg2.arf") ogip = OGIPLike("ogip", observation=xmm_pha, response=xmm_rmf, arf_file=xmm_arf) pts = PointSource('mysource',0,0,spectral_shape=mod_comb)  Let’s have a look at the loaded spectrum… [11]:  fig = ogip.view_count_spectrum() fig.set_size_inches(13,10) ax = fig.get_axes()[0] ax.set_xlim(left=0.3,right=14.) ax.set_xlabel('Energy [keV]',fontsize=28) ax.set_ylabel('Rate [counts s$^{-1}$keV$^{-1}\$]',fontsize=28)


Here we set up the likelihood and fit the data

[12]:

ogip.remove_rebinning()

ogip.set_active_measurements('0.5-10.')

ogip.rebin_on_source(20)

model = Model(pts)

jl = JointLikelihood(model,DataList(ogip))

result = jl.fit()

Best fit values:


result unit
parameter
mysource.spectrum.main.composite.K_2 (9.8742 +/- 0.0033) x 10^-1 1 / (cm2 keV s)
mysource.spectrum.main.composite.kT_2 5.019 +/- 0.005 keV

Correlation matrix:


 1 -0.3 -0.3 1

Values of -log(likelihood) at the minimum:


-log(likelihood)
ogip 9682.695868
total 9682.695868

Values of statistical measures:


statistical measures
AIC 19369.398043
BIC 19380.497261

The fitted values are within 1% of the input ones from XSPEC, showing that the native 3ML implementation is accurate. We can visualize the fit in the following way

[13]:

fig = display_spectrum_model_counts(jl,data_color='blue',model_color='red', min_rate=5e-4)

fig.set_size_inches(13,10)