search

sirocco-rt / sirocco

This is the repository for Sirocco, the radiative transfer code used to model winds in AGN and other systems
GNU General Public License v3.0
26 stars 24 forks source link

Adding charge exchange #447

Open Higginbottom opened 5 years ago

Higginbottom commented 5 years ago

Just pottering around making plots of abundances in python and cloudy for matching IP, T and SED. Its been a while since I did this, and I note that charge exchange is now the only thing I have to turn off in cloudy to get a really super match….

Everything on in both cloudy and python (H,He, C,O from top to bottom) Solid is python, dotted is cloudy.

image001

turn off charge exchange in cloudy…

image002

Seems like these rates should just be quite simple things to put in. The relevant paper is: http://adsabs.harvard.edu/abs/1996ApJS..106..205K

which contains tabulated rate coefficients. There may well be issues with heating and cooling rates (as with dielectric recombination) and the rates are not that accurate anyway. We would be closer to cloudy, but still not 'correct' - but nobody is.

Higginbottom commented 4 years ago

This is the paper to deal with the heating effects https://ui.adsabs.harvard.edu/abs/1999ApJ...516L.107K/abstract

Higginbottom commented 4 years ago

Taking a sneaky peak at cloudy - it looks like they use the data from table 1 in Kingdon and Ferland plus either the total from table 2 if there is more than one level listed, or just the only level if there is just one. Table 1 and 2 give the 'recombination' rates, which is the rate for A+ + H -> A + H+, i.e. ions hitting a neutral hydrogen, and ionising the hydrogen whilst recombining themselves. There is also a (small) table of 'ionization' rates which is (presumably) A + H+ -> A+ + H - i.e. ionised hydrogen hitting an ion and stealing an electron. In most cases, the resulting ion is excited, and so if we simply implement all this, we are likely to improve our ionization calculations, but leave yet another leak - like the dielectric recombination situation where we have excited ions where we don't take account of the relaxation.

Higginbottom commented 4 years ago

Reading the paper on heating effects, it is interesting that they say the heating effect is strongest in regions where H and H+ co-exist and also is more important where the SED is hard. So this may go some way to mitigating some of our problems with regions of partial H ionization giving odd temperatures. Maybe...

Higginbottom commented 4 years ago

OK, so I have now implemented this. I have used the rates from https://ui.adsabs.harvard.edu/abs/1996ApJS..106..205K/abstract combined with the energy exchange data from https://ui.adsabs.harvard.edu/abs/1999ApJ...516L.107K/abstract and a few extra corrected rates dug out of cloudy. The relevant file in cloudy is atmdat_char_tran.cpp

The energy transfer is almost always a heating effect - so I am only implementing this as a heating effect - any negative energy changes are just offset against the heating effect when it is calculated.

The routines are all in charge_exchange.c and the new data file is charge_exchange.dat.

I have not included charge transfer for Helium - this data is scattered across many papers, and is not nearly as good as for hydrogen. There are also various corrections in cloudy, that I have not laboriously copied over, it feels a little too much like copying!!

Anyway, here is a comparison for Oxygen with charge exchange on - all the graph below have crosses for cloudy and lines for python.

oxygen_plot

and even iron looks better

Fe_abund_comp

and here are the heating rates (for fixed T=10000K)

aaaa_heat

Looks good eh? I'm going to run regression tests, and comment everything then potentially pull the changes in later...

Higginbottom commented 4 years ago

I should note that this method does have one wrinkle.

The rate depends on both the number density of hydrogen and that of the ion in question. It is not possible to create a term in the rate matrix that depends on two densities, so the density of hydrogen has to be fixed in the same way as the electron density prior to trying to solve the matrix.. So in the case where hydrogen is not super dominant in terms of overall numbers of particles, one could envision problems,

Higginbottom commented 4 years ago

I had to decide where to compute and add in the heating effect. It made most sense to do this at the end of matrix_ion - since the ionization process is only included in a matrix calculation. This is up for debate however.

Higginbottom commented 4 years ago

Run the regression tests - only difference other than noice is in agn_short:

agn_short A little extra feature - OK, so all seems fine.

jhmatthews commented 4 years ago

This is really great, Nick.

Just to check, I'm guessing this is not included for macro-atoms - is it included for simple-atoms in indivisible packet mode? i.e. what would happen to the metals if I use line transfer mode macro_atoms_thermal_trapping and the h10_standard80 dataset.

Higginbottom commented 4 years ago

you guess correctly!

All it does is work out the ‘correct’ rates in the matrix ionization calculation, and then add in the heating implied by those rates. There are no ‘protective’ statements around these rates in the matrix ion calculation.

Do we need to worry about macro atom mode - this is ‘kind of’ like adiabatic heating - it isnt anything to do with photons…. Showing my ignorance of what really goes on in macro atoms again here…

Nick

On 29 May 2020, at 10:52, James Matthews notifications@github.com wrote:

This is really great, Nick.

Just to check, I'm guessing this is not included for macro-atoms - is it included for simple-atoms in indivisible packet mode? i.e. what would happen to the metals if I use line transfer mode macro_atoms_thermal_trapping and the h10_standard80 dataset.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/agnwinds/python/issues/447#issuecomment-635883841, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAZMZPMUUEWUBZK4TVUSN73RT6AUTANCNFSM4GBYWFMQ.

jhmatthews commented 4 years ago

We might need to think about this a bit, since one can still go through the same matrix ionization mode for macro-atoms, but the elements being treated as macro-atoms (which normally includes H) have a protective statement so that their ion and level populations come from macro_pops. So I'm not sure how that will work in tandem with the fixing of Hydrogen.

Also, regarding the number density of hydrogen, it seems to me that this might be problematic, if we only require an accuracy of 0.03 or whatever FRACTIONAL_ERROR is on Hydrogen.

Higginbottom commented 4 years ago

Cool - I’d certainly welcome any input on this.

I’m a little anxious myself about fixing hydrogen - but I cant see any way around it. Originally I thought that the matrix scheme would allow one to do this ‘properly’ but of course whilst any rate can depend on the density of any ion, it cant depend on two ions.

Good to discuss this on Monday..

N

On 29 May 2020, at 11:04, James Matthews notifications@github.com wrote:

We might need to think about this a bit, since one can still go through the same matrix ionization mode for macro-atoms, but the elements being treated as macro-atoms (which normally includes H) have a protective statement so that their ion and level populations come from macro_pops. So I'm not sure how that will work in tandem with the fixing of Hydrogen.

Also, regarding the number density of hydrogen, it seems to me that this might be problematic, if we only require an accuracy of 0.03 or whatever FRACTIONAL_ERROR is on Hydrogen.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/agnwinds/python/issues/447#issuecomment-635889173, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAZMZPPTFT6HEIUOHDDGQ2TRT6CB7ANCNFSM4GBYWFMQ.