radiation enhancements

[BenWibking]

We need to modify the radiation-matter exchange to obtain the specific heat capacity $c_v$ from Microphysics.

My mistake, I forgot this was already done here: https://github.com/quokka-astro/quokka/blob/7767642c624ac8ea3765e7ba8107b19078bb2052/src/radiation_system.hpp#L991

[BenWibking]

@psharda @chongchonghe can you coordinate on this?

[markkrumholz]

I actually think we want to have a broader discussion of the interaction of multigroup radiation with microphysics; this should be on the agenda the next time we talk as a group. My background thinking is as follows: right now we do our calculation assuming that the matter source function is the Planck function and that matter has constant c_v. However, in the future we want to be able to generalize this considerably. One generalization, which Ben is proposing here, is to have a specific heat that is not constant, and comes from microphysics. However, there are several other generalizations we need for realistic astrophysical systems. The two specific examples I have in mind are:

• A more realistic three-temperature system for thermal radiation where we do not assume that the gas and the dust have the same temperature, but instead that the dust temperature is in an equilibrium set by the balance between radiative heating and cooling and collisional exchange. In this case the matter temperature is still the only variable to be iterated to convergence, but the emissivity isn't just B_nu(T_gas), it's a more complex function of T_gas.
• Ionization, where the source function is determined by the matter temperature and ionization state, and looks like j_nu = (mean photon energy) (recombination recombination rate coefficient) n_e * n_ion. In this case the source function depends on T_gas and the chemical abundances, which need to be iterated to convergence together.

My thought on how to implement this is that we want to have an enum that specifies how each radiation group is to be treated; the options that I have outlined so far are "2Tconst" (what we're doing now), "2T" (Ben's proposal to get c_v from microphysics and not assume it is constant), "3T" (dust and gas not assumed to be in equilibrium), and "ionization". It is possible that some of these could be collapsed together, since for example 2Tconst is the limiting case of 2T for d(c_V)/dT -> 0, and 2T is the limiting case of 3T for dust-gas collision coefficient -> infinity. Each of these cases generates different entries in the matrix, but it retains the key property it has now which lets us invert it in O(N) time, which is that the different radiation groups talk to each other only through the matter, so each row of the matrix has only two (or if we're iterating on ionization fraction as well, three) non-zero entries. Thus we can use the same basic solver we are using now, and the only thing that really has to change in the code is the bit that generates the matrix entries.

This is just an idea, though, and in general now that @chongchonghe seems to have the 2Tconst case working, I think it is time to come up with a game plan for how this looks going forward.

The other complication we need to discuss is how this interacts with the changes @BenWibking has been talking about to improve our recovery of good behavior for RHD waves in the high-optical depth limit, where our current operator-split treatment doesn't seem to do so well. I think this is relatively simple and it's just an extra RHS source term, but we need to work through it.

[BenWibking]

A few comments:

1. We currently allow $c_v$ to be user-defined -- it does not have to be constant. However, unless the user specifies the function via a template specialization, the code assumes an ideal gas with a constant gamma. Since we now have a full EOS module in Microphysics, we can just use Microphysics to compute $c_v$ (except for a handful of radiation test problems that require an unphysical form of this).

2. This paper (https://arxiv.org/abs/1403.1874) suggested that the differences between 2T and 3T were not significant. Are there test problems where we can demonstrate it makes a difference?

3. The last complication is discussed here: https://github.com/quokka-astro/quokka/issues/33.

[markkrumholz]

Responses:

1. I agree that we can use microphysics to compute c_v, and that we want to do it this way in the future. It's just a question of how we want to stage the changes, and whether we want to do this separately or as part of a broader overhaul of radiation with "complex" microphysics.

2. Not sure why you're citing that paper, since they don't do 3T, and to the extent that they do discuss it they're just quoting the analysis from Krumholz & Thompson (2013). However, that is in the context of driving galactic winds in ULIRGs, an environment where the gas temperature basically doesn't matter because thermal pressure is completely unimportant. For star formation problems where you actually do care about the temperature and thermal pressure because it affects fragmentation, it definitely matters in at least some situations. There is an extensive discussion of this in the method paper by Bate & Keto (2015), and lots in subsequent papers by Bate and by the Starforge group. There are also several test problems in the Bate & Keto paper.

3. Thanks for the reminder. I think we want to discuss how that interacts with a more general coupling between radiation and gas in microphysics.

[BenWibking]

1. I forgot this was already done. See the edited issue description.

2. I was thinking of these sentences: "This κ ∝ T^2 scaling approximately holds for dust in thermal equilibrium with a blackbody radiation field for T ~ 150K (Semenov et al. 2003). This implicitly assumes the characteristic temperature of the radiation and T are equal. Since we evolve the radiation energy density, we can define a characteristic radiation temperature Tr = (Er/ar) 1/4 and check a posteriori if Tr ≃ T . This is generally a good approximation although modest differences are present in some regions, which motivated us to perform alternative simulations with κR,P ∝ T_r^2. The two sets of simulations agreed very closely with each other so we report only the κR,P ∝ T^2 here." In particular, are there opacity tables that do not assume Tr = T that are available? Ok, I see the discussion in Bate & Keto about this. It does not seem great to just switch between Planck and Rosseland means at T_d ~ 100 K, but it's better than nothing.

3. Yes, agreed.

[markkrumholz]

Dropping 1 and 3 and focusing on 2: in a ULIRG it is indeed the case that near-equality is a good assumption, but this is because the density is very high. A rough rule of thumb is that, at Solar metallicity, dust and gas come into equality at a density between around 10^4 and 10^5 cm^-3. In a ULIRG the density is almost everywhere above this threshold, which is why equality is a pretty good assumption. But this is clearly not the case in the Milky Way. This gets to be even more of an issue at low metallicity, because lower dust abundance raises the density required to for equilibration. So if we want to do low-metallicity star formation well, we really do need 3T.

As for opacity tables and what is available, the answer is "sort of" and "it's complicated". There are a few distinct effects occuring in different density and temperature regimes, some of which we will be able to deal with well and some where we'll have more problems. As a start, let me point out that the single most important temperature dependence of opacity comes from the frequency distribution of the radiation field, and that if we have multiple thermal radiation groups we automatically get this even if the opacities themselves have no explicit temperature dependence. Thus for example the thing that drives the kappa ~ T^2 scaling that holds up to ~150 K is that the properties of the dust are basically fixed, so kappa_nu does not depend on T at all, and all of the temperature dependence is coming from the fact that radiation is moving to shorter wavelengths where the quantum efficiency of the grains is higher. So we can get this right just by having a kappa_nu that is independent of dust temperature and letting the radiation solver move energy between photon groups with different nu.

At higher temperatures and densities things get more complicated, which is what Bate & Keto are trying to capture. Here there are two distinct effects: the photon frequency distribution changes, but on top of that the grain properties themselves start to change. For example, the big thing tha happens at 150 K is that ices start evaporating off grains, which lowers their cross section, and this roughly compensates for the quantum efficiency effect, which is why we get kappa ~ const instead of kappa increasing with T. Life is complicated here because the properties of the grains in reality depend on three distinct things: grain temperature (which is always close to radiation temperature) which affects things like the rate of ice sublimation off the grain, gas temperature which affects things like the rate of ice deposition onto the grain and, at higher temperature, the rate of grain sputtering, and density, which affects all the gas-grain coupling rates. So there is no single source of which I am aware that contains tables for all of this stuff and does it all self-consistently. We are going to have to make approximations, and they are only going to work in certain regimes. So for example if we want to do galaxy-scale simulations like @aditivijayan has been doing, and not in the ULIRG regime, the most sensible thing to do is probably just to ignore changes in grain properties and do something like take Bruce Draine's tables of kappa_nu and integrate them over our frequency bins to get an effective opacity in each bin. If we want to do star formation where the gas gets considerably denser and warmer, we might do something like Bate & Keto. Really doing things right probably requires a fully self-consistent grain model.

[BenWibking]

Ah ok, that's the context I was missing. Thanks. Maybe we should make a separate issue for 3T/dust?

[markkrumholz]

Happy to have it be a separate issue on github. I mostly wanted to put it on the agenda for real-time discussion when we next meet, because there are some high-level architecture decisions that it would be good to hash out, for example my proposed approach of "flagging" different photon groups based on what type of source terms they generate in the radiation matrix but then using the same basic solver for all of them.

[BenWibking]

Since we agreed to move this discussion into a separate planning document for radiation enhancements, I'm going to close this issue for now.