dust/gas heat transfer cooling rate can be unphysical

[brittonsmith]

Setting either h2_on_dust or dust_chemistry to 1 activates a heating/cooling process to represent the heat exchange between gas and dust through inelastic collisions. This term has a factor of (Tgas - Tdust) to make the heat exchange proportional to the temperature difference between gas and dust. At very high gas temperatures, this can lead to enhanced cooling rates that are likely unphysical. Below is a sample cooling function for a gas at solar metallicity with heating from a UV background enabled.

Here is the same cooling function with h2_on_dust=1.

The increased cooling at temperatures above a few million K is entirely this term. In reality, dust would probably not exist long in an environment this hot. However, since we do not calculate dust destruction (although this is coming eventually), Grackle computes a cooling rate assuming the dust is still present. There are valid use cases that include the effects of dust without actively evolving the dust population (e.g., assuming the dust to gas ratio simply scales with metallicity). For example, if you just want to capture H2 formation on dust. I think it's worth trying to limit this unphysical behavior for this case by applying a temperature ceiling above which we do not include this cooling term (just for this passive dust model). I think a temperature or 1e4 K or 1e5 K would be reasonable.

[scog1234]

For what it's worth, I include a high temperature cutoff like this in the treatment of dust cooling in my chemical module, in order to avoid the unphysical high T cooling you point out here. What I have found works well is multiplying the cooling term by an exponential factor exp(-T/T_cut) with some large value for T_cut (I use 3e4 K), so that the rate falls off gradually at high T rather than sharply cutting off.

[brittonsmith]

Thanks, @scog1234, this is a very sensible solution. I'll submit a PR for this.

[mabruzzo]

When you make this PR, would you mind working off the branch in #486? In that PR I factored out the dust heating/cooling functions into a separate file

ASIDE: that PR doesn’t actually have any dependencies. We could merge it in now if that helps

[brittonsmith]

@mabruzzo, my thought was that I could apply this directly to the calculating of the gasgr rate in rate_functions.c. I think that would be the most minimally invasive change since I wouldn't then have to apply the exponential dampener in both cool1_multi and the grain temperature calculation. I was hoping to do this in the main branch as its a relevant bug in some other work. I would hope that this could be ported over to newchem-cpp easily. We'd then probably have to apply it to gasgr2 separately, but that might be safer.

[brittonsmith]

Resolved in PR #503.