atomic number vs ion density

[ddudt]

Presently, we do not allow the user to specify the ion density profile directly. Instead, we input the atomic_number $Z_{eff}$ and then compute the ion density as:

$n_i = \frac{n_e}{Z_{eff}}$

However, the effective atomic number is commonly computed as:

$Z_{eff} = \frac{\Sigma n_i Z_i^2}{n_e}$

This will result in an incorrect ion density profile. What we really mean is the "average" atomic number defined as:

$Z_{avg} = \frac{\Sigma n_i Z_i}{n_e}$

I see two solutions to resolve this issue:

1. Update the documentation for the existing atomic_number profile and "Zeff" compute quantity to clarify that it is really the "average" atomic number.
2. Add a new equilibrium profile for ion_density. This would be an optional profile that can be used instead of atomic_number, similar to how we handle the option for either a current or iota profile.

I suggest that we implement both solutions.

[ddudt]

We require the atomic_number profile as an input to enforce that the plasma is quasi-neutral. So we do not want to add an ion_density profile input option.

[ddudt]

References indicating that our definition of $Z_{eff}$ is incorrect:

• Fusion Wiki
• TRANSP docs
• Wesson's Tokamaks textbook, equation 2.16.6

[ddudt]

The only place in the existing code where $Z_{eff}$ is used is in the bootstrap current computations. SIMSOPT also assumes that $n_i = n_e / Z_{eff}$ and notes that $Z_{eff}$ is "the average impurity charge", consistent with our usage. But Redl et al. does not define $Z_{eff}$ and Landreman et al. assumes $Z_{eff}=1$ in their results.

If the definition of $Z_{eff}$ in the bootstrap current equations is the "average" atomic number, then this issue will be resolved by clarifying our documentation. But if the Redl formula uses the definition I cited above, then we might need to input the atomic number and ion density profiles independently.

@landreman can you comment on this? Or @f0uriest can you discuss it with him?

@dpanici can you look into how V3FIT uses their profile inputs?

[f0uriest]

OK, talked to matt and looked at some papers and I'm pretty sure the Redl formula uses

$$Z_{eff} = 1/n_e \sum_i n_i Z_i^2$$

Consider a plasma with 3 species: electrons, main ions (H/D/T) and some impurity species.
We have quasineutrality:

$$n_e = n_D + Z n_Z$$

and the total pressure:

$$p = T_e n_e + T_i (n_D + n_Z)$$

In the code, we have a qty $n_i = n_e/Z_{eff}$, but we only ever really use it to calculate the pressure, so what we really want is something like $n_i = n_D + n_Z$ as an "effective" ion density. In most practical cases, $n_Z << n_D$, and $Z_{eff} \sim 1-2$ so for calculating pressure we could probably always just take $n_i = n_D = n_e$?

[ddudt]

Thanks for looking into this @f0uriest.

I think this means that we need to add a new profile input for the "effective ion density" $n_i = n_D + n_Z$, since we generally can't compute that from the $n_e$ and $Z_{eff}$ profiles. But to make the change less painful and backwards compatible for Hydrogen plasmas, we could have the following default values:

• If the atomic_number profile is not given, default to $Z_{eff}(\rho) = 1$
• If the ion_density profile is not given and $Z_{eff}(\rho) = 1$, then $n_i = n_e$. If $Z_{eff}(\rho) \neq 1$, we have to throw an error.

[f0uriest]

I'm still worried that if you give both $n_e$ and $n_i$ then there's no real way to enforce quasineutrality during optimization. In pretty much any fusion plasma, the impurity fraction is very small, maybe a few percent at most. So I think it would usually be pretty safe to always take $n_i = n_e$ for the purposes of calculating pressure, at least at the level of detail that we need for things currently in the code. The error from this is likely less than any error you get by lumping all impurities together rather than treating them separately.

[ddudt]

As it stands, the code is giving incorrect results for the self-consistent bootstrap current. It might be a small error in most cases, but it is still a bug that needs to be fixed so that this functionality works for $Z_{eff}>1$. The only way I see to resolve that is to give the $n_e$ and $n_i$ profiles independently.

Keep in mind that the kinetic profiles are only ever used for the bootstrap current calculations (at least in the present code). If a user is specifying inputs for these profiles, they are likely physically informed from experimental data or the output of a profile prediction code. We can still throw a warning whenever $n_e \neq n_i$ to make the user aware that quasi-neutrality may not be enforced.

More importantly, just because we can't guarantee quasineutrality doesn't mean that it isn't satisfied. In reality with $Z_{eff}>1$, $n_e \neq n_i$ yet quasineutrality holds, and we need to be able to handle those cases.