Dealing with metal fields

[brittonsmith]

The addition of non-equilibrium metal species and dust evolution models will require us to make some decisions about how metal-related fields are defined and some other expectations of behavior.

The Current Behavior

The Grackle main branch essentially only models gas-phase chemistry right now. We only allow for a passive dust model where the dust field is not considered to contribute to the total density. Specifically, we have the total gas density, density, defined as:
density = metal_density + $\sum_i$ primordial_density_i,
where primordial_density_i represents the mass densities of HI, HII, HeI, etc.
This affords us the luxury of guaranteeing that density remains constant over the course of a call to solve_chemistry

The Single-Grain Dust Model

With the new single-grain dust model, we will now be tracking a dust_density field that will exchange mass with the metal_density field as dust grains grow or are destroyed. Now we face a cross-roads with two possible paths:

Option 1: density = metal_density + dust_density + $\sum_i$ primordial_density_i
This has the benefit of continuing to guarantee that density remains constant, but has the disadvantage of requiring several minute alterations all over the plus, like changing how we define metallicity and how we define the primordial component of the density field.

Option 2: density = metal_density + $\sum_i$ primordial_density_i
Here, we have the opposite situation. After speaking with @aqua-hl I think this is better option. It also allows us to maintain that the density field is the one used for hydrodynamics.

The Multi-Grain Dust Model

Using the new multi-grain model will require the user to also supply injection metal fields (see #444) to trace the map metals coming from specific creation pathway to yields of gas-phase metals and dust species. This introduces another instance where the path forward is not obvious. The injection pathways give the fraction of metal mass of going into seven elements and up to 13 different dust species (see #465). However, in no configuration do the gas and dust yields fraction for the total metal yield. Below, I show the total gas and dust yields for each injection pathway. Note, if fewer density species are used, the dust-phase and absolute totals are even lower.

local_ISM:

• gas-phase total: 0.33537661799999996
• dust-phase total: 0.5061552
• total: 0.841531818

ccsn13:

• gas-phase total: 0.698820896
• dust-phase total: 0.09255480452
• total: 0.79137570052

ccsn20:

• gas-phase total: 0.8700132180000001
• dust-phase total: 0.055497320659999995
• total: 0.9255105386600001

ccsn25:

• gas-phase total: 0.809653515
• dust-phase total: 0.11757425859999998
• total: 0.9272277736000001

ccsn30:

• gas-phase total: 0.9064142779999999
• dust-phase total: 0.039527683765792514
• total: 0.9459419617657924

fsn13:

• gas-phase total: 0.8244579428270001
• dust-phase total: 0.17554200000000802
• total: 0.9999999428270081

fsn15:

• gas-phase total: 0.8107713175090001
• dust-phase total: 0.1892290000000369
• total: 1.000000317509037

fsn50:

• gas-phase total: 0.9998899782668
• dust-phase total: 0.0001098505113741244
• total: 0.9999998287781742

fsn80:

• gas-phase total: 0.9913192637702
• dust-phase total: 0.00868025000275811
• total: 0.999999513772958

pisn170:

• gas-phase total: 0.92304398
• dust-phase total: 0.013309690847284955
• total: 0.9363536708472849

pisn200:

• gas-phase total: 0.922013672
• dust-phase total: 0.0009824352838940825
• total: 0.922996107283894

y19:

• gas-phase total: 0.25894751
• dust-phase total: 0.5
• total: 0.75894751

How do we track the addition metal mass not being mapped to a specific species? Not fully understanding what is needed by the multi-grain growth functionality, I think we can either put untracked metal mass into the metal_density field or we need an extra metal field. This brings me to the final (and possibly identical) point.

Tracking Gas-Phase Metals

The non-equilibrium metal fields that we will be able to follow will of course be a subset of all possible species. We will likely need to supplement the non-equilibrium metal cooling with metallicity dependent tables. At high temperature, these tables will have to account for the cooing from all metals (including those we are tracking) as we have not specifically implemented high temperature transitions. At low temperature, they will only complement the non-equilibrium cooling. Potentially, the low temperature regime could be done another way, but we'll need to think about it.

Minimally, we need to understand what metal_density will still mean. It should either represent the sum total of all metals (including the ones we track) or it should represent only the metals we are NOT tracking. I am not confident I have completely mapped out the ramifications of either choice yet.

[mabruzzo]

There's a typo in your post:

However, in no configuration do the gas and dust yields fraction for the total metal yield.

To clarify, are the numbers in your table tracking the fractions of the injection-pathway-metal-density field? (I believe that these are the primary fractions tracked within the data-tables in Cholla)

[aqua-hl]

On the “Tracking Gas-Phase Metals” question, I think the same reasoning we used for Option 2 in the single-grain dust model applies here as well.

If we redefine metal_density to mean “only the metals we are NOT tracking” (i.e., let it stand alone as an “other metals” field once we start tracking specific species like C and O), we’d run into the same problem we discussed for the dust case: every piece of downstream code that currently assumes metal_density represents the total metal content would need to be updated. Anything computing metallicity, cooling rates, or using metal_density in logical checks would have to be revisited and changed to sum metal_density together with whatever tracked species are active, and that coupling would grow every time we add a new tracked metal.
I have a working version where I went down this path and made metal_density represent only the untracked metals, and in practice I found there were quite a lot of places that needed to be changed to keep everything consistent, more than I initially expected. That experience makes me more wary of this direction.
By analogy with the density argument, I’d lean toward keeping metal_density as the sum total of all metals (including the ones we track individually). That way metal_density keeps a stable, consistent meaning across configurations, and the rest of the code doesn’t need to branch on which specific metals happen to be tracked in a given run. The tracked species fields would then be treated as a breakdown/subset of metal_density, similar to how primordial species relate to the total gas density.

Or, alternatively, maybe we just introduce a new field called total_metal_density, that way metal_density can keep whatever meaning we land on, and we have a dedicated field that unambiguously represents the sum of all metals for the code paths that need it.

[mabruzzo]

If we redefine metal_density to mean “only the metals we are NOT tracking” (i.e., let it stand alone as an “other metals” field once we start tracking specific species like C and O), we’d run into the same problem we discussed for the dust case: every piece of downstream code that currently assumes metal_density represents the total metal content would need to be updated.

When you say "downstream code," are you talking about external simulation codes? I'm not sure I follow what they'll need to change...

I'm also not sure what you mean by the dust case. I don't think we currently require a total dust density field if we are tracking individual species. Would there be any value to tracking such a field? (I can't think of one)

Anything computing metallicity, cooling rates, or using metal_density in logical checks would have to be revisited and changed to sum metal_density together with whatever tracked species are active, and that coupling would grow every time we add a new tracked metal.

Right, so I'm not totally sold on one direction or another here... Can we talk through the points you are making here (I suspect some of the highlighted code probably needs to change one way or the other):

• for "using metal_density in logical checks," are you just describing make_consistent? Or something else? For make consistent, the coupling probably needs to grow anyway with either approach? Right? Plus make_consistent has always been a bit of a numerical "kludge"
• we use metallicity for cooling rates (more on that in a moment) and to maybe skip some logic. Frankly, I'm not entirely sold on using metallicity alone for the latter point, but you point is definitely valid here
• in terms of cooling function -- we need to give this some serious thought. Any long-term solution "correct" solution I can think of probably needs metal_density minus all explicitly tracked species. (Or do you have something else in mind?)

I guess this will also complicate calculation of mean molecular weight

[mabruzzo]

It's worth mentioning that mass transfer between gas and dust is going to introduce a lot of complexity when it comes to density. We're going to need to stop using HydrogenFractionByMass in these configurations and we'll need to totally change our Deuterium accounting.

Perhaps most importantly, we probably need to be very careful about how exactly we choose to update mass density. I have some concerns about numerical stability on this point

[aqua-hl]

If we redefine metal_density to mean “only the metals we are NOT tracking” (i.e., let it stand alone as an “other metals” field once we start tracking specific species like C and O), we’d run into the same problem we discussed for the dust case: every piece of downstream code that currently assumes metal_density represents the total metal content would need to be updated.

When you say "downstream code," are you talking about external simulation codes? I'm not sure I follow what they'll need to change...

I meant Grackle-internal code paths, not external simulation codes. The public field interface doesn't have to change either way, what I was worried about is the set of callers inside src/clib/ that read metal_density: cooling-rate lookups, calc_temp*, make_consistent, the mean-molecular-weight computation, and so on. In my working branch where I redefined metal_density as "untracked metals only," every one of those sites needed to be audited and several needed to start summing in the tracked species to reconstruct the total. That coupling grows every time we add another tracked element.

I'm also not sure what you mean by the dust case. I don't think we currently require a total dust density field if we are tracking individual species. Would there be any value to tracking such a field? (I can't think of one)

You're right that we don't carry a total_dust_density alongside per-species dust. What I was reaching for is there was a time we tried to let density represent gas + metal + dust, but turns out there are quite a lot to change if we do that, and we just decided to keep gas density and dust density separate.

Anything computing metallicity, cooling rates, or using metal_density in logical checks would have to be revisited and changed to sum metal_density together with whatever tracked species are active, and that coupling would grow every time we add a new tracked metal.

Right, so I'm not totally sold on one direction or another here... Can we talk through the points you are making here (I suspect some of the highlighted code probably needs to change one way or the other):

• for "using metal_density in logical checks," are you just describing make_consistent? Or something else? For make consistent, the coupling probably needs to grow anyway with either approach? Right? Plus make_consistent has always been a bit of a numerical "kludge"

Mostly just make_consistent. I agreed it's already a numerical kludge and the coupling probably has to grow either way.

• in terms of cooling function -- we need to give this some serious thought. Any long-term solution "correct" solution I can think of probably needs metal_density minus all explicitly tracked species. (Or do you have something else in mind?)

On cooling rates, I think you're right. Long-term-correct cooling path needs metal_density − sum(tracked species) at the point where Cloudy tables (parameterized by total Z) are combined with network-resolved C/O cooling, otherwise we double-count. So the subtraction happens regardless of which convention we pick for the field itself. Given that, I'm thinking between:

Keep metal_density = total, and compute metal_density − sum(tracked) locally inside the cooling routine.
Introduce a dedicated untracked_metal_density field so the subtraction is explicit in the data model, and metal_density stays total.

[aqua-hl]

Perhaps most importantly, we probably need to be very careful about how exactly we choose to update mass density. I have some concerns about numerical stability on this point

All fair. Whatever we decide here, the gas↔dust mass bookkeeping has to be done carefully, and those are genuinely the headaches.

[brittonsmith]

I'm only going to add a though/realization that I had a couple weeks ago. I've shared this in person, but it's worth recording for posterity. The process of preparing PR #536 changed my original opinion to aligning more with the discussion here, namely that metal_density should track all gas-phase metals, including the ones we follow explicitly. This will make it far simpler to create supplementary cooling tables since we can use the tables to represent not simply the species we are not following, but the coolants we are following. In the past, I have made metals-only cooling tables by doing a component by component subtraction with a metal-free table. We can do the same here by first identifying the processes we are solving for explicitly, then removing them from the full run. In this case, metal_density will need to represent all metals.