Choosing the value of initial_f_ion (i.e. the ionized hydrogen fraction at r=1Rp)

[dlinssen]

Hi Leonardo and community,

I have noticed that the density and velocity structure of the wind can change quite a bit depending on the initial_f_ion parameter of the hydrogen.ion_fraction() function, via the mean molecular weight structure (and thus mu_bar). At first I thought that since initial_f_ion only fixes the value at the lower boundary, the value wouldn't matter too much, since the rest of the radial domain would be self-consistently solved and thus be relatively insensitive to the value at the boundary. However, when I run a generic planet, the mean molecular weight structure is different throughout a very large portion of the radial domain (see attached figure). The density structure changes by a factor ~2 depending on the value of initial_f_ion. In the figure, I've also added a simulation where I run p-winds's output structure through Cloudy, and use the mean molecular weight structure reported by Cloudy to calculate the mu_bar parameter and generate a new p-winds profile based on this value (until mu_bar converges in this way). This approach should give us a completely self-consistent mu_bar and outflow structure. I'm wondering if something similar would be possible with p-winds alone, without invoking Cloudy. Can p-winds somehow evaluate if the chosen initial_f_ion is "self-consistent" with the flow structure? And if not, do you have some insight in how to choose a reasonable value for initial_f_ion beforehand? As the Cloudy simulation shows, the flow cannot always be assumed to be completely atomic at the optical radius. I hope anyone has some ideas here!

Thanks! Dion

[ladsantos]

Thanks for the issue, Dion! This is giving me some food for thought. I ran a test where I calculate the H ion fraction using different values of initial_f_ion (0 and 1), and they yield similar profiles of velocities and densities (save a factor or ~2 or less, as you pointed out). But what bugs me is that the profile of H ionization fraction is wildly different (see plot below) and the He triplet profile is also somewhat different; hence why the mean molecular weight profile is also very different from self-consistent models.

I've been meaning to implement a more self-consistent approach to calculate the outflow structure, following, e.g., the formulation described in Allan+2023. It may take me some time to figure it out, though.

In the interim, I will search for a "post-processing" solution for this that wouldn't involve completely re-writing the parker module.

[ladsantos]

FYI, I started developing the fluid dynamics formulation (based on Allan+2023) in the branch dev. Feel free to contribute, if you would like!