ksl130430_x-sections - PHOTOIONIZATION X-SECTIONS

[jhmatthews]

Current Status: Understood
Priority: Highest -- need for Nick's paper
Version discovered in: - 75, but it will be in any spectrum
Brief statement of the problem: For Topbase xsections, a crossection of 0 is returned from sigma_phot_topbase above the maximum of the photon frequencies for which data is provided
How to reproduce: Any model in which one expects significant X-ray absorption, including the fiducial model for Nick's paper
Discussion: TopBase does not provide an extrapolation formula for x-sections outside of the energy range for which data is provided. When we calculate photoionization x-sections from topbase we perform a linear interpolation between the data provided, but we did not add our on extrapolation for frequencies beyond that given in topbase.

There is a problem because TopBase X-sections, which are now the preferred type of x-section in Python do not extent to very high energies

TopBase simply does not provide data to very high frequencies, probably because it was set up for stellar atmospheres. (Note that we use TopBase because it has partial ionization x-sections and because the other information it provides is convenient for macro-atoms). There are several potential solutions:

• Add an extrapolation formula. Probably the best way to doe this would be to determine the power law index of the photoabsorption at the maxium energy (nu**-alpha) provided and use this to extrapolate
• Drop back to Verner x-sections for Nick's paper. This may be a little dangerous, but ksl believes (but did not check) that the Verner X-sections do not have this problem.
• Long term, move to Chianti X-sections which almost certainly do not have this problem, but we need to do this in such a way that macro atoms can still be calculated.

130506 - ksl looked for photoionization x-sections in chianti and they have "disappeared" in the latest releases. Stuart confirmed that he really had gotten them from chianti. A query to the chianti folks resulted in a response from Ken Dere: (Ken Dere [email protected]) who said"

that part of things never went very far but I had obtained the data from:

http://amdpp.phys.strath.ac.uk/DATA/PI/

This does have a lot of photo-ionization data.

JM130625: I have looked into the strathclyde data. The documentation is somewhat unclear, (for example, no units on energy that I can find), which makes drawing a conclusion at the moment different. I'm also not sure I understand their splitting of data sets into IC and LS coupling. I also noticed this sentence

''The inner-shell data initial states are a subset of those computed for the Updated Opacities (which considered excited initial complexes as well) but cover more elements (all up to Zn, plus Kr, Mo, Xe) and in more detail (using non-orthogonal relaxed orbitals).''

Personally, I think we should just extrapolate topbase photo-XSs outside of python, with a script which will work on topbase files in general.

Please also consult the X-ray attenuation page

[jhmatthews]

@kslong suggested just doing an extrapolation by extending the gradient in loglog space near the maximum- this works fine for most x-sections, but there are a few problematic ones.

Here's an example XS of topbase

Figure 2

which is the XS for

PhotTopS 8 1 501 8 0.17605

i.e. an excited state of OVIII

I can just extrapolate ones like this to a very low value but I'm not sure if that is the right thing to do.

[jhmatthews]

I've got a version of the extrapolated atomic data with topbase, and I'm doing a few runs with it to see what the X-ray spectra look like. You can find it, along with the masterfile on

~/Dropbox/Python/atomic_data/

or also by just pulling the changes from the data branch, or cloning it if you haven't done so. Alternatively, go to https://github.com/agnwinds/python/tree/data, click on Zip and then stick that in your python directory- it has an updated Setup_Py_Dir script in there.

Now in contact with Adam Foster and Randall Smith regarding obtaining new XSs.

[jhmatthews]

Here is a comparison for the fiducial model, from the 85 degree sightline, between standard_extrap and standard73. I've cleaned up these comments a little as it was difficult to follow the thread.

Figure 6

[jhmatthews]

I've started looking at the AtomDB data sent to me by Randall Smith and Adam Foster. If anyone else fancies having a lookg there is a tarball of a load of fits files here:
http://hea-www.cfa.harvard.edu/~afoster/xfer/APED_XSTAR.tar.bz2

I'm still working on understanding the data, but thought I would upload as I went along so anyone who is interested can take a look. Decided to look at Oxygen exclusively for the moment, as its produced some interesting issues. So here is just a plot of all the OV Xsections, for AtomDB data:

Figure 7

And the OV Xsections in my topbase extrapolate data set, standard_extrap:
Figure 8

Like I said, only just started on this so not sure entirely what to glean from this yet. Initial impressions are that by using topbase we are missing structure around the keV area, but the more striking observation of the AtomDB data is that it runs out at just above 2KeV in this case...which is actually worse than topbase.

I also took a look at the O8 Xsections.

Figure 9

Which tend to imply that the weird looking Xsections seen a few plots above in figure 2 are definitely a bit dodgy.

Will write more soon!

[jhmatthews]

A quick point on this- myself and Nick were just discussing, and he pointed out that the XSTAR/AtomDB data columns I've been sent have a 'Reference' associated with them- e.g. for oxygen 5 these are XSTAR Type 49 and XSTAR Type 53. For oxygen 8 these are XSTAR Type 53.

And, lo and behold if you go to the XSTAR website and check what these types are...

http://heasarc.gsfc.nasa.gov/docs/software/xstar/docs/html/node164.html

49) opacity project pi x-sections
Photoionization cross section from TOPbase averaged over resonances. Photon energies are in Ry with respect to the subshell ionization threshold and cross sections are in Mb. Just like 53, except for energy scale.. Every line contains: reals (2_np; x(i),y(i),i=1,np) ; integers (8; N, L, 2_J, Z, lev+, nion+, nlev, nion) ; characters (0) ; (ion nlev correspond to the initial state and ion+ lev+ correspond to the state to which that ionizes.)

53) opacity project pi x-sections
Photoionization cross section from TOPbase averaged over resonances. Photon energies are in Ry with respect to the first ionization threshold and cross sections are in Mb. Every line contains: reals (2_np; x(i),y(i),i=1,np) ; integers (8; N, L, 2_J, Z, lev+, nion+, nlev, nion) ; characters (0) ; (ion nlev correspond to the initial state and ion+ lev+ correspond to the state to which that ionizes.)

So they are TOPbase Xsections!

[jhmatthews]

I've created a page for this here

[Higginbottom]

Just to add an extra comment here: I have noted when profiling python in instruments that during the ionization cycles, the main time sink is doing integrations over VERNER cross sections. This is done when computing recombination rates for ions that do not have badnell rates (via the milne relaion) and also computing the denominator in our pairwise photoionization rate correction factor.
Digging deeper, it is the execution of the pow command in the verner equation that takes all the time.
We should consider having a consistent photoionization rate format for all ions, tabulating verner data if necessary. This could have a significant effect on the execution time.
FWIW, I have been able to reduce the execution time by a noticable amount on my mac by computing the pow commands via an exponential form - i.e. x^y = exp(y*ln(x))

[kslong]

Hmm, I don't doubt your result that a different formulation of pow produces a different result, alhough I have to say I find it quite surprising. I wonder if this is information known on the web

Anyway, I agree that we ought to look into the question of how to represent the verner data. Probably we ought to set up a separate google page to lay out a strategy for where we want to go with atomic data. (James has a page set up, but it is more about the problems) There is a lot of material that is not well documented about what is done. When I wrote atomic data, I did not appreciate that we might want to have a generic way of describing the photoionization x-sections, so I used the word Topbase to describe a specific kind of data. If we go the suggested route, Topbase would refer to way of digitizing data as opposed to a source. This gives plenty of opportunity for confusion that we should work out how to deal with.

[Higginbottom]

The problems with the pow command is certainly mentioned on the web - and the exp/ln formulation was taken from a discussion there. Other ideas are approximate version or at the most extreme - machine code - but if we want python to remain portable.....

Just to throw an idea out - we could consider having a binary file that we pre-prepare which contains all the atomic data using a seperate code. We could produce what we think is the best data - extrapolated - interpolated etc and then one would just read in one file, and turn on what one wants....

[kslong]

I am not in favor of binary files for atomic data, and do not see the need. Binary files can easily be different on different operating systems.

[Higginbottom]

I have incorporated and tested a little subroutine that populates a structure of type topbase_phot with sampled verner cross sections.
Using this structure means that existing code can be used for integrations and interpolations, however as Knox points out, this is slightly confusing. Of course, this interpolation could be done offline, and the data read in as an interpolated file. We should discuss it - the speed up is very significant. A simple test working out ionization fractions for a benchmark model dropped from 171s to 33s.
I have tested the implementation by comparison with cloudy and the plot is just the same as with non-tabulated verner data.

I've not committed the change to the dev branch yet - would you like me to?

[kslong]

I think committing to the dev branch is fine, but this may be only a short term fix. I currently favor creating some kind of helper routine to put the VFKY data into the format we would
like to read in with getatomic_data. At any rate we need to continue this discussion.

Knox

[Higginbottom]

I agree totally that this is a temporary fix - and hopefully James's grand unifiied PI cross section scheme will render it obsolete. However, it halves the execution time for the proga models, so I think it is a handy stopgap!!

[jhmatthews]

You can now take a look at every cross section plotted with this link- Which should be useful for reference

https://www.dropbox.com/sh/mp9q7w9tzxvsdx2/FNNoSpNPdB

I will host these on my website at some point but at the moment the servers are down. There's a readme which explains the format and colours.

Edit: there's now a form here which would do that

The reason I did this was because we want to start using the standard_extrap data set, and so the first order check I did was to check the gradients of the extrapolation, and for the script to warn me if a gradient is >-3. This happens in a lot of cases. Some of these are probably real (I don't understand the details of what determines the gradient) and some not. For example, if you look at some of the CNO cross sections there are examples (such as XS_8_3_141.png) where there is a slight upwards curve at the end of the topbase data which then causes the extrapolation to be with a gradient of close to -1 or even 0 in some cases on the log-log scale.

[sirocco-rt]

The suggested solution is as follows, as discussed on 27 Sept:

We are going to use topbase cross sections but with the following tests and adaptations:

• All these improvements will prioritise by abundance and whether we are interested in that element (e.g. start with O, N, C, H, and go from there)
• We will gradually correct the thresholds so they agree with AtomDB and/or Verner.
• The extrapolated data set will be used from now on, but we will improve and test it. JM now has a Topbase XS query form coded up here. Action for JM to look into the XSections with less negative gradients, assess if any are important/abundant, and if they are then look in more detail at related literature
• We will try and encorporate inner shell edges from AtomDB into topbase, especially for important states
• We need to document this well, and the data branch on github should be used for atomic data tracking

[jhmatthews]

I'm finally sorting this out. The procedure is pretty klugey, but here's what I'm doing

• identify odd XSections by eye
• Check none of them are the ground state
• If the gradient is above ~3 (I used 2.9 to give a slight buffer) and it is in the 'odd looking' list, then fudge it by setting gradient to -3
• if not, use the gradient at the last point in the data
• extrapolate up to 100Kev

Here's an example of a fixed problematic XSection:

compared to what it looked like previously:

This will now become the standard77 data set. I will test and then make a release. Ideally, anyone using the atomic data should be synced up with the data branch, so that way they can just pull in the changes. Once I've updated and tested, we can close this issue, remembering that our data is not perfect and the threshold energies are still slightly off.

We may also want to just get Python to print out the recombination coefficients with and without the extrapolation, just so we have peace of mind regarding integration of these XSections.

[jhmatthews]

I've compared fiducial models with standard73 and standard77 atomic data sets, and they look identical in the UV. I will double check that the X-ray edge disappears as above and then this issue is closed with the caveat that thresholds have not been adjusted.

The individual cross sections I've corrected are shown here - you can see pretty well the fudge made.

@Higginbottom - you should pull in the changes from the data branch and use standard77

[jhmatthews]

Stop Press- don't use this quite yet @Higginbottom . Some slightly odd things happening in the X-ray spectra. Remember, UV spectra look identical here, but I must be doing something wrong. Remember I have 3 atomic data sets:

standard73
standard_extrap: My original extrapolated set without any fudging of odd XSections
standard77: the most recent one with the fixing of odd looking XSections

standard73 v standard_extrap

standard77 v standard_extrap

it looks like I've messed something up as it would seem that I've managed to create more bf opacity in the region that didn't get affected by the extrapolation originally...

[jhmatthews]

Good news, this was a simple error in the script on my part, so the standard77 data set now looks like this (as it should):

Nick can now go ahead and use this data.

I consider this issue closed now, with the following caveats:

• thresholds have not been adjusted
• we should be aware that some XSections have been fudged with a gradient of -3, but these are never ground states
• we should be aware that there are some XSections which have been extrapolated but perhaps shouldn't be. As these are partial
• luckily, none of this affects the fiducial model at all except the postive outcome of getting rid of the unwanted x-ray edge

I will document what I've done properly and post somewhere appropriate, and also make the scripts available.

[kslong]

That's good news. Knox