Plasmapy: Add average ionization functionality

[lemmatum]

We should allow for average ionization as an argument to functions like plasma_frequency, instead of just assuming integer values of ionization based on what particle is passed.

Related to #182

[lemmatum]

We should decide whether we want this by 0.1 or not.

[namurphy]

In addition to allowing for average ionization functionality, we should also have a way to represent the full charge state distribution (e.g., what fraction of iron atoms is in each of of the Fe 0+, Fe 1+, Fe 2+, ... states).

A colleague and I have done before this for a prototype non-equilibrium ionization modeling package for solar eruptions. [The applications we had in mind were low density astrophysical plasmas that rapidly heated or cooled such that ionization and recombination were not able to keep up with the fast temperature changes.] I ended up using a dictionary of dictionaries of dictionaries, which makes me cringe a little bit when I think about it, so I think we can represent this in a way that is more reasonable and intuitive and less cumbersome. We'd also need the charge state distribution to predict spectra. I've discussed this a bit with @wtbarnes who is developing fiasco. This extra issue is beyond what we should plan on including for 0.1 though!

[namurphy]

With respect to how to modify the atomic subpackage to account for this, would it be enough to allow for syntax like: "Fe-56 2.7+"?

[wtbarnes]

Ionization equilibrium is one of the few things I've actually implemented in fiasco thus far. fiasco actually provides two ways of getting at the ionization equilibria,

1. CHIANTI ships ionization fractions (as a function of temperature) as part of the database so this can be accessed as a property of an Ion. This just pulls data straight out of the database and interpolates it to a given temperature.

>>> import fiasco
>>> import astropy.units as u
>>> import numpy as np
>>> ion = fiasco.Ion('fe_15', np.logspace(5,8,100)*u.K) #note by fe_15 I mean Fe XV or Fe +14
>>> ion.ioneq

1. Calculate the ionization fraction of every ion of an element by calculating the ionization recombination rates as a function of temperature and solving the ionization balance equations,

>>> el = fiasco.Element('iron', np.logspace(5,8,100)*u.K)
>>> ioneq = el.equilibrium_ionization()
>>> ioneq[:,14] # gives ionization fraction of Fe 15 (or +14)

(@namurphy I've also found a pretty nice way of doing the non-equilibrium ionization using an implicit method for a 1D array of temperature and density using an implicit method though this has not been merged into fiasco.)

[lemmatum]

I'm tempted to add simple average ionization handling for v0.1. It should be relatively easy, just adding a default argument Z=np.nan with a snippet like this for defaulting to particle charge if no average ionization is provided by the user.

if np.isnan(Z):
    try:
        Z = charge_state(ion)
    except ValueError:
        Z = 1
else:
    # some test to check if given ionization is above max number of electrons for the particle.

This code can then be made into its own function and invoked in each function where average ionization is relevant.

For the more complex charge state distribution idea (beyond 0.1)... one way to handle charge state distributions would be to create a chargeDistribution class which is associated with a particle and carries around information on all the relevant charge states of that particle, along with the population densities/fractions of those charge states. We could then use something like single-dispatch decorators to add functionality for handling this new class to already existing functions. The charge state densities/fractions can be either input manually or generated with a method, like the one in fiasco.

[StanczakDominik]

Just popping in to say that I like this idea. It would clear up a bunch of code in #191 and likely elsewhere.

[lemmatum]

@StanczakDominik on that note, we should consider doing the same for ion masse(s), especially with the whole reduced mass business in transport.py. I think the particles framework is starting to get fragile and in the way...

[StanczakDominik]

Yeah, perhaps it's wiser to rework that a bit before moving on with further changes in code. @namurphy, whatcha think about this?

[samurai688]

Ooh, this is the sort of mild issue that I have strong opinions about, like measuring Te in electronvolts. Which is beautiful and wonderful and gives deep connection to the atomic physics, as we recall. I apologize in advance, but here comes a mini-rant about Zbar!

As @lemmatum points out, plasma physicists have a deep-seated, primal urge to introduce average ionization values for their equations. You will see these all over the place, usually represented by the letter Z with a bar over it, or in computerese, Zbar. Also swimming around out there are with a bar over it, Z_eff / Z_effective, and so forth. I am CONSIDERABLY less happy about this than the electronvolts business.

The most common one of these I see is the arithmetic mean ionization state, which I'll refer to as Zbar. Together with the assumption of quasineutrality, this allows you to write ne = Zbar ni. This is powerful in a lot of contexts. Particularly if you are making a fluid model of a plasma, introducing Zbar in this way allows to represent a lot of the plasma behavior and you can still think of things as more or less a single fluid, and use all your traditional fluid moves.

What don't I like about this? For one thing, it evokes the idea of ions with non-integer charges, which is unphysical. You will never find a 0.5-ionized ion, nor a 0.3-ionized one, etc., unless we want to start getting philosophical about things like excited/Rydberg states. In atomic physics, the physics is set from accounting of all the different energy states, no way around it. This is also true in plasma physics... but you can kinda sorta get away with it a lot of times.

Let's say you see a nice physical equation in the formulary that's got a plain Z in it, but you've had a bad upbringing and you're thinking about a plasma as a thing that has a Zbar. Why not just whack the Zbar in there? Often, if the rigorous result contains a summation that's linear in Z, you'll get the answer exactly right. This is true for the Debye length I think.

Now let's say you want to figure out some Bremsstrahlung radiation. Uh oh, that one's got a Z^2. Lulled into a false sense of security, maybe you use the Zbar again. To make matters worse, people have defined a special Zbar (what I would call a Z_effective) that gets the answer right for this one, and will go around writing equations in textbooks with their new Z-thing. You might find such an equation and not notice the difference. You usually still won't have errors that are that bad. Or at least, you'll be able to sweep them under the rug by blaming it on something else, such as the experimentalists.

Okay, I think we see where I'm going with this. /rant. In summary, I think Zbar is a powerful approximation for plasma physics and useful in a lot of contexts. But with great power, comes great responsibility. I think a beautiful world would be one in which we all use programs like future-plasmapy, which would contain elegant structures to take care of all the accounting for us, and think about Zbars as what they are, descriptive statistics of a more complex plasma beast.

For Plasmapy, my opinion is that indeed we absolutely should have Zbar or similar Z_effective allowed in the plasma functions where reasonable, keep them very well-documented, and have consistent notation. I also think we should we absolutely NOT allow assorted Zbars to leak into the atomic physics part of the code.

[lemmatum]

Yes, I mentioned this elsewhere, though I seem to have failed to mention it here. I want zbar or zeffective for parameters that are macroscopic and to keep integer Z for microscopic (particle level stuff).

We can include machinery on the macro level for defaulting to integer Z based on particle if Zbar is not provided. We can also add functionality for calculating Zbar and Zeffective for a given charge state distribution.

I agree that these should be documented to be as precise as possible.

[StanczakDominik]

That was a really good writeup @samurai688, I'll be way more mindful of these in the future. Thanks!

I like the idea to use integer Z for actual atomic physics functionality. I'm not yet sure how the part of the code that interprets Zbars in terms of int Z should work just yet, though.

[lemmatum]

We should decide on this soon as it is holding up #191 .
I take it everyone is fine with using Zbar for macroscopic and Z for microscopic? (Making sure to distinguish these names and clearly document why we have this difference.)

[samurai688]

Yep, sounds great! Yep, and even for our Zbars I think we should strive to say what kind of Zbar we mean, whether we want the user to put in Zbar = ne / ni, or some other one, for example I just encountered Zeff = sum_i(n_i Z_i^2) / n_e. And I didn't document it, oops. Maybe I should practice what I preach :)

[namurphy]

I was also starting to get an unsettling feeling about what happens when average ionization functionality runs into Z² stuff, but I hadn't thought about it deeply. Then again, I'm an astrophysicist where 10 is the same order of magnitude as 5, 5 is the same order of magnitude as 3, and 3 is the same order of magnitude as 1, so 10 is the same order of magnitude as 1... 😹 In any case, @samurai688's comment above is important.

We should definitely write this up in the developer documentation, and probably also a high level summary of it in the narrative documentation for the atomic subpackage.

It will be important to think about how we name the variables since we want to make the code as readable as possible. For example, if we have Zbar and Zeff meaning the mean and rms charge, then we won't be able to distinguish between them when reading the code. If we did Zeff or Z_eff for the mean and Zrms or Z_rms for the root mean square then that would be a bit more readable. Instead of Z_eff, another possibility would be Z_mean, but Z_eff is more widely used. I'm less inclined to use Zbar because the bar in it directly refers to the notation rather than directly to the physical meaning.

Do we think we should allow notation like 'Fe 5.4+' and have that refer to the mean charge?

[StanczakDominik]

I'm less inclined to use Zbar because the bar in it directly refers to the notation rather than directly to the physical meaning.

Good eye for detail. Agreed here.

So, Z_eff and Z_rms? Those seem like the best options available.

[lemmatum]

Despite z_effective being more widely used, I am leaning towards z_mean. It does what is says on the tin - it's an arithmetic mean for ionization.

I don't like the idea of Fe 5.4+ notation; to me it looks like it is mixing atomic data and fluid data, especially since our notation is in atomic terms... atomic, particle=, etc. If it were in macroscopic/fluid terms then I could get behind it.

[namurphy]

I also would personally much prefer Z_mean since it says what it is and avoids ambiguity.

My vote would also be against allowing Fe 5.4+ notation, for the reasons @lemmatum mentioned and because of the chance of ambiguity. [I think one of the reasons I brought this up in the first place is because I'm starting to enjoy regular expressions and this sounded like fun!]

[lemmatum]

So separate required parameter Z_mean where it makes sense or shall we make it optional and default to Z_mean = Z_atomic with a physicswarning?

[samurai688]

I agree with what you guys have said, I like the names Z_mean and Z_rms as they are descriptive of what they are, and I like Z_effective or Z_eff for other cases. I will join in the hate on the Fe 5.4+ notation.

That's an interesting question @lemmatum, hmm... well, so we've got the particle proto-class that is currently represented by the strings. In theory, I think all the pure-Z functionality could go through that... So I think I like your idea, make any _mean or _eff parameters optional, use them if they're input, and use the particle if they're not input. Maybe without warning? Muahahaha.

You are right though @namurphy, I suppose this totally doesn't matter for astrophysics. Even in laboratory plasmas, I think one can usually get away with using the wrong Z_effective, or at least have enough other weird stuff going on that it won't be the leading-order wrong thing. :)

[namurphy]

You are right though @namurphy, I suppose this totally doesn't matter for astrophysics. Even in laboratory plasmas, I think one can usually get away with using the wrong Z_effective, or at least have enough other weird stuff going on that it won't be the leading-order wrong thing. :)

True enough, but it is still good that we're trying to do this as correctly as possible. There will undoubtedly be some applications where the difference between Z_mean and Z_rms will be important.

The comment about "leading-order wrong thing" just reminded me of my first draft of Coulomb_logarithm...

def Coulomb_logarithm():
    return 10

...which is perhaps the most simple, elegant, and readable code that I have ever written!

[namurphy]

Another thing that we should keep in mind more generally is that a lot of atomic data is accurate to ~10% to ~20% (e.g., ionization and recombination rates).

[lemmatum]

I'm pretty sure one of the lab plasma experiments I am working on analyzing/publishing uses z_rms for inferring hot electron production (temperature and total energy) based on measured hot x-rays, and could cause significant errors if the wrong Z is used.

[lemmatum]

Does #201 close this or do we need to do further work?

[StanczakDominik]

I keep forgetting to close issues after merging.