Element diffusion is artificially creating new species when diffusion_use_full_net = .false.

[simonguichandut]

We have been experimenting with element diffusion in neutron star atmospheres. Starting with a model of a pure iron layer at the bottom, and a homogeneous layer of carbon and hydrogen above, we would simply expect the carbon and hydrogen to separate out. This is indeed what happens when diffusion_use_full_net = .true..

However when .false., we find that other species in the network are being created in non-negligible amounts. The network is cno_extras_plus_fe56.net, but burning is explicitely turned off with max_abar_for_burning=-1.

Here's an animation of the abundance evolution from pgstar:

Notice new mass fractions appearing at the initial iron boundary. We can also plot the evolution of the total mass of each species:

While the total mass stays constant (black line), we see that iron and carbon-12 are being destroyed, meanwhile every other CNO species is being created. The initial mass of the carbon/hydrogen layer is ~6e-16 Msun, so 1e-17 mass variations are quite large. Also note that while helium is in the network, it is not being created. This is definitely not CNO burning!

Could this be a bug in the way MESA "unwraps" the different classes it considers for diffusion? Apologies if I've missed an important control or something even more trivial. But something seems wrong!

To Reproduce
This is being run on a recent clone of this repo. I'm including a minimal working directory in this zip file:
diffusion_bug.zip

• inlist_diffusion
• inlist_pgstar makes the abundance plot
• history_columns.list (just need add_total_mass for the plot)
• plot.py makes the mass change plot above
• env_c12_h1.mod is the starting model. From the test suite's make_env model, I just accreted a c12/h1 mixture.

[Debraheem]

Thanks for documenting this so well! I was able to reproduce your plot with the files you shared. I also checked with diffusion_use_full_net = .true. which yields:

[evbauer]

I do think you're right that there's a bug here, but it looks like there's a fairly easy workaround that should be sufficient for what you're doing. If you aren't familiar with the way diffusion classes and representatives are set up, I'd recommend taking a quick look at the defaults here.

Diffusion lumps things into classes according to a hierarchy of atomic weights, and you'll notice that the defaults are quite minimal in terms of the number of classes. In particular, c12 doesn't even get it's own class, and the result is that it will get lumped in with the "class" represented by o16. Similarly, everything heavier than o16 will get lumped in with fe56 under these defaults.

I think the "bug" here arises because diffusion has some fixup routines that take small errors in mass conservation and does some renormalization to make sure that mass fractions always add up to 1. When things are lumped together into classes, the fixup corrections could be applied to everything in the net represented by a particular class, and so you might get some non-conservation of species.

However, you can alter your setup to have a more reasonable set of classes without having to go all the way to using the full net, and it looks like things go much better in that case. Here's something that looks like it works well for me:

    diffusion_num_classes = 9

    diffusion_class_representative(1) = 'h1'
    diffusion_class_representative(2) = 'he3'
    diffusion_class_representative(3) = 'he4'
    diffusion_class_representative(4) = 'c12'
    diffusion_class_representative(5) = 'n14'
    diffusion_class_representative(6) = 'o16'
    diffusion_class_representative(7) = 'ne20'
    diffusion_class_representative(8) = 'mg24'
    diffusion_class_representative(9) = 'fe56'


    diffusion_class_A_max(1) = 2
    diffusion_class_A_max(2) = 3
    diffusion_class_A_max(3) = 4
    diffusion_class_A_max(4) = 12
    diffusion_class_A_max(5) = 14
    diffusion_class_A_max(6) = 16
    diffusion_class_A_max(7) = 20
    diffusion_class_A_max(8) = 24
    diffusion_class_A_max(9) = 10000

[simonguichandut]

Thanks @evbauer, this does work! I guess it's better to carefully setup the classes, rather than add some overhead code to avoid dumping mass into species with initial xa=1e-99. So I think this can be closed?

[evbauer]

Great! Thanks for confirming that this is an adequate solution.