Photons absorbed in macro atom model of a dense CV wind

[kslong]

Photons are not normally supposed to be absorbed in macro atom models, and yet in a dense CV wind that is exactly what is happening. These photons are being reported as absorbed and in 87d they appear carry a small amount of energy, but in fact the value that is being reported is not the original energy but the the final energy carried by the photons that is being reported.

foo.pf.txt

[kslong]

This is part of a print out of some variables atter an ionization cycle.

The lines labelled stat refer to status of photons after the transport, e.g escaped 2, absorbed = 6. What is reported is the number of photons, the final weight of those photons, and the initial weight of the photons.

The lines labelled rad, refer to emitted photons. but here the numbers are associated with the origin (and whether they are macro atom photons are not.

The main point of this is to show that the absorbed photons actually represent a large fraction of the total luminosity of the system.

!!python: photon transport completed in: 98.738525 seconds
xxx Sun Sep 24 14:59:25 2023    106.6 NOK                  Photon transport completed
XXX stat        0            0         0.000e+00       0.000e+00
XXX stat        1            0         0.000e+00       0.000e+00
XXX stat        2       282086         2.247e+35       2.753e+35
XXX stat        3            0         0.000e+00       0.000e+00
XXX stat        4            0         0.000e+00       0.000e+00
XXX stat        5            0         0.000e+00       0.000e+00
XXX stat        6       134580         6.060e+28       1.637e+35
XXX stat        7            0         0.000e+00       0.000e+00
XXX stat        8            0         0.000e+00       0.000e+00
XXX stat        9            0         0.000e+00       0.000e+00
XXX stat       10            0         0.000e+00       0.000e+00
XXX stat       11            0         0.000e+00       0.000e+00
XXX stat       12            0         0.000e+00       0.000e+00
XXX rad        0        4.278e+32       4.284e+32
XXX rad        1        0.000e+00       0.000e+00
XXX rad        2        1.781e+35       1.785e+35
XXX rad        3        0.000e+00       0.000e+00
XXX rad        4        0.000e+00       0.000e+00
XXX rad        5        0.000e+00       0.000e+00
XXX rad        6        0.000e+00       0.000e+00
XXX rad        7        0.000e+00       0.000e+00
XXX rad        8        0.000e+00       0.000e+00
XXX rad        9        0.000e+00       0.000e+00
XXX rad       10        1.142e+32       2.199e+32
XXX rad       11        0.000e+00       0.000e+00
XXX rad       12        4.602e+34       9.611e+34
XXX rad       13        0.000e+00       0.000e+00
XXX rad       14        0.000e+00       0.000e+00
XXX rad       15        0.000e+00       0.000e+00
XXX rad       16        0.000e+00       0.000e+00
XXX rad       17        0.000e+00       0.000e+00
XXX rad       18        0.000e+00       0.000e+00
XXX rad       19        0.000e+00       0.000e+00

[kslong]

@jhmatthews @lazygun37

The problem is occurring here

Hardware watchpoint 9: pp.w

Old value = 3.7634096718185621e+26
New value = 7.4169601326339181e+25
scatter (p=0x7fffffffd750, nres=0x7fffffffd6f4, nnscat=0x7fffffffd6f0) at resonate.c:1102
1102              macro_gov (p, nres, 2, &which_out);   //routine to deal with kpkt
2: prob_kpkt = 0.19708085963040745
3: p->w_orig/p->w = 166.94810355415294
4: *p = {x = {-2962297871.4855661, 1874136098.8262272, 647767978.13935328}, lmn = {0.39856551291945708, -0.64088882677881231, -0.65605414686695429}, freq = 16545003481571256, freq_orig = 968201694504242.75, w = 7.4169601326339181e+25, w_orig = 1.2382474282799913e+28, tau = 0, istat = P_SCAT, frame = F_LOCAL, nscat = 6, nrscat = 0, nmacro = 5, nres = 200132, line_res = -3, nnscat = 1, grid = 60, origin = PTYPE_STAR_MATOM, origin_orig = PTYPE_STAR, np = 1, path = 2940407395.672895, ds = 57694.811732947303}
(gdb) where
#0  scatter (p=0x7fffffffd750, nres=0x7fffffffd6f4, nnscat=0x7fffffffd6f0) at resonate.c:1102
#1  0x000055555557b100 in trans_phot_single (w=0x55559f896630, p=0x7fff872e70a8, iextract=0) at trans_phot.c:430
#2  0x000055555557a56c in trans_phot (w=0x55559f896630, p=0x7fff872e7010, iextract=0) at trans_phot.c:131
#3  0x00005555555eb4a6 in calculate_ionization (restart_stat=0) at run.c:225
#4  0x0000555555559b70 in main (argc=2, argv=0x7fffffffe368) at python.c:653

The problem can be avoided by turning off "upweighting", which here is downweighting photons regularly.

[kslong]

The "upweigthing" for bf is carried out a few lines earlier in the code.

else if (phot_top[*nres - NLINES - 1].macro_info == FALSE || geo.macro_simple == TRUE)
      {
        /* Simple ion case.  It's a bf interaction in a calculation involving macro-atoms, 
           but this photoionization x-section is is associated with one of the levels
           a simple atom, not one that is part of a full blown macro atom.

           (Alternatively we are treating all atoms in a simplfied mode)

           Need to make decision about making a k-packet. Get the fraction of the energy
           that goes into the electron rather than being stored as ionisation energy: this
           fraction gives the selection probability for creating a k packet. It's given by 
           (1 - ionization edge frequecy / photon frequency)

           Note:  In some cases, one can obtain a negative prob_kpkt below.  This happens because
           we use the  mean frequency along the path to calculate continuum opacities
           in calculate_ds, whereas here we have the exact frequency at a particular point.
           This leads to a situation where the co-moving frequency can be less than
           the edge frequency.  See issues #436 and #867.

           If the error is large, it suggests that the length of that a photon is allowed
           to travel in a single step is too large.

         */

        prob_kpkt = 1. - (phot_top[*nres - NLINES - 1].freq[0] / freq_comoving);

        if (prob_kpkt < 0)
        {
          /* only report an error for a negative prob_kpkt if it's large-ish in magnitude. see #436 discussion */
          if (prob_kpkt < -1e-2)
          {
            Error ("scatter: kpkt probability (%8.4e) < 0 for phot_top %d, zeroing\n", prob_kpkt, *nres - NLINES - 1);
            Log ("scatter: photon edge frequency: %8.4e, comoving (observer) frequency %8.4e %8.4e\n", phot_top[*nres - NLINES - 1].freq[0],
                 freq_comoving, p_orig.freq);
          }
          prob_kpkt = 0.0;
        }

As written it always reduces the photon weight. How is the energy supposed to come out @jhmatthews

[jhmatthews]

See the docs here: https://agnwinds.readthedocs.io/en/dev/physics/macroatoms.html?highlight=upweight#bound-free-continua-of-simple-atoms - there's another place where we upweight.

Note also that we record and report the amount of energy being extracted/deposited in the bound-free simple ion pool. I think the two bits of code are in kpkt() and scatter() -- although one of the comments is a bit confusing there and seems to have got mixed up.

[kslong]

@jhmatthews The situation described in Readthedocs is not what is happening. In this case, the problem is that the photons are being downweighted so much from multiple scatterings that they eventually count as absorbed. Note that I have tried Christian's suggestion for running a model to convergence (75%) without upweighting, and then running another cycle with it on. While we don't loose 50% of the flux as we do when you run with upweighting on, we do lose 25%, which I think should be concerning (even if we make not upweighting the default.)

Is there anywhere that we have a full description of how this is all supposed to work; the readthedocs page really does not explain this anywhere. If this does not exist, it would be extremely useful if one were created.

[jhmatthews]

I posted the Readthedocs page because you asked above how the energy is supposed to come out, not because it explains this issue which has only just come up in this particular form. The cautionary note about numerical problems is indeed a different issue (although related).

In terms of how all this is supposed to work, the description on RtD, whilst doubtless flawed, is my best attempt to explain that. Perhaps there are some things missing but I can't think what exactly.

[kslong]

There was a suggestion made by @ssim that this problem might be due to the fact that we do not use the Milne relation to calculate recombination rates might result in the problem we are seing, because the recombination rates might be higher than predicted. In that case, one might eliminate the problm if we used the Milne relation instead. @jhmatthews suggested we could test this by removing the separate recombination rates from the program. So I modified the data file we have been using h20_hetop_standard80.dat, removing these lines '''' data/atomic/recomb.dat data/atomic/chianti_dr.dat data/atomic/chianti_rr.dat data/atomic/Badnell_GS_RR.dat data/atomic/di_dere.dat '''' but this had no effect. We are still losing a lot of energy into the ion pool

The logs say that the Milne relation is being used, so presumably this is a good test of this hypothesis.

[kslong]

We document the flow of energy in the log files for the entire wind. For the original model and the one which is supposed to follow the Milne relation this is what one gets:

!! report_bf_simple_ionpool: Total flow into: 1.0706e+35 and out of: 3.0572e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 2.0153e+35 and out of: 9.7539e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9328e+35 and out of: 5.5755e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9107e+35 and out of: 3.7056e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9225e+35 and out of: 3.0034e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9375e+35 and out of: 2.5415e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9348e+35 and out of: 2.3325e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9404e+35 and out of: 2.1314e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9451e+35 and out of: 2.0195e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9428e+35 and out of: 1.9931e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9571e+35 and out of: 1.8773e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9604e+35 and out of: 1.9484e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9579e+35 and out of: 2.0849e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9628e+35 and out of: 2.0314e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9689e+35 and out of: 2.1613e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9715e+35 and out of: 2.0468e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9681e+35 and out of: 2.1715e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9849e+35 and out of: 2.1706e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9851e+35 and out of: 2.1656e+34 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.9806e+35 and out of: 2.5205e+34 bf_simple ion pool

In this case the energy in is typically 10 x the energy out. Note there is one line for each of the 20 cycles

BUT if one runs a macro atom model, without any macro atoms you get something quite different

!! report_bf_simple_ionpool: Total flow into: 1.7142e+36 and out of: 1.6677e+36 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 9.2296e+35 and out of: 7.4214e+35 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.2051e+36 and out of: 1.0200e+36 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 4.9550e+38 and out of: 4.9610e+38 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 2.8208e+40 and out of: 3.1459e+40 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.0955e+41 and out of: 5.4736e+40 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 2.9381e+39 and out of: 3.1784e+39 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 6.0569e+37 and out of: 6.5355e+37 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 5.5339e+36 and out of: 6.0763e+36 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 3.6363e+38 and out of: 3.6920e+38 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 2.7500e+37 and out of: 3.2789e+37 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.0501e+37 and out of: 1.1638e+37 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 6.9131e+37 and out of: 6.7843e+37 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 2.1363e+37 and out of: 3.3468e+37 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 1.5875e+36 and out of: 2.9579e+36 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 2.3165e+36 and out of: 2.2220e+36 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 7.4902e+36 and out of: 7.3923e+36 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 2.2931e+37 and out of: 2.2221e+37 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 2.3686e+37 and out of: 2.3999e+37 bf_simple ion pool
!! report_bf_simple_ionpool: Total flow into: 3.1149e+37 and out of: 3.1298e+37 bf_simple ion pool

The totals are higher now, and there is variability. But in this case the energy in is equal to the energy out.

@jhmatthews @ssim Is this expected? Is it possible that something is wrong when a simple atom is photoionized but a macro atom is selected for recombination.

[kslong]

@jhmatthews On branch bug1024_lost, I added more tracking of what bf process are doing. The file shows (I think) what fb process have absorption and emission. The first column is the fb no, the second is the elment involved, the third is the ion state and the 4th and 5 columns are the number of times one absorbed by that fb process and the number of decays by that process.

Some of it makes sense, we see lots of re-combinations to H and He. What does not make sense is the huge number of interactions with NV and OVI, since both are almost absent in the wind.

foo.txt

[jhmatthews]

I agree that's confusing. I guess it's possible that there are some sizable minority of cells where there is a decent amount of NV and OVI, since I think this file is a cumulative one across the whole wind? It's not clear to me exactly what we mean here by "number of times one is absorbed by fb process" i.e. whether it means r-packet activations of a simple-atom via that bound-free process? If so, the chance of it happening is set just by the bound-free opacity, I think.

The post before about macro-atoms is less surprising to me, because removing the macro-atom data is changing how the H and He continua are treated, so I might well expect the problem to increase (the downweighting can also take place in the H and He continua if they are treated as simple ions, rather than being done "self consistently" within the macro-atom).

[kslong]

The numbers here are, as James guessed the total for the wind, but as far as I can tell the OVI and NV concentrations are very very small. If one does a diff between the bug1024_lost branch, it is easy to see where I did the counting.

I guess my question about the difference between results with and without macro atoms is for the H and He macro atom case, do these create a situation where necessarily we are creating an imbalance. Photons that undergo bound free absorption necessarily lose the ionization energy, but when they come back through a macro atom recombination, they may not get it back. I know the upweighting which occurs early is expected to compensate for that, but it does not seem to.

[kslong]

Well, my previous comment was incorrect. Indeed OVI and NV are the dominant versions of these species, as the disk in this model is very hot 10**5 K. So it rather looks like the following is happening. BF transitions are dominated by ionization of ions like NV and OVI, but they are recombining mainly into H and He. Photons are losing ionization energy when they are absorbed, but the boosting that occurs is insufficient to make up for the energy that they lose.

The question that needs answering it seems to me now is why the upweighting does not, on average at least, balance this all out.

[ssim]

Hi - maybe we can discuss more if we all talk later. But it sounds like what we're getting is that photoionzation is mainly absorbing photons via the metals, but recombination is mainly going through hydrogen and helium? If that's true as a statement about photons then it must mean that there's an inconsistency between the actual ion populations and the rates of photon destruction / emission somewhere.

Do we know if this issue occurs if we don't include any metals at all? I.e. if we use only H+He macro atoms (and nothing else) does the same problem occur? Or is this some interplay between the two.

One thing that comes to mind is that we probably do expect to lose energy in a scheme with the upweighting by ionisation energy - because that ignores any energy that comes via the radiative bound bound cascade in 'simple' atoms.

On 3 Oct 2023, at 22:17, Knox S. Long @.***> wrote:

This message is from an external sender. Please take care when responding, clicking links or opening attachments.

Well, my previous comment was incorrect. Indeed OVI and NV are the dominant versions of these species, as the disk in this model is very hot 10**5 K. So it rather looks like the following is happening. BF transitions are dominated by ionization of ions like NV and OVI, but they are recombining mainly into H and He. Photons are losing ionization energy when they are absorbed, but the boosting that occurs is insufficient to make up for the energy that they lose.

The question that needs answering it seems to me now is why the upweighting does not, on average at least, balance this all out.

— Reply to this email directly, view it on GitHubhttps://github.com/agnwinds/python/issues/1024#issuecomment-1745739439, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AAGA6LNHCI65T2WQ3BF2EVDX5R6E7AVCNFSM6AAAAAA5FFBMPCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTONBVG4ZTSNBTHE. You are receiving this because you were mentioned.Message ID: @.***>