Inconsistency in UVB rates with self-shielding

[brittonsmith]

Originally reported by Andrew Emerick (Bitbucket: aemerick, GitHub: aemerick)

---

I've uncovered a fairly significant issue when including H_2 formation with the self-shielding prescription turned on that leads to runaway H_2 formation at number densities above about n = 10 cm^-3 (the exact point is somewhat metallicity dependent).

The issue boils down to the fact that the rate coefficients from the UVB are split into independent bands (k24 - k31). Although we account for shielding of HI, HeI, and HeII radiation in (k24, k25, k26) with the Rahmati et. al. prescription, and account for H2 self-shielding against LW radiation (k31), the rest of the bands are not modified in any fashion.

This becomes a problem in self-shielding regions that start forming molecular hydrogen. Some of this H2I becomes ionized to H2II through rate k29, generating free electrons that go into making H-, which in turn makes more H2. The H2II quickly collides with the abundant neutral HI, making H2I and HII. The net result of this is a continual source of free electrons (H2I->H2II) that drives further H2 formation via the gas phase. This spikes the electron fraction to ~percent levels, and the H2I fraction to ~50-70% (by mass).

The reality, however, is that direct H2I ionization (15.2 eV) should be shielded significantly by the neutral HI. In addition (though not related to this problem) is the H2II dissociation radiation (30 eV) should also be shielded, mostly by HeI shielding.

I propose implementing shielding of these bands (rates k29, k30, and k28) using the same scaling provided by the Rahmati et. al. presectption, attenuating the 15.2 eV band using the same reduction factor as HI, and attenuating the 30 eV band using the same reduction factor as HeI shielding.

I'm including two plots showing the evolution of all H mass fractions over time in a one-zone cooling cell model, with an initial number density of 10, low metal mass fraction (2% solar), including H2 self-shielding from LW, and with ``self_shielding_method = 3''. The top is what happens with our current settings, the bottom shows what happens when including the shielding as I proposed above.

---

• Bitbucket: https://bitbucket.org/grackle/grackle/issue/17

[brittonsmith]

Original comment by Simon Glover (Bitbucket: scog, GitHub: scog)

---

I think this makes sense as a short-term fix, but we might want to document this somewhere as something to revisit in the future.

[brittonsmith]

Original comment by Britton Smith (Bitbucket: brittonsmith, GitHub: brittonsmith)

---

Nice catch. Your proposal seems reasonable to me, but maybe @scog or @gbryan have other opinions?

[brittonsmith]

I just happened to be looking at solve_rate_cool.F and it looks like Andrew's proposed solution was implemented for self_shielding_method=2,3. Specifically, I see this:

!     Scale H2 direct ionization radiation
            if (k29 .lt. tiny8) then
              k29shield(i) = 0._DKIND
            else
              k29shield(i) = k29shield(i)*f_shield_H(i)
            endif

!
!     Apply same equations to HeI (assumes HeI closely follows HI)
!

            if (k26 .lt. tiny8) then
               k26shield(i) = 0._DKIND
            else
               k26shield(i) = k26shield(i)*f_shield_He(i)
            endif

!     Scale H2+ dissociation radiation
            if (k28 .lt. tiny8) then
                k28shield(i) = 0.0_DKIND
            else
                k28shield(i) = k28shield(i)*f_shield_He(i)
            endif

            if (k30 .lt. tiny8) then                 
                k30shield(i) = 0.0_DKIND
            else
                k30shield(i) = k30shield(i)*f_shield_He(i)
            endif

@gregbryan, maybe Andrew did something other than this that never made it back into the code?

@leonardromano, I imagine you are using a version of the code that contains this fix, but maybe not? Maybe you can comment on the solutions you've adopted.

[gregbryan]

Thanks @brittonsmith -- I also came to that conclusion shortly after our call. I think these are the changes that resolved the issue (discussed briefly in https://arxiv.org/pdf/1807.07182.pdf). If @leonardromano can check if this addresses the issue he mentioned during the call today, we should add documentation reflecting this change.

[leonardromano]

@brittonsmith My use case was somewhat different, in that in my simulations I had an interstellar radiation field present, that was typically of order 1 to 100 in units of the Habing radiation field, so in the regime where shielding was relevant, the UVB was only marginal. Therefore, the solution that worked for me, was to simply attenuate the "RT" rates for HI and HeI before I pass them to grackle, and leave the UVB rates untouched (as some sort of floor value for the radiation field). This worked for me, but in general I think it doesn't really solve the issue. @gregbryan I will check and let you know, whether or not this did the trick.

[leonardromano]

@gregbryan I checked and tried to reproduce the problem, but I couldn't really reproduce it. It's possible that I was previously using an older version of Grackle, so it might have been just that.

Anyways, I checked with a number of different setups (self_shielding_method=0, 1, 2 & 3; with & without RT rates; 1 run with artificially high metallicity) and I could not find any obvious issues, so maybe it really isn't an issue anymore.

[gregbryan]

Thanks @leonardromano! @brittonsmith - should we close this since it seems like the proposed fix was implemented?

[brittonsmith]

Yes! Let's close it. Thanks for looking into this, @leonardromano.

[leonardromano]

Hi Britton,

As I mentioned during the last meeting, I would be willing to check if the issue as actually been resolved. I have been on vacation in the meantime, but I might start looking into it soon. Has this test already been done in the meantime and that’s why you closed the issue?

Best, Leonard

Am 30.11.2022 um 18:17 schrieb Britton Smith @.***>:

Closed #17 as completed.

— Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you were mentioned.