Plasma treatment enhancements

[Gaetanosaure]

****************************************** EDIT 16/07/2026 *************************************************************

This PR has been split into several smaller PRs:

#2144 : reviewed and merged
#2145 : under review

---

Here are the changes proposed in this pull request. They all focus on the plasma modelling in cantera.

1. Solved the EEDF computation issue: the solver now converges to the right EEDF. As you may have noticed, the example on the website gets the right trend but the values are off (figures 1a and 1c). This is because the effective cross-sections and elastic cross-sections were treated in the same way whereas they should not be: effective = elastic + inelastic. This ended up in counting the inelastic cross-sections twice, thus lowering the EEDF value on the high energies with respect to the correct solution. As can be seen on figure 1b and 1d, this is now fixed.

Figure 1 – Example of EEDF calculation on the cantera website. a) current implementation - logarithmic x axis, b) Proposed pull request - logarithmic x axis, c) a) current implementation - linear x axis, d) Proposed pull request – linear x axis.

2. Corrected the physics of plasma heating. In the current version of the code, intrinsincHeating() returns the Joule effect plus the power of the discharge going into elastic collisions. This is incorrect since the work of the electromagnetic force is the Joule effect. In the proposed pull request, intrinsincHeating() only returns the Joule power.

3. Entering grid parameters instead of having to enter the grid itself in the YAML: this is much more compact and practical than having to copy-paste the grid by hand, especially when the grid reaches several hundreds of points (see figure 2 for the parameters).

4. Automatic grid adaptation to the problem faced. This feature allows for the adaptation of the maximum value of the energy grid so that in cases where the user is not sure up to which energies the EEDF should be resolved the code solves it for them (with a slight computational overhead). Both this feature and the previous one are controlled by a block electron-energy-distribution added in the YAML phase node looking roughly like this (depending on the user’s choices):

Figure 2 - New YAML phase node architecture for the EEDF calculation parameters.

The parameters to be chosen by the user are the following:
-> Initial max_energy_level: in eV, the energy level at which the energy grid will be initialized. This number can be led to change if the energy grid adaptation is activated.
-> Initial number of energy grid cells: the number of cells the user wishes to have in the energy grid. This number is not modified by grid adaptation in its current implementation but the name is chosen as is to keep compatibility with a potential later update.
-> Energy_levels_distribution: governs how the energy levels are spaced between themselves. Three options are currently available: Linear, Quadratic and Geometric. For the Geometric case, the user can also choose the geometric factor using the field geometric_grid_ratio.
-> energy_grid_adaptation block: used if the user wishes the grid to be dynamically adapted during the computation. If "enabled" is set to true, the solver checks the amount of decades lost between the head and the tail of the EEDF and compares it to "max_decay_decades" and "min_decay_decades" to decide whether to extend or shrink the current grid by a factor "update_factor" at the next iteration. The maximum amount of iterations allowed to the grid adaptation routine is chosen via "max_iterations".

5. Implemented a new format for the YAML files at the cross-roads between yours and ours, avoiding duplicates but solving the current problem of cross-section datasets inconsistency, all the while keeping readability. Hopefully you will find it satisfactory. Here is its presentation.
   The idea of this YAML format for electron collision reactions is to separate the data needed for the chemical solver from the data needed for the EEDF solver, the ultimate goal being to allow for a minimal user error while manipulating the cross-section dataset (which remains as close as possible to the original LXcat reference). However, you made a good point in the previous discussions by stressing that the reactions should not duplicate information from the cross-sections. Therefore, now all the information is stored in the cross-sections, which are themselves all stored in the electron-collisions YAML node for safety and clarity purposes. To link the reactions involving electron collisions to their data, a field “name” was added to the cross-sections, which now typically look like that:

This field is constructed from the cross-section data taken from the LXCat database as follows: databaseName_target_kind_product_threshold to create a unique identifier that is then declared in the appropriate reaction in the reaction YAML node if this reaction if wished by the user to be taken into account into the chemistry. Reactions involving electron-collisions now typically look like:

6. Added a reaction class for vibrational relaxation reactions, including detailed relaxation and relaxation modelled using the vibrational energy equation. Four possible models are currently available : one for detailed vibration named “multi-state-resolved” (taking reactions having a base rate constant in the form of the reactions given in the Zhong mechanism and a scaling using the harmonic approximation of the SSH theory), and three options implementing the relaxation rates of the mean vibrational energy for mean vibrational energy models: “constant”, “castela” and “starikovskiy”. More details about these models can be found in the documentation of the header of the new class VibrationalRelaxationRate.h. The main point is that all of these models implement the following equation for the reaction rates:

However, depending on the chosen model, all parameters are not required and some default values are applied.
7. Managed to go around the PlasmaReactor by implementing the mean vibrational energy equation model using fictitious species as energy reservoirs. This avoids having to declare a new class for each reactor in which one wishes to use a PlasmaPhase. Typically, this is done by, instead of declaring N2(v1), N2(v2), … N2(v45) , declaring one unique species N2(v), which enthalpy is shifted from that of ground N2 by a chosen energy E_res. This allows to collapse the 45 species equations for the vibrational levels of N2 into one species equation for N2(v). The reactions supposed to produced a vibrational level v_i are then declared as such in the YAML file:
N2 + Electron => Electron + alpha*N2(v) + (1-alpha)*N2,
Where alpha is defined as E_shift_v_i / E_res.
To avoid interfering with chemistry we recommend taking a high value with respect to dissociation energy for instance for E_res (in order to get small N2(v) concentration). To avoid mistakes and to enhance readability it is also recommended to take the same E_res for each vibrational species declared.

The changes proposed here are validated on figure 3 against the test case computed in the two following papers:

Harry Gillingham, Thomas Lesaffre, Nicolas Barléon,
A semi-analytical model for Plasma-Assisted Combustion: Application to NH3-H2-air mixtures,
Combustion and Flame,
Volume 283,
2026,
114528,
ISSN 0010-2180,
https://doi.org/10.1016/j.combustflame.2025.114528.

and

Xingqian Mao, Hongtao Zhong, Ning Liu, Ziyu Wang, Yiguang Ju,
Ignition enhancement and NOx formation of NH3/air mixtures by non-equilibrium plasma discharge,
Combustion and Flame,
Volume 259,
2024,
113140,
ISSN 0010-2180,
https://doi.org/10.1016/j.combustflame.2023.113140.

_Figure 3 - - Validation of the new cantera plasma implementations against the Mao test case._

AI Statement (required)

• Limited use of generative AI.
  Standard or boilerplate code snippets were generated with AI and manually reviewed;
  all design, logic, and implementation decisions were made by the contributor.
  Examples: IDE code-completions or brief LLM queries for common patterns.

Closes #

[speth]

Hi @Gaetanosaure, thank you for sharing this PR. It looks like there are several interesting things here. Are there some changed input files that need to be committed and included with this PR? You mention a number of changes to the input format, but I don't see any changes to any input files demonstrating them. I assume this is also the reason for some of the CI failures.

[Gaetanosaure]

Hi @speth , indeed this is a very good point, sorry! I made some tests in local but forgot to upload the changed files into the pull request. I will change the formats of every file that needs to be modified and also change the lxcat2yaml routine and let you know once it is done.

[speth]

Hi @Gaetanosaure, the process for updating the example data is a little bit tricky -- for the first step, you should leave the data/example_data submodule at the current commit (i.e., do an interactive rebase and remove commit 228d673), but make a companion pull request to that submodule's repo, https://github.com/cantera/cantera-example-data. There are detailed instructions in that repo's README file. If you have any questions or problems, please let me know.

[Gaetanosaure]

@speth if I did everything well you should see a pull request with the new YAML files on data/example_data.

Supposedly i think what i did should have removed the commit 228d673, but let me know if that is not the case.

All tests now pass with my scons test, please note that since i retained compatibility with the currently implemented format to avoid interfering with current users some tests print in the console a deprecation warning because they run on the old format.

I still have to modify the lxcat2yaml script and then I believe the changes should be comprehensive. Let me know if you see any remaining problems.

[Gaetanosaure]

@speth I have started working on the lxcat2yaml tool and just realised that in our case it is a bit risky to generate the reactions directly since the products of the reactions are very much subject to modelling choices of the user (ex: N2(A, v1-4) being treated as N2(A), N2(w) being treated as N2(B), ... therefore what i propose to do is to generate the complete electron-collisions dataset through this tool, but not to generate directly into the yaml file the corresponding reactions. Rather, i could generate a block of commented reactions in the YAML file where all reactions' product is directly the one from LXCat and let the user change this name and uncomment the reaction at will.
Another possibility for this would be to provide the lxcat2yaml function a list of correspondences wished by the user, telling the code exactly to which species of the phase they wish LXCat species names to be associated, the routine returning an error in the case of a missing match.

What do you think?

[Gaetanosaure]

Hi @speth, I have now added the lxcat2yaml converter and solved the failing tests.

Please tell me if there is still something I should do to update the PR or whether what remains to be done is on your side.

[BangShiuh]

@Gaetanosaure The changes of lxcat2yaml.py is not in the PR description. I think changing the converter can be another PR. There are backward compatibility issues that we need to carefully address. We could use two different converter for now so we can ensure the previous tests still serve its purpose.

[Gaetanosaure]

@BangShiuh @speth yes this can be considered!
However please take into account that, if I am not mistaken, the former format currently cannot distinguish between effective and elastic collisions, and as such cannot provide the correct EEDF. For now, the solver can still load the former format for compatibility but these formats must be used with care.
It is however completely possible to add a key « —deprecated » to lxcat2yaml so that it does the conversion to the former format if you believe it to be necessary, but my opinion is that keeping the compatibility to read new mechanisms is sufficient and that new conversions should be done exclusively towards the latest format for both correctness reasons but also to reduce potential users’ confusion.

Please let me know what you both think!

[BangShiuh]

Thank you for working on this. I have briefly reviewed the code. There are some issues and questions need to be addressed. Some designs such as the collapsing of inelastic requires description in the PR details.

[BangShiuh]

We fix the first point cross section data in lxcat2yaml.py

    if energy_levels[0] > threshold:
        # Use Flowlist again to ensure correct YAML format
        energy_levels = Flowlist([threshold, *energy_levels])
        cross_sections = Flowlist([0.0, *cross_sections])
    else:
        cross_sections[0] = 0.0

The threshold energy is added as the first point. In case the first point is less than or equal to the threshold, we set the cross section to zero. So if you use threshold instead, the number can be different. lxcat data do have some issues and that is one reason we use lxcat2yaml to fix those issues.

[NicolasBarleon]

In case energy_levels[0] <= threshold, then m_energyLevels[0] may not necessarily corresponds to the threshold energy. There are some strange examples such as in the Phelps database for N2

EXCITATION
N2 -> N2(v5)
1.470000e+0
SPECIES: e / N2
PROCESS: E + N2 -> E + N2(v5), Excitation
PARAM.: E = 1.47 eV, complete set
COMMENT: N2 VIBRATIONAL (V=5) SCHULZ 64.
UPDATED: 2010-02-25 15:17:35
COLUMNS: Energy (eV) | Cross section (m2)

0.000000e+0 0.000000e+0
1.470000e+0 0.000000e+0
2.100000e+0 0.000000e+0
2.200000e+0 8.250000e-21
2.300000e+0 1.230000e-20
2.400000e+0 1.530000e-20
2.500000e+0 1.440000e-20
2.600000e+0 3.450000e-21
2.700000e+0 2.250000e-22
2.750000e+0 3.450000e-21
2.800000e+0 5.400000e-21
2.900000e+0 6.600000e-21
3.000000e+0 2.175000e-21
3.100000e+0 1.050000e-21
3.200000e+0 3.150000e-21
3.300000e+0 1.035000e-21
3.400000e+0 0.000000e+0
1.000000e+3 0.000000e+0

[NicolasBarleon]

It is probably safer to keep m_threshold there.

[BangShiuh]

We can remove the excess trivial point. The energy of the first point is supposed to match the threshold value. I think we will need to choose either using the threshold energy or using the first point energy. If we choose to use the threshold energy, we can remove the first point if the energy is smaller than the threshold energy.

[NicolasBarleon]

In most cases, I would say yes.

However, this does not seem to be generic. For instance, you may refer to most of the cross-section databases for N2 vibrational excitation. Using the energy threshold is straightforward to implement and helps avoid any mistakes. For this reason, I would recommand to keep this option.

[BangShiuh]

same issue as the definition of the first point and threshold.

[NicolasBarleon]

Same comment as before, m_energyLevels[0] might not work every times.

[BangShiuh]

similar issue here for threshold energy and the first cross section point.

[BangShiuh]

No reason to touch this. This is only different in style. Same thing for other places. Too many changes on the format can make the review difficult.

[Gaetanosaure]

Indeed, this bit of code is exactly back to the current implementation on the official dev github which is:

void PlasmaPhase::normalizeElectronEnergyDistribution() {
Eigen::ArrayXd eps32 = m_electronEnergyLevels.pow(3./2.);
double norm = 2./3. * numericalQuadrature(m_quadratureMethod,
m_electronEnergyDist, eps32);
if (norm < 0.0) {
throw CanteraError("PlasmaPhase::normalizeElectronEnergyDistribution",
"The norm is negative. This might be caused by bad "
"electron energy distribution");
}
m_electronEnergyDist /= norm;
}

I know I made some tests of different functions at some point to understand whether that could be a cause for EEDF computation errors but eventually put it back to normal.

Sorry for the confusion due to my handling of github, but I confirm this bit of code is as in the original, but for some "2.0" that replaced some "2." as I just saw due to me retyping the code again when switching it back to original after some tests.

Same thing might have happened in other places in the code, therefore I apologize if this lengthened your review time as this was not my intention...

[BangShiuh]

I don't understand this. Can you explain the logic of collapsing the inelastic collision? What is the purpose of having this function? If now you move the collision data to the section collision, should the attributes of the collisions move with the data? Why do we still have collision kind under a reaction?

[Gaetanosaure]

Thank you for this question which shows I should probably have clarified things more on this bit of code.

1. What is the purpose of having this function?
   The idea of this function is to try to infer automatically the process kind and to compare it to the kind declared in the yaml data to check for potential inconsistencies that could have appeared while handling the data. As having the right kind is essential to the EEDF calculation, if the inferred kind and the declared kind do not match, an error is raised so that the user gets a look at it and checks for consistency.

2. What is this collapse that is mentionned?

It is one of the two exceptions that I found necessary to the previous check. The first one is that the inferring algorithm cannot distinguish between effective and elastic. Therefore, if one of the two is inferred and one of the two is declared, then the data is accepted, using the entry declared in the YAML.
The second exception is related to this collapse. Rotational species are not always interesting to be put in the phase, therefore a useful method to add them is to declare them simply as an ground species addition.
Thus, N2 + Electron => N2(rot) + Electron can be written in the reaction node as: N2 + Electron => N2 + Electron, all the while being referenced to the cross-section related to the generation of N2(rot) to keep EEDF consistency. But in that case, the reaction will be inferred as effective by the automatic kind detector, so this specific case needed to be allowed to pass.

3. Why is there still a collision kind under a reaction?

In reality there is not but thank you for pointing out that inconsistency in my comment which I forgot to update. There is no kind under the reaction, the kind is under the electron-collisions data in the collision corresponding to the reaction. I will update this in my next commit to make things clearer to the code reader.

[BangShiuh]

@Gaetanosaure Can you provide an example how the vibrational reactions and vibrational collisions are used? Do you link vibrational collision to different reactions or just one reaction?

[Gaetanosaure]

@BangShiuh I have provided you with an example test and all the necessary data to run it in my latest e-mail which I sent on the same day I opened the pull request. If I understood well what you are looking for, you should be able to see the structure on this test case! Please let me know if you need something more.

[BangShiuh]

Thanks. I just saw the email. I will take a look!

[speth]

Thanks for all of this work, @Gaetanosaure. I think there are some very valuable contributions here. However, this is a bit of a monster PR, which would be easier to review and revise if we could split it up into a couple discrete chunks. I'd like to try the following:

1. Items 1-4, which are all clear fixes / improvements and should be relatively straightforward to review and merge.
2. The YAML format changes. While I'm inclined to agree with the proposal to keep the full set of cross sections in a single YAML section, I think this format still needs some refinements. This should be considered as an independent PR.
3. The vibrational relaxation reaction type: I haven't yet looked at the implementation, but conceptually this sounds like a good approach, and I'm glad to see that there is an option that avoids the need for another new Reactor variant. You could include this in the first PR or separately.

If you want to make this split yourself, please go ahead. Or if you'd like, I can split the parts here and push a separate branch that you can use for the additional PR(s). Let me know either way.

[Gaetanosaure]

Hi @speth,

Thank you for this proposal of work division.
It sounds good to me, however there is a slight tweak to do in your organisation because the point 1 relies on the YAML format variation to operate successfully: "elastic" and "effective" collisions need to be properly separated, and this is not possible with the current YAML format.

I would therefore suggest to do as follows:

1. Integrate 2, 3 and 4
2. Integrate 1 and YAML format variations (point 5)
3. Integrate vibrational treatments (points 6 and 7)

If you are happy with this you can proceed with the splitting as you suggested above.
Should you need any help please tell me. Sorry for responding late this time, I will check regularly for any comments.

[speth]

Hi @Gaetanosaure, I've implemented the branch split-up:

• https://github.com/speth/cantera/tree/plasma-heating-grid (010001f) covers items 2, 3, and 4, starting from the current main branch.
• https://github.com/speth/cantera/tree/vibrational-relaxation (c60d579) covers items 5 and 6, starting from the current main branch. This one is independent of the first, so we can review/revise this simultaneously.
• plasma-collision-yaml covers items 1 and 7, and is stacked on top of plasma-heating-grid (5f31dde). We'll handle this one after the first.

I would suggest we close this PR and you can open new ones based on these branches.

I think we'll still have a little bit of complexity with the example-data, since there are changes to both the phase definition in the first PR, then changes to the reaction definitions in the 3rd, so we'll probably need to merge these simultaneously to avoid having to juggle two example-data updates. For the moment, we'll just not worry about the examples failing to run.

[Gaetanosaure]

@speth That is perfect thank you for doing this split.

I will now close the pull request and open two concurrent new ones based on plasma-heating-grid and vibrational relaxation

[Gaetanosaure]

Closing this pull request, the two new ones are coming soon.