From 7890b50e6e1c4f52c3969b739c7260f68c9c877d Mon Sep 17 00:00:00 2001
From: Remi Lehe <remi.lehe@normalesup.org>
Date: Thu, 10 Oct 2024 16:50:46 -0700
Subject: [PATCH 01/25] Activate DMSC for ionization

---
 .../inputs_base_1d_picmi.py                   | 58 +++++++++++--------
 1 file changed, 33 insertions(+), 25 deletions(-)

diff --git a/Examples/Physics_applications/capacitive_discharge/inputs_base_1d_picmi.py b/Examples/Physics_applications/capacitive_discharge/inputs_base_1d_picmi.py
index 3de88f3b3cb..5fd19bc3b4b 100644
--- a/Examples/Physics_applications/capacitive_discharge/inputs_base_1d_picmi.py
+++ b/Examples/Physics_applications/capacitive_discharge/inputs_base_1d_picmi.py
@@ -268,32 +268,40 @@ def setup_run(self):
         #######################################################################
 
         cross_sec_direc = "../../../../warpx-data/MCC_cross_sections/He/"
-        electron_colls = picmi.MCCCollisions(
-            name="coll_elec",
-            species=self.electrons,
-            background_density=self.gas_density,
-            background_temperature=self.gas_temp,
-            background_mass=self.ions.mass,
-            ndt=self.mcc_subcycling_steps,
-            scattering_processes={
-                "elastic": {
-                    "cross_section": cross_sec_direc + "electron_scattering.dat"
-                },
-                "excitation1": {
-                    "cross_section": cross_sec_direc + "excitation_1.dat",
-                    "energy": 19.82,
-                },
-                "excitation2": {
-                    "cross_section": cross_sec_direc + "excitation_2.dat",
-                    "energy": 20.61,
-                },
-                "ionization": {
-                    "cross_section": cross_sec_direc + "ionization.dat",
-                    "energy": 24.55,
-                    "species": self.ions,
-                },
+        electron_scattering_processes = {
+            "elastic": {"cross_section": cross_sec_direc + "electron_scattering.dat"},
+            "excitation1": {
+                "cross_section": cross_sec_direc + "excitation_1.dat",
+                "energy": 19.82,
             },
-        )
+            "excitation2": {
+                "cross_section": cross_sec_direc + "excitation_2.dat",
+                "energy": 20.61,
+            },
+            "ionization": {
+                "cross_section": cross_sec_direc + "ionization.dat",
+                "energy": 24.55,
+                "species": self.ions,
+            },
+        }
+
+        if self.dsmc:
+            electron_colls = picmi.DSMCCollisions(
+                name="coll_elec",
+                species=[self.ions, self.neutrals],
+                ndt=5,
+                scattering_processes=electron_scattering_processes,
+            )
+        else:
+            electron_colls = picmi.MCCCollisions(
+                name="coll_elec",
+                species=self.electrons,
+                background_density=self.gas_density,
+                background_temperature=self.gas_temp,
+                background_mass=self.ions.mass,
+                ndt=self.mcc_subcycling_steps,
+                scattering_processes=electron_scattering_processes,
+            )
 
         ion_scattering_processes = {
             "elastic": {"cross_section": cross_sec_direc + "ion_scattering.dat"},

From 9b11492514db6934ea6430e1ee7cadbca95fbf48 Mon Sep 17 00:00:00 2001
From: Remi Lehe <remi.lehe@normalesup.org>
Date: Thu, 10 Oct 2024 17:02:33 -0700
Subject: [PATCH 02/25] Read name of the ionization species

---
 .../Collision/BinaryCollision/DSMC/DSMCFunc.H |  3 +++
 .../BinaryCollision/DSMC/DSMCFunc.cpp         | 19 ++++++++++++++++++-
 2 files changed, 21 insertions(+), 1 deletion(-)

diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H
index a692d2cbb9e..4c813ff1b5d 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H
@@ -205,7 +205,10 @@ public:
 
 private:
     amrex::Vector<ScatteringProcess> m_scattering_processes;
+    amrex::Vector<ScatteringProcess> m_ionization_processes;
     amrex::Gpu::DeviceVector<ScatteringProcess::Executor> m_scattering_processes_exe;
+    amrex::Gpu::DeviceVector<ScatteringProcess::Executor> m_ionization_processes_exe;
+
     bool m_isSameSpecies;
 
     Executor m_exe;
diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp
index e40a4e9822c..299b356404a 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp
@@ -52,7 +52,24 @@ DSMCFunc::DSMCFunc (
         WARPX_ALWAYS_ASSERT_WITH_MESSAGE(process.type() != ScatteringProcessType::INVALID,
                                         "Cannot add an unknown scattering process type");
 
-        m_scattering_processes.push_back(std::move(process));
+        // if the scattering process is ionization get the secondary species
+        // only one ionization process is supported, the vector
+        // m_ionization_processes is only used to make it simple to calculate
+        // the maximum collision frequency with the same function used for
+        // particle conserving processes
+        if (process.type() == ScatteringProcessType::IONIZATION) {
+            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(!ionization_flag,
+                                             "Background MCC only supports a single ionization process");
+            ionization_flag = true;
+
+            std::string secondary_species;
+            pp_collision_name.get("ionization_species", secondary_species);
+            m_species_names.push_back(secondary_species);
+
+            m_ionization_processes.push_back(std::move(process));
+        } else {
+            m_scattering_processes.push_back(std::move(process));
+        }
     }
 
     const int process_count = static_cast<int>(m_scattering_processes.size());

From 1ac2c44c9a16e5c1ad5c4697fe43679ae10912e4 Mon Sep 17 00:00:00 2001
From: Remi Lehe <remi.lehe@normalesup.org>
Date: Thu, 10 Oct 2024 17:23:36 -0700
Subject: [PATCH 03/25] Create ionization executors

---
 .../Collision/BinaryCollision/DSMC/DSMCFunc.H  |  3 +++
 .../BinaryCollision/DSMC/DSMCFunc.cpp          | 18 +++++++++++++-----
 2 files changed, 16 insertions(+), 5 deletions(-)

diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H
index 4c813ff1b5d..1e296414b95 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H
@@ -209,6 +209,9 @@ private:
     amrex::Gpu::DeviceVector<ScatteringProcess::Executor> m_scattering_processes_exe;
     amrex::Gpu::DeviceVector<ScatteringProcess::Executor> m_ionization_processes_exe;
 
+    bool ionization_flag = false;
+
+    amrex::Vector<std::string> m_species_names; // TODO: is this the right place to store this?
     bool m_isSameSpecies;
 
     Executor m_exe;
diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp
index 299b356404a..cdcbf61173b 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp
@@ -72,26 +72,34 @@ DSMCFunc::DSMCFunc (
         }
     }
 
-    const int process_count = static_cast<int>(m_scattering_processes.size());
-
     // Store ScatteringProcess::Executor(s).
 #ifdef AMREX_USE_GPU
     amrex::Gpu::HostVector<ScatteringProcess::Executor> h_scattering_processes_exe;
+    amrex::Gpu::HostVector<ScatteringProcess::Executor> h_ionization_processes_exe;
     for (auto const& p : m_scattering_processes) {
         h_scattering_processes_exe.push_back(p.executor());
     }
-    m_scattering_processes_exe.resize(process_count);
+    for (auto const& p : m_ionization_processes) {
+        h_ionization_processes_exe.push_back(p.executor());
+    }
+    m_scattering_processes_exe.resize(h_scattering_processes_exe.size());
+    m_ionization_processes_exe.resize(h_ionization_processes_exe.size());
     amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_scattering_processes_exe.begin(),
-                        h_scattering_processes_exe.end(), m_scattering_processes_exe.begin());
+                          h_scattering_processes_exe.end(), m_scattering_processes_exe.begin());
+    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_ionization_processes_exe.begin(),
+                          h_ionization_processes_exe.end(), m_ionization_processes_exe.begin());
     amrex::Gpu::streamSynchronize();
 #else
     for (auto const& p : m_scattering_processes) {
         m_scattering_processes_exe.push_back(p.executor());
     }
+    for (auto const& p : m_ionization_processes) {
+        m_ionization_processes_exe.push_back(p.executor());
+    }
 #endif
 
     // Link executor to appropriate ScatteringProcess executors
     m_exe.m_scattering_processes_data = m_scattering_processes_exe.data();
-    m_exe.m_process_count = process_count;
+    m_exe.m_process_count = static_cast<int>(m_scattering_processes_exe.size());
     m_exe.m_isSameSpecies = m_isSameSpecies;
 }

From 014fdedbb3ca2fb4795176dc101329224094b0af Mon Sep 17 00:00:00 2001
From: Remi Lehe <remi.lehe@normalesup.org>
Date: Thu, 10 Oct 2024 18:28:13 -0700
Subject: [PATCH 04/25] Move ionization into the same vector as the other
 processes

---
 .../DSMC/CollisionFilterFunc.H                |  1 +
 .../Collision/BinaryCollision/DSMC/DSMCFunc.H |  2 --
 .../BinaryCollision/DSMC/DSMCFunc.cpp         | 20 ++-----------------
 3 files changed, 3 insertions(+), 20 deletions(-)

diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H
index 46b228b049e..1283a086012 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H
@@ -65,6 +65,7 @@ void CollisionPairFilter (const amrex::ParticleReal u1x, const amrex::ParticleRe
 
     // Evaluate the cross-section for each scattering process to determine
     // the total collision probability.
+    // TODO: Why only 4?
     AMREX_ASSERT_WITH_MESSAGE(
         (process_count < 4), "Too many scattering processes in DSMC routine."
     );
diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H
index 1e296414b95..191cf698df1 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H
@@ -205,9 +205,7 @@ public:
 
 private:
     amrex::Vector<ScatteringProcess> m_scattering_processes;
-    amrex::Vector<ScatteringProcess> m_ionization_processes;
     amrex::Gpu::DeviceVector<ScatteringProcess::Executor> m_scattering_processes_exe;
-    amrex::Gpu::DeviceVector<ScatteringProcess::Executor> m_ionization_processes_exe;
 
     bool ionization_flag = false;
 
diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp
index cdcbf61173b..e4d17a97e19 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp
@@ -53,10 +53,7 @@ DSMCFunc::DSMCFunc (
                                         "Cannot add an unknown scattering process type");
 
         // if the scattering process is ionization get the secondary species
-        // only one ionization process is supported, the vector
-        // m_ionization_processes is only used to make it simple to calculate
-        // the maximum collision frequency with the same function used for
-        // particle conserving processes
+        // only one ionization process is supported
         if (process.type() == ScatteringProcessType::IONIZATION) {
             WARPX_ALWAYS_ASSERT_WITH_MESSAGE(!ionization_flag,
                                              "Background MCC only supports a single ionization process");
@@ -65,37 +62,24 @@ DSMCFunc::DSMCFunc (
             std::string secondary_species;
             pp_collision_name.get("ionization_species", secondary_species);
             m_species_names.push_back(secondary_species);
-
-            m_ionization_processes.push_back(std::move(process));
-        } else {
-            m_scattering_processes.push_back(std::move(process));
         }
+        m_scattering_processes.push_back(std::move(process));
     }
 
     // Store ScatteringProcess::Executor(s).
 #ifdef AMREX_USE_GPU
     amrex::Gpu::HostVector<ScatteringProcess::Executor> h_scattering_processes_exe;
-    amrex::Gpu::HostVector<ScatteringProcess::Executor> h_ionization_processes_exe;
     for (auto const& p : m_scattering_processes) {
         h_scattering_processes_exe.push_back(p.executor());
     }
-    for (auto const& p : m_ionization_processes) {
-        h_ionization_processes_exe.push_back(p.executor());
-    }
     m_scattering_processes_exe.resize(h_scattering_processes_exe.size());
-    m_ionization_processes_exe.resize(h_ionization_processes_exe.size());
     amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_scattering_processes_exe.begin(),
                           h_scattering_processes_exe.end(), m_scattering_processes_exe.begin());
-    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_ionization_processes_exe.begin(),
-                          h_ionization_processes_exe.end(), m_ionization_processes_exe.begin());
     amrex::Gpu::streamSynchronize();
 #else
     for (auto const& p : m_scattering_processes) {
         m_scattering_processes_exe.push_back(p.executor());
     }
-    for (auto const& p : m_ionization_processes) {
-        m_ionization_processes_exe.push_back(p.executor());
-    }
 #endif
 
     // Link executor to appropriate ScatteringProcess executors

From 202fab0c2b58ea7e22482917ce6f777a75022ca7 Mon Sep 17 00:00:00 2001
From: Remi Lehe <remi.lehe@normalesup.org>
Date: Thu, 10 Oct 2024 18:44:27 -0700
Subject: [PATCH 05/25] Update number of collisions

---
 .../BinaryCollision/DSMC/CollisionFilterFunc.H        | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H
index 1283a086012..d05cc4f8c4f 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H
@@ -65,12 +65,15 @@ void CollisionPairFilter (const amrex::ParticleReal u1x, const amrex::ParticleRe
 
     // Evaluate the cross-section for each scattering process to determine
     // the total collision probability.
-    // TODO: Why only 4?
+    // TODO: Why only 4 previously?
+    // Note: this is brittle because it is not AMREX_ALWAYS_ASSERT_WITH_MESSAGE
+    // In practice, a user can specify e.g. 10 scattering processes and not get an error.
+    // We should check this in the constructor of DSMCFunc ahead of time.
     AMREX_ASSERT_WITH_MESSAGE(
-        (process_count < 4), "Too many scattering processes in DSMC routine."
+        (process_count < 5), "Too many scattering processes in DSMC routine."
     );
-    int coll_type[4] = {0, 0, 0, 0};
-    amrex::ParticleReal sigma_sums[4] = {0._prt, 0._prt, 0._prt, 0._prt};
+    int coll_type[5] = {0, 0, 0, 0, 0};
+    amrex::ParticleReal sigma_sums[5] = {0._prt, 0._prt, 0._prt, 0._prt, 0._prt};
     for (int ii = 0; ii < process_count; ii++) {
         auto const& scattering_process = scattering_processes[ii];
         coll_type[ii] = int(scattering_process.m_type);

From f6a3d809d8dcc66469db214d9d4a03e2f3a735b7 Mon Sep 17 00:00:00 2001
From: Remi Lehe <remi.lehe@normalesup.org>
Date: Fri, 11 Oct 2024 10:15:54 -0700
Subject: [PATCH 06/25] Fix test: use electrons

---
 .../capacitive_discharge/inputs_base_1d_picmi.py               | 2 +-
 Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp   | 3 +++
 2 files changed, 4 insertions(+), 1 deletion(-)

diff --git a/Examples/Physics_applications/capacitive_discharge/inputs_base_1d_picmi.py b/Examples/Physics_applications/capacitive_discharge/inputs_base_1d_picmi.py
index 5fd19bc3b4b..d5e0db7213c 100644
--- a/Examples/Physics_applications/capacitive_discharge/inputs_base_1d_picmi.py
+++ b/Examples/Physics_applications/capacitive_discharge/inputs_base_1d_picmi.py
@@ -288,7 +288,7 @@ def setup_run(self):
         if self.dsmc:
             electron_colls = picmi.DSMCCollisions(
                 name="coll_elec",
-                species=[self.ions, self.neutrals],
+                species=[self.electrons, self.neutrals],
                 ndt=5,
                 scattering_processes=electron_scattering_processes,
             )
diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp
index e4d17a97e19..3fa1c9a74b1 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp
@@ -59,6 +59,9 @@ DSMCFunc::DSMCFunc (
                                              "Background MCC only supports a single ionization process");
             ionization_flag = true;
 
+            // TODO: Add a check that the first species is the electron species
+            // (This should be done for impact ionization with MCC too)
+
             std::string secondary_species;
             pp_collision_name.get("ionization_species", secondary_species);
             m_species_names.push_back(secondary_species);

From b291097fdc2ee982e5b9349c2fead5f99cb61b75 Mon Sep 17 00:00:00 2001
From: Remi Lehe <remi.lehe@normalesup.org>
Date: Fri, 11 Oct 2024 10:49:27 -0700
Subject: [PATCH 07/25] Attempt to generalize DSMC

---
 .../DSMC/SplitAndScatterFunc.cpp              | 34 +++++++++++++++----
 1 file changed, 27 insertions(+), 7 deletions(-)

diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp
index de8de9b505d..d145475119e 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp
@@ -15,17 +15,37 @@ SplitAndScatterFunc::SplitAndScatterFunc (const std::string& collision_name,
 {
     const amrex::ParmParse pp_collision_name(collision_name);
 
+    // Check if ionization is one of the scattering processes
+    bool ionization_flag = false;
+    amrex::Vector<std::string> scattering_process_names;
+    pp_collision_name.queryarr("scattering_processes", scattering_process_names);
+    for (const auto& scattering_process : scattering_process_names) {
+        if (scattering_process.find("excitation") != std::string::npos) {
+           ionization_flag = true;
+        }
+    }
+
     if (m_collision_type == CollisionType::DSMC)
     {
-        // here we can add logic to deal with cases where products are created,
-        // for example with impact ionization
-        m_num_product_species = 2;
-        m_num_products_host.push_back(1);
-        m_num_products_host.push_back(1);
+        if (ionization_flag) {
+            // Product species include the ion
+            // TODO: as an alternative, we could use the runtime attribute `ionization_level` for this species
+            m_num_product_species = 3;
+            m_num_products_host.push_back(2); // electron species:
+            // potentially 2 products per reaction: the scattered incoming electron, and the new electron from ionization
+            m_num_products_host.push_back(1);
+            m_num_products_host.push_back(1); // corresponds to the ionized species
+        } else {
+            m_num_product_species = 2;
+            m_num_products_host.push_back(1);
+            m_num_products_host.push_back(1);
+        }
+
 #ifndef AMREX_USE_GPU
         // On CPU, the device vector can be filled immediately
-        m_num_products_device.push_back(1);
-        m_num_products_device.push_back(1);
+        for (int i = 0; i < m_num_product_species; i++) {
+            m_num_products_device.push_back(m_num_products_host[i]);
+        }
 #endif
     }
     else

From 12a5de47c454deacf56b6a4ffcde59e5821623fc Mon Sep 17 00:00:00 2001
From: Remi Lehe <remi.lehe@normalesup.org>
Date: Fri, 11 Oct 2024 11:37:09 -0700
Subject: [PATCH 08/25] Outline implementation

---
 .../Collision/BinaryCollision/DSMC/DSMCFunc.H |  2 +-
 .../DSMC/SplitAndScatterFunc.H                | 37 +++++++++++++++++--
 .../DSMC/SplitAndScatterFunc.cpp              |  3 +-
 3 files changed, 35 insertions(+), 7 deletions(-)

diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H
index 191cf698df1..d3ef57e4c92 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H
@@ -207,7 +207,7 @@ private:
     amrex::Vector<ScatteringProcess> m_scattering_processes;
     amrex::Gpu::DeviceVector<ScatteringProcess::Executor> m_scattering_processes_exe;
 
-    bool ionization_flag = false;
+    bool ionization_flag = false; // Not sure we actually need to store this
 
     amrex::Vector<std::string> m_species_names; // TODO: is this the right place to store this?
     bool m_isSameSpecies;
diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
index 473199a6b21..d7167a96c87 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
@@ -85,7 +85,17 @@ public:
         amrex::Vector<int> num_added_vec(m_num_product_species);
         for (int i = 0; i < m_num_product_species; i++)
         {
-            // How many particles of product species i are created.
+            // How many particles of product species i are potentially created.
+            // ("potentially" because the actual number of particles created depends on the
+            // type of reaction in the DSMC collision, which is stored in the mask. For
+            // instance, an IONIZATION reaction will create 2 electrons (the scattered
+            // electron, which will be split off from the incident electron, and the
+            // ionized electron), while an ELASTIC reaction will only 1 electron (the
+            // scattered electron). Therefore, here we allocate enough memory to store
+            // all potential particles that can be created (see also the constructor of
+            // SplitAndScatterFunc, where `m_num_products_host` is set), but set an
+            // invalid ID for the particles that should not actually be created, so that
+            // they are removed in the subsequent call to ReDistribute.
             const index_type num_added = total * m_num_products_host[i];
             num_added_vec[i] = static_cast<int>(num_added);
             tile_products[i]->resize(products_np[i] + num_added);
@@ -120,15 +130,14 @@ public:
 #endif
 
         const int* AMREX_RESTRICT p_num_products_device = m_num_products_device.data();
+        const bool ionization_flag = m_ionization_flag;
 
         amrex::ParallelForRNG(n_total_pairs,
         [=] AMREX_GPU_DEVICE (int i, amrex::RandomEngine const& engine) noexcept
         {
             if (mask[i])
             {
-                // for now we ignore the possibility of having actual reaction
-                // products - only duplicating (splitting) of the colliding
-                // particles is supported.
+                // TODO: by default: set all IDs to be invalid
 
                 const auto product1_index = products_np_data[0] +
                                            (p_offsets[i]*p_num_products_device[0] + 0);
@@ -146,6 +155,23 @@ public:
                 // Set the weight of the new particles to p_pair_reaction_weight[i]
                 soa_products_data[1].m_rdata[PIdx::w][product2_index] = p_pair_reaction_weight[i];
 
+                if (ionization_flag) {
+                    // The electron created from ionization is stored as the second particle in product 1
+                    // Copy its properties from species 2 (momentum, position, etc.)
+
+
+                    // The ion created from ionization is stored as the first particle in product 3
+                    // Copy its properties from species 2 (momentum, position, etc.)
+                    const auto product3_index = products_np_data[2] +
+                                            (p_offsets[i]*p_num_products_device[2] + 0);
+
+                    // Make a copy of the particle from species 3
+                    copy_species2[2](soa_products_data[2], soa_3, static_cast<int>(p_pair_indices_3[i]),
+                                    static_cast<int>(product3_index), engine);
+                    // Set the weight of the new particles to p_pair_reaction_weight[i]
+                    soa_products_data[2].m_rdata[PIdx::w][product3_index] = p_pair_reaction_weight[i];
+                }
+
                 // Set the child particle properties appropriately
                 auto& ux1 = soa_products_data[0].m_rdata[PIdx::ux][product1_index];
                 auto& uy1 = soa_products_data[0].m_rdata[PIdx::uy][product1_index];
@@ -202,6 +228,8 @@ public:
                     else {
                         amrex::Abort("Uneven mass charge-exchange not implemented yet.");
                     }
+                } else if (mask[i] == int(ScatteringProcessType::IONIZATION)) {
+                    // Same as in doBackgroundIonization
                 }
                 else {
                     amrex::Abort("Unknown scattering process.");
@@ -244,6 +272,7 @@ public:
 private:
     // How many different type of species the collision produces
     int m_num_product_species;
+    bool m_ionization_flag = false;
     // Vectors of size m_num_product_species storing how many particles of a given species are
     // produced by a collision event. These vectors are duplicated (one version for host and one
     // for device) which is necessary with GPUs but redundant on CPU.
diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp
index d145475119e..b28871bd375 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp
@@ -16,12 +16,11 @@ SplitAndScatterFunc::SplitAndScatterFunc (const std::string& collision_name,
     const amrex::ParmParse pp_collision_name(collision_name);
 
     // Check if ionization is one of the scattering processes
-    bool ionization_flag = false;
     amrex::Vector<std::string> scattering_process_names;
     pp_collision_name.queryarr("scattering_processes", scattering_process_names);
     for (const auto& scattering_process : scattering_process_names) {
         if (scattering_process.find("excitation") != std::string::npos) {
-           ionization_flag = true;
+           m_ionization_flag = true;
         }
     }
 

From b5ebe9248281a8e3584e29c8ee3f19c4b56ca6fa Mon Sep 17 00:00:00 2001
From: Remi Lehe <remi.lehe@normalesup.org>
Date: Fri, 11 Oct 2024 11:43:36 -0700
Subject: [PATCH 09/25] Continue outline

---
 .../Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H   | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
index d7167a96c87..f580df48bda 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
@@ -158,17 +158,18 @@ public:
                 if (ionization_flag) {
                     // The electron created from ionization is stored as the second particle in product 1
                     // Copy its properties from species 2 (momentum, position, etc.)
+                    // TODO
+                    // Set the weight of the new particles to p_pair_reaction_weight[i]
 
 
                     // The ion created from ionization is stored as the first particle in product 3
                     // Copy its properties from species 2 (momentum, position, etc.)
                     const auto product3_index = products_np_data[2] +
                                             (p_offsets[i]*p_num_products_device[2] + 0);
+                    // TODO
 
-                    // Make a copy of the particle from species 3
-                    copy_species2[2](soa_products_data[2], soa_3, static_cast<int>(p_pair_indices_3[i]),
-                                    static_cast<int>(product3_index), engine);
                     // Set the weight of the new particles to p_pair_reaction_weight[i]
+                    soa_products_data[0].m_rdata[PIdx::w][product1_index+1] = p_pair_reaction_weight[i];
                     soa_products_data[2].m_rdata[PIdx::w][product3_index] = p_pair_reaction_weight[i];
                 }
 

From c02497773ff9dab52e5b7c1fd71431a1616c16b5 Mon Sep 17 00:00:00 2001
From: Remi Lehe <remi.lehe@normalesup.org>
Date: Fri, 11 Oct 2024 12:29:30 -0700
Subject: [PATCH 10/25] Implement creation of products of ionization

---
 .../BinaryCollision/BinaryCollision.H         | 15 ++++++++
 .../Collision/BinaryCollision/DSMC/DSMCFunc.H |  3 --
 .../BinaryCollision/DSMC/DSMCFunc.cpp         |  7 ++--
 .../DSMC/SplitAndScatterFunc.H                | 35 +++++++++++--------
 .../DSMC/SplitAndScatterFunc.cpp              |  2 +-
 5 files changed, 40 insertions(+), 22 deletions(-)

diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H
index b64c6d4b1fa..a8cb3e4970a 100644
--- a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H
+++ b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H
@@ -106,6 +106,21 @@ public:
         // Therefore, we insert the colliding species at the beginning of `m_product_species`
         if (collision_type == CollisionType::DSMC) {
             m_product_species.insert( m_product_species.begin(), m_species_names.begin(), m_species_names.end() );
+            // Check if ionization is one of the scattering processes
+            // TODO: This is reused in several places ; we should put it in a separate function
+            amrex::Vector<std::string> scattering_process_names;
+            bool ionization_flag = false;
+            pp_collision_name.queryarr("scattering_processes", scattering_process_names);
+            for (const auto& scattering_process : scattering_process_names) {
+                if (scattering_process.find("excitation") != std::string::npos) {
+                    ionization_flag = true;
+                }
+            }
+            if (ionization_flag) {
+                std::string ionization_species;
+                pp_collision_name.get("ionization_species", ionization_species);
+                m_product_species.push_back(ionization_species);
+            }
         }
         m_have_product_species = !m_product_species.empty();
 
diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H
index d3ef57e4c92..a26bb599447 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H
@@ -207,9 +207,6 @@ private:
     amrex::Vector<ScatteringProcess> m_scattering_processes;
     amrex::Gpu::DeviceVector<ScatteringProcess::Executor> m_scattering_processes_exe;
 
-    bool ionization_flag = false; // Not sure we actually need to store this
-
-    amrex::Vector<std::string> m_species_names; // TODO: is this the right place to store this?
     bool m_isSameSpecies;
 
     Executor m_exe;
diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp
index 3fa1c9a74b1..63373cce2b9 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp
@@ -32,6 +32,7 @@ DSMCFunc::DSMCFunc (
 
     // create a vector of ScatteringProcess objects from each scattering
     // process name
+    bool ionization_flag = false;
     for (const auto& scattering_process : scattering_process_names) {
         const std::string kw_cross_section = scattering_process + "_cross_section";
         std::string cross_section_file;
@@ -61,10 +62,8 @@ DSMCFunc::DSMCFunc (
 
             // TODO: Add a check that the first species is the electron species
             // (This should be done for impact ionization with MCC too)
-
-            std::string secondary_species;
-            pp_collision_name.get("ionization_species", secondary_species);
-            m_species_names.push_back(secondary_species);
+            // And add a check that the ionization species has the same mass
+            // (and a positive charge), compared to the target species
         }
         m_scattering_processes.push_back(std::move(process));
     }
diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
index f580df48bda..1171da2f1c7 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
@@ -137,8 +137,6 @@ public:
         {
             if (mask[i])
             {
-                // TODO: by default: set all IDs to be invalid
-
                 const auto product1_index = products_np_data[0] +
                                            (p_offsets[i]*p_num_products_device[0] + 0);
                 // Make a copy of the particle from species 1
@@ -147,6 +145,7 @@ public:
                 // Set the weight of the new particles to p_pair_reaction_weight[i]
                 soa_products_data[0].m_rdata[PIdx::w][product1_index] = p_pair_reaction_weight[i];
 
+                // TODO: when the reaction is ionization, this still creates a second particle of the target...
                 const auto product2_index = products_np_data[1] +
                                            (p_offsets[i]*p_num_products_device[1] + 0);
                 // Make a copy of the particle from species 2
@@ -156,21 +155,29 @@ public:
                 soa_products_data[1].m_rdata[PIdx::w][product2_index] = p_pair_reaction_weight[i];
 
                 if (ionization_flag) {
-                    // The electron created from ionization is stored as the second particle in product 1
-                    // Copy its properties from species 2 (momentum, position, etc.)
-                    // TODO
-                    // Set the weight of the new particles to p_pair_reaction_weight[i]
-
-
-                    // The ion created from ionization is stored as the first particle in product 3
-                    // Copy its properties from species 2 (momentum, position, etc.)
                     const auto product3_index = products_np_data[2] +
                                             (p_offsets[i]*p_num_products_device[2] + 0);
-                    // TODO
+                    if (mask[i] == int(ScatteringProcessType::IONIZATION)) {
+                        // The electron created from ionization is stored as the second particle in product 1
+                        // Copy its properties from the target atoms, i.e. species 2 (momentum, position, etc.)
+                        copy_species2[0](soa_products_data[0], soa_2, static_cast<int>(p_pair_indices_2[i]),
+                                        static_cast<int>(product1_index+1), engine);
+                        // Set the weight of the new particles to p_pair_reaction_weight[i]
+                        soa_products_data[0].m_rdata[PIdx::w][product1_index+1] = p_pair_reaction_weight[i];
 
-                    // Set the weight of the new particles to p_pair_reaction_weight[i]
-                    soa_products_data[0].m_rdata[PIdx::w][product1_index+1] = p_pair_reaction_weight[i];
-                    soa_products_data[2].m_rdata[PIdx::w][product3_index] = p_pair_reaction_weight[i];
+                        // The ion created from ionization is stored as the first particle in product 3
+                        // Copy its properties from the target atoms, i.e. species 2 (momentum, position, etc.)
+                        copy_species2[2](soa_products_data[2], soa_2, static_cast<int>(p_pair_indices_2[i]),
+                                        static_cast<int>(product3_index), engine);
+                        // Set the weight of the new particles to p_pair_reaction_weight[i]
+                        soa_products_data[2].m_rdata[PIdx::w][product3_index] = p_pair_reaction_weight[i];
+                    } else {
+                        // Set IDs as invalid for slots that have been allocated for potential
+                        // ionization events, so that the corresponding particles get removed in the next
+                        // call to `ReDistribute` or `RemoveInvalidParticles`
+                        soa_products_data[0].idcpu(product1_index+1) = amrex::ParticleIdCpus::Invalid;
+                        soa_products_data[2].idcpu(product3_index) = amrex::ParticleIdCpus::Invalid;
+                    }
                 }
 
                 // Set the child particle properties appropriately
diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp
index b28871bd375..070aa9caf58 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp
@@ -26,7 +26,7 @@ SplitAndScatterFunc::SplitAndScatterFunc (const std::string& collision_name,
 
     if (m_collision_type == CollisionType::DSMC)
     {
-        if (ionization_flag) {
+        if (m_ionization_flag) {
             // Product species include the ion
             // TODO: as an alternative, we could use the runtime attribute `ionization_level` for this species
             m_num_product_species = 3;

From cfc9a9b6173d2a6c727e83b35985b30c43ce0018 Mon Sep 17 00:00:00 2001
From: Remi Lehe <remi.lehe@normalesup.org>
Date: Fri, 11 Oct 2024 12:48:53 -0700
Subject: [PATCH 11/25] Modify velocity of the electron

---
 .../BinaryCollision/DSMC/SplitAndScatterFunc.H    | 15 ++++++++++++++-
 1 file changed, 14 insertions(+), 1 deletion(-)

diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
index 1171da2f1c7..0e26c6da59c 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
@@ -237,7 +237,20 @@ public:
                         amrex::Abort("Uneven mass charge-exchange not implemented yet.");
                     }
                 } else if (mask[i] == int(ScatteringProcessType::IONIZATION)) {
-                    // Same as in doBackgroundIonization
+                    // calculate kinetic energy of the incident electron
+                    double E_coll;
+                    double energy_cost = 0; // TODO: update this: get from scattering process
+                    const ParticleReal u_coll2 = ux1*ux1 + uy1*uy1 + uz1*uz1;
+                    ParticleUtils::getEnergy(u_coll2, m1, E_coll);
+                    // each of the two electron (the scattered e) gets half the energy (could change this later)
+                    const auto E_out = static_cast<amrex::ParticleReal>((E_coll - energy_cost) / 2.0_prt * PhysConst::q_e);
+                    constexpr auto c2 = PhysConst::c * PhysConst::c;
+                    const auto mc2 = m1*c2;
+                    const amrex::ParticleReal up = sqrt(E_out * (E_out + 2.0_prt*mc2) / c2) / m1;
+
+                    // isotropically scatter electrons
+                    ParticleUtils::RandomizeVelocity(ux, uy, uz, up, engine);
+                    ParticleUtils::RandomizeVelocity(e_ux, e_uy, e_uz, up, engine);
                 }
                 else {
                     amrex::Abort("Unknown scattering process.");

From d29b23f5d4d456342c1ad74e076fe0f11739817f Mon Sep 17 00:00:00 2001
From: Remi Lehe <remi.lehe@normalesup.org>
Date: Fri, 11 Oct 2024 12:51:12 -0700
Subject: [PATCH 12/25] Fix compilation errors

---
 .../BinaryCollision/DSMC/SplitAndScatterFunc.H        | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
index 0e26c6da59c..d46b7b56bac 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
@@ -240,7 +240,7 @@ public:
                     // calculate kinetic energy of the incident electron
                     double E_coll;
                     double energy_cost = 0; // TODO: update this: get from scattering process
-                    const ParticleReal u_coll2 = ux1*ux1 + uy1*uy1 + uz1*uz1;
+                    const amrex::ParticleReal u_coll2 = ux1*ux1 + uy1*uy1 + uz1*uz1;
                     ParticleUtils::getEnergy(u_coll2, m1, E_coll);
                     // each of the two electron (the scattered e) gets half the energy (could change this later)
                     const auto E_out = static_cast<amrex::ParticleReal>((E_coll - energy_cost) / 2.0_prt * PhysConst::q_e);
@@ -249,8 +249,13 @@ public:
                     const amrex::ParticleReal up = sqrt(E_out * (E_out + 2.0_prt*mc2) / c2) / m1;
 
                     // isotropically scatter electrons
-                    ParticleUtils::RandomizeVelocity(ux, uy, uz, up, engine);
-                    ParticleUtils::RandomizeVelocity(e_ux, e_uy, e_uz, up, engine);
+                    // - The incident electron
+                    ParticleUtils::RandomizeVelocity(ux1, uy1, uz1, up, engine);
+                    // - The newly created electron
+                    auto& e_ux1 = soa_products_data[0].m_rdata[PIdx::ux][product1_index+1];
+                    auto& e_uy1 = soa_products_data[0].m_rdata[PIdx::uy][product1_index+1];
+                    auto& e_uz1 = soa_products_data[0].m_rdata[PIdx::uz][product1_index+1];
+                    ParticleUtils::RandomizeVelocity(e_ux1, e_uy1, e_uz1, up, engine);
                 }
                 else {
                     amrex::Abort("Unknown scattering process.");

From 3b2001a7ec5909e8604d5fa53e355b6c8dd51bdd Mon Sep 17 00:00:00 2001
From: Remi Lehe <remi.lehe@normalesup.org>
Date: Thu, 12 Dec 2024 10:33:56 -0800
Subject: [PATCH 13/25] Update
 Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp

---
 .../Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp      | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp
index 070aa9caf58..43e3351451e 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp
@@ -19,7 +19,7 @@ SplitAndScatterFunc::SplitAndScatterFunc (const std::string& collision_name,
     amrex::Vector<std::string> scattering_process_names;
     pp_collision_name.queryarr("scattering_processes", scattering_process_names);
     for (const auto& scattering_process : scattering_process_names) {
-        if (scattering_process.find("excitation") != std::string::npos) {
+        if (scattering_process.find("ionization") != std::string::npos) {
            m_ionization_flag = true;
         }
     }

From aa9ab75982715d47c2eeb50f7e0d2293254060f4 Mon Sep 17 00:00:00 2001
From: Remi Lehe <remi.lehe@normalesup.org>
Date: Mon, 16 Dec 2024 13:51:33 -0800
Subject: [PATCH 14/25] Update
 Source/Particles/Collision/BinaryCollision/BinaryCollision.H

---
 Source/Particles/Collision/BinaryCollision/BinaryCollision.H | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H
index a8cb3e4970a..67fe9443f27 100644
--- a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H
+++ b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H
@@ -112,7 +112,7 @@ public:
             bool ionization_flag = false;
             pp_collision_name.queryarr("scattering_processes", scattering_process_names);
             for (const auto& scattering_process : scattering_process_names) {
-                if (scattering_process.find("excitation") != std::string::npos) {
+                if (scattering_process.find("ionization") != std::string::npos) {
                     ionization_flag = true;
                 }
             }

From 1b1f1581e3210eefd3abe45fd5b34d05247d1a4c Mon Sep 17 00:00:00 2001
From: Olga Shapoval <oshapoval@lbl.gov>
Date: Thu, 19 Dec 2024 17:44:06 -0800
Subject: [PATCH 15/25] Added ionization_dsmc CI test

---
 Examples/Tests/CMakeLists.txt                 |   1 +
 Examples/Tests/ionization_dsmc/CMakeLists.txt |  13 +
 .../analysis_ionization_dsmc_3d.py            | 230 ++++++++++++++++++
 .../inputs_test_3d_ionization_dsmc            | 109 +++++++++
 .../test_3d_ionization_dsmc.json              |  33 +++
 5 files changed, 386 insertions(+)
 create mode 100644 Examples/Tests/ionization_dsmc/CMakeLists.txt
 create mode 100755 Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py
 create mode 100644 Examples/Tests/ionization_dsmc/inputs_test_3d_ionization_dsmc
 create mode 100644 Regression/Checksum/benchmarks_json/test_3d_ionization_dsmc.json

diff --git a/Examples/Tests/CMakeLists.txt b/Examples/Tests/CMakeLists.txt
index 6fea9368e78..8605453a621 100644
--- a/Examples/Tests/CMakeLists.txt
+++ b/Examples/Tests/CMakeLists.txt
@@ -25,6 +25,7 @@ add_subdirectory(gaussian_beam)
 add_subdirectory(implicit)
 add_subdirectory(initial_distribution)
 add_subdirectory(initial_plasma_profile)
+add_subdirectory(ionization_dsmc)
 add_subdirectory(field_ionization)
 add_subdirectory(ion_stopping)
 add_subdirectory(langmuir)
diff --git a/Examples/Tests/ionization_dsmc/CMakeLists.txt b/Examples/Tests/ionization_dsmc/CMakeLists.txt
new file mode 100644
index 00000000000..085902915ff
--- /dev/null
+++ b/Examples/Tests/ionization_dsmc/CMakeLists.txt
@@ -0,0 +1,13 @@
+# Add tests (alphabetical order) ##############################################
+#
+
+add_warpx_test(
+    test_3d_ionization_dsmc  # name
+    3  # dims
+    2  # nprocs
+    inputs_test_3d_ionization_dsmc  # inputs
+    analysis_ionization_dsmc_3d.py  # analysis
+    diags/diag1/  # output
+    OFF  # dependency
+)
+
diff --git a/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py b/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py
new file mode 100755
index 00000000000..7e3801481f1
--- /dev/null
+++ b/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py
@@ -0,0 +1,230 @@
+#!/usr/bin/env python3
+# DSMC ionization test script: 
+#   - compares WarpX simulation results with theoretical model predictions.
+import os
+import sys
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.stats import qmc
+import argparse
+from openpmd_viewer import OpenPMDTimeSeries
+import tqdm
+from pywarpx import picmi
+constants = picmi.constants
+
+sys.path.insert(1, "../../../../warpx/Regression/Checksum/")
+from checksumAPI import evaluate_checksum
+
+filename = sys.argv[1]
+testname = os.path.split(os.getcwd())[1]
+ts = OpenPMDTimeSeries(filename)
+
+q_e = constants.q_e
+m_p = constants.m_p
+m_e = constants.m_e
+k_B = constants.kb
+ep0 = constants.ep0
+clight = constants.c
+
+plasma_density = 1e14
+neutral_density = 1e20
+dt = 1e-9
+electron_temp = 10
+neutral_temp = 300
+max_steps = 250
+max_time = max_steps * dt
+
+L = [0.1] * 3
+
+sigma_iz_file = '../../../../warpx-data/MCC_cross_sections/Xe/ionization.dat'
+iz_data = np.loadtxt(sigma_iz_file)
+
+energy_eV = iz_data[:, 0]
+sigma_m2 = iz_data[:, 1]
+iz_energy = energy_eV[0]
+
+
+def get_Te(ts):
+    T_e = []
+    for iteration in tqdm.tqdm(ts.iterations):   
+        ux, uy, uz = ts.get_particle(["ux", "uy", "uz"], species="electrons", iteration=iteration)
+        v_std_x = np.std(ux * clight)
+        v_std_y = np.std(uy * clight)
+        v_std_z = np.std(uz * clight)
+        v_std = (v_std_x + v_std_y + v_std_z) / 3
+        T_e.append( m_e * v_std**2 / q_e)
+    return T_e        
+
+def get_density(ts):
+    number_data = np.genfromtxt("./diags/counts.txt")
+    Te = get_Te(ts)
+    total_volume = L[0] * L[1] * L[2]
+    electron_weight = number_data[:, 8]
+    neutral_weight = number_data[:, 9]
+    ne = electron_weight / total_volume
+    nn = neutral_weight / total_volume
+    print(np.size(ne), np.size(Te))
+
+    return [ne, nn, ne*Te]
+
+def compute_rate_coefficients(temperatures_eV, energy_eV, sigma_m2, num_samples = 1024):
+    # """Integrate cross sections over maxwellian VDF to obtain reaction rate coefficients
+    # Given electron energy in eV (`energy_eV`) and reaction cross sections at those energies (`sigma_m2`),
+    # this function computes the reaction rate coefficient $k(T_e)$ for maxwellian electrons at
+    # a provided list of electron temperatures `temperatures_eV`.
+
+    # The rate coefficient is given by
+
+    # $$
+    # k(T_e) = \int \sigma(E) E dE
+    # $$
+
+    # where the energy $E$ is drawn from a maxwellian distribution function with zero speed and temperature $T_e$.
+    # We solve this using a quasi-monte carlo approach, by drawing a large number of low-discrepancy samples from
+    # the appropriate distribution and obtaining the average of $\sigma(E) E$.
+    # """
+    thermal_speed_scale = np.sqrt(q_e / m_e)
+    k = np.zeros(temperatures_eV.size)
+    
+    # obtain low-discrepancy samples of normal dist
+    dist = qmc.MultivariateNormalQMC(np.zeros(3), np.eye(3))
+    v = dist.random(num_samples)
+
+    for (i, T) in enumerate(temperatures_eV):
+        # scale velocities to proper temperature
+        # compute energies corresponding to each sampled velocity vector
+        speed_squared = (v[:, 0]**2 + v[:, 1]**2 + v[:, 2]**2) * T
+        e = 0.5 * speed_squared
+        speed = np.sqrt(speed_squared) * thermal_speed_scale
+        # get cross section by interpolating on table
+        sigma = np.interp(e, energy_eV, sigma_m2, left = 0)
+        k[i] = np.mean(sigma * speed)
+    
+    return k
+
+def rhs(state, params):
+    # """ Compute the right-hand side of ODE system that solves the global model described below.
+    # The global model solves for the evolution of plasma density ($n_e$), neutral density ($n_n), 
+    # and electron temperature ($T_e$) in the presence of ionization.
+    # The model equations consist of a continuity equation for electrons and neutrals,
+    # combined with an energy equation for electrons.
+    
+    # $$
+    # \frac{\partial n_e}{\partial t} = \dot{n}
+    # \frac{\partial n_n}{\partial t} = -\dot{n}
+    # \frac{3}{2}\frac{\partial n_e T_e}{\partial t} = \dot{n} \epsilon_{iz},
+    # $$
+    
+    # where 
+    
+    # $$
+    # \dot{n} = n_n n_e k_{iz}(T_e),
+    # $$
+
+    # $k_iz$ is the ionization rate coefficient as a function of electron temperature in eV, and $\epsilon_{iz}$ is the ionization energy cost in eV.
+    # """
+    # unpack parameters
+    E_iz, Te_table, kiz_table = params
+    ne, nn, energy = state[0], state[1], state[2]
+
+    # compute rhs
+    Te = energy / ne / 1.5
+    rate_coeff = np.interp(Te, Te_table, kiz_table, left = 0)
+    ndot = ne * nn * rate_coeff
+    f = np.empty(3)
+    
+    # fill in ionization rate
+    f[0] =  ndot   # d(ne)/dt
+
+    # at present, the dsmc solver does not deplete electron energy but does deplete neutrals
+    # the opposite is true for the MCC solver
+    f[1] = -ndot         # d(nn)/dt
+    f[2] = 0             # d(ne*eps) / dt
+
+    return f
+    
+def solve_theory_model():
+    # integrate cross-section table to get rate coefficients
+    Te_table = np.linspace(0, 2 * electron_temp, 256)
+    kiz_table = compute_rate_coefficients(Te_table, energy_eV, sigma_m2, num_samples = 4096)
+
+    # set up system
+    num_steps = max_steps+1
+    state_vec = np.zeros((num_steps, 3))
+    state_vec[0, :] = np.array([plasma_density, neutral_density, 1.5 * plasma_density * electron_temp])
+    t = np.linspace(0, max_time, num_steps)
+    params = (iz_energy, Te_table, kiz_table)
+    
+    # solve the system (use RK2 integration)
+    for i in range(1, num_steps):
+        u = state_vec[i-1, :]
+        k1 = rhs(u, params)
+        k2 = rhs(u + k1*dt, params)
+        state_vec[i, :] = u + (k1 + k2)*dt/2
+
+    # return result
+    ne = state_vec[:, 0]
+    nn = state_vec[:, 1]
+    Te = state_vec[:, 2] / (1.5 * ne)
+    print(np.size(ne), np.size(Te))
+    return t, [ne, nn, ne*Te]
+
+
+t_warpx = np.loadtxt("diags/counts.txt")[:,1]
+data_warpx = get_density(ts)
+
+t_theory, data_theory = solve_theory_model()
+
+fig, axs = plt.subplots(1, 3, figsize=(12, 4))
+
+# Plot 1
+method = "dsmc"
+labels = ["$n_e$ [m$^{-3}$]", "$n_n$ [m${-3}$]", "$n_e T_e$ [eVm$^{-3}$]"]
+titles = ["Plasma density", "Neutral density", "Normalized electron temperature"]
+        
+for (i, (title, label, field_warpx, field_theory)) in enumerate(zip(titles, labels, data_warpx, data_theory)):
+    axs[i].set_ylabel(label)
+    axs[i].set_title(title)
+    axs[i].plot(t_warpx, field_warpx, label = "WarpX (" + method + ')' )
+    axs[i].plot(t_theory, field_theory, label = "theory", color = 'red', ls = '--' )
+
+    axs[i].legend()
+plt.tight_layout()
+plt.savefig('ionization_dsmc__density_Te.png', dpi=150)    
+
+
+def check_tolerance(array, tolerance):
+    assert np.all(
+        array <= tolerance
+    ), f"Test did not pass: one or more elements exceed the tolerance of {tolerance}."
+    print("All elements of are within the tolerance.")
+
+# For now, setting very high tolerances on purpose. This needs to be modified once all features are implemented in the PR.
+tolerances = [1, 1e-6, 1e-2]
+def check_tolerance(array, tolerance):
+    assert np.all(
+        array <= tolerance
+    ), f"Test did not pass: one or more elements exceed the tolerance of {tolerance}."
+    print("All elements of are within the tolerance.")
+
+# For now, setting very high tolerances on purpose. This needs to be modified once all features are implemented in the PR.
+tolerance = 1 #5e-4
+
+plt.figure()
+labels = ["Plasma density $(n_e$)", "Neutral density $(n_n$)", "Normalized electron temperature $(T_en_e$)"]
+for (i, (label, field_warpx, field_theory, tolerance)) in enumerate(zip( labels, data_warpx, data_theory, tolerances)):
+
+    relative_error = np.array(abs((data_warpx[i] - data_theory[i][::5])/data_theory[i][::5]))
+    plt.plot(t_warpx, relative_error, label = label)
+    plt.ylabel('Relative error')
+    plt.xlabel('Time [s]')
+    plt.legend()
+    check_tolerance(relative_error, tolerance)
+plt.savefig('./relative_error_density_Te.png', dpi=150) 
+
+# compare checksums
+evaluate_checksum(
+    test_name=os.path.split(os.getcwd())[1],
+    output_file=sys.argv[1],
+    output_format="openpmd",
+)
\ No newline at end of file
diff --git a/Examples/Tests/ionization_dsmc/inputs_test_3d_ionization_dsmc b/Examples/Tests/ionization_dsmc/inputs_test_3d_ionization_dsmc
new file mode 100644
index 00000000000..8e430635787
--- /dev/null
+++ b/Examples/Tests/ionization_dsmc/inputs_test_3d_ionization_dsmc
@@ -0,0 +1,109 @@
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+max_step = 250
+amr.n_cell = 8 8 8
+amr.max_level = 0
+
+warpx.do_electrostatic = labframe
+algo.particle_shape = 1
+amrex.verbose = 1
+geometry.coord_sys = 0
+#################################
+############ CONSTANTS #############
+#################################
+my_constants.Te = 10
+warpx.const_dt = 1e-09
+
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+geometry.dims = 3
+geometry.prob_hi = 0.1 0.1 0.1
+geometry.prob_lo = 0 0 0
+amr.max_grid_size = 8
+
+#################################
+###### BOUNDARY CONDITIONS ######
+#################################
+geometry.is_periodic = 1 1 1
+boundary.field_hi = periodic periodic periodic
+boundary.field_lo = periodic periodic periodic
+boundary.particle_hi = periodic periodic periodic
+boundary.particle_lo = periodic periodic periodic
+
+#################################
+############ PLASMA #############
+#################################
+particles.species_names = electrons ions neutrals
+
+electrons.charge = -q_e
+electrons.density = 1e+14
+electrons.initialize_self_fields = 0
+electrons.injection_style = nuniformpercell
+electrons.mass = m_e
+electrons.momentum_distribution_type = gaussian
+electrons.num_particles_per_cell_each_dim = 4 4 4
+electrons.profile = constant
+electrons.ux_m = 0.0
+electrons.ux_th =  sqrt(q_e * Te / m_e) / clight
+electrons.uy_m = 0.0
+electrons.uy_th =  sqrt(q_e * Te / m_e)/ clight
+electrons.uz_m = 0.0
+electrons.uz_th =  sqrt(q_e * Te / m_e)/ clight
+
+ions.charge = q_e
+ions.density = 1e+14
+ions.initialize_self_fields = 0
+ions.injection_style = nuniformpercell
+ions.mass = 2.196035502270312e-25
+ions.momentum_distribution_type = gaussian
+ions.num_particles_per_cell_each_dim = 4 4 4
+ions.profile = constant
+ions.ux_m = 0.0
+ions.ux_th = 4.5810168302300867e-07
+ions.uy_m = 0.0
+ions.uy_th = 4.5810168302300867e-07
+ions.uz_m = 0.0
+ions.uz_th = 4.5810168302300867e-07
+
+neutrals.charge = 0
+neutrals.density = 1e+20
+neutrals.initialize_self_fields = 0
+neutrals.injection_style = nuniformpercell
+neutrals.mass = 2.196035502270312e-25
+neutrals.momentum_distribution_type = gaussian
+neutrals.num_particles_per_cell_each_dim = 4 4 4
+neutrals.profile = constant
+neutrals.ux_m = 0.0
+neutrals.ux_th = 4.5810168302300867e-07
+neutrals.uy_m = 0.0
+neutrals.uy_th = 4.5810168302300867e-07
+neutrals.uz_m = 0.0
+neutrals.uz_th = 4.5810168302300867e-07
+
+collisions.collision_names = coll_elec
+
+coll_elec.ionization_cross_section = ../../../../warpx-data/MCC_cross_sections/Xe/ionization.dat
+coll_elec.ionization_energy = 12.1298431
+coll_elec.ionization_species = ions
+coll_elec.ndt = 1
+coll_elec.scattering_processes = ionization
+coll_elec.species = electrons neutrals
+coll_elec.type = dsmc
+
+#################################
+############ DIAGNOSTICS #############
+#################################
+diagnostics.diags_names = diag1
+warpx.reduced_diags_names = counts
+counts.intervals = 5
+counts.path = diags/
+counts.type = ParticleNumber
+
+# Diagnostics
+diag1.intervals = 5
+diag1.diag_type = Full
+diag1.electrons.variables = ux uy uz
+diag1.neutrals.variables = ux uy uz
+diag1.format = openpmd
diff --git a/Regression/Checksum/benchmarks_json/test_3d_ionization_dsmc.json b/Regression/Checksum/benchmarks_json/test_3d_ionization_dsmc.json
new file mode 100644
index 00000000000..6681b90f1c9
--- /dev/null
+++ b/Regression/Checksum/benchmarks_json/test_3d_ionization_dsmc.json
@@ -0,0 +1,33 @@
+{
+  "lev=0": {
+      "Bx": 0.0,
+      "By": 0.0,
+      "Bz": 0.0,
+      "Ex": 33418.976671097655,
+      "Ey": 28556.630099887625,
+      "Ez": 33908.990756787396,
+      "jx": 351.709319628182,
+      "jy": 389.97631892921305,
+      "jz": 366.54101928747855
+    },
+    "electrons": {
+      "particle_momentum_x": 4.0749126130424784e-20,
+      "particle_momentum_y": 4.086095034273236e-20,
+      "particle_momentum_z": 4.0905976427530184e-20
+    },
+    "ions": {
+      "particle_position_x": 2503.961424209735,
+      "particle_position_y": 2503.383309677753,
+      "particle_position_z": 2505.603131072416,
+      "particle_momentum_x": 1.2028342698649854e-18,
+      "particle_momentum_y": 1.2008761224148308e-18,
+      "particle_momentum_z": 1.2009924309769205e-18,
+      "particle_weight": 152685546875.00003
+    },
+    "neutrals": {
+      "particle_momentum_x": 1.204336106475819e-18,
+      "particle_momentum_y": 1.2126021881554375e-18,
+      "particle_momentum_z": 1.2056811544009856e-18
+    }
+  }
+  
\ No newline at end of file

From 0b7be809f6167e9d9b925f4d6579c033dd718150 Mon Sep 17 00:00:00 2001
From: Olga Shapoval <oshapoval@lbl.gov>
Date: Thu, 19 Dec 2024 17:48:42 -0800
Subject: [PATCH 16/25] Clean-up

---
 Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py b/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py
index 7e3801481f1..dc0bbee4c0f 100755
--- a/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py
+++ b/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py
@@ -207,9 +207,6 @@ def check_tolerance(array, tolerance):
     ), f"Test did not pass: one or more elements exceed the tolerance of {tolerance}."
     print("All elements of are within the tolerance.")
 
-# For now, setting very high tolerances on purpose. This needs to be modified once all features are implemented in the PR.
-tolerance = 1 #5e-4
-
 plt.figure()
 labels = ["Plasma density $(n_e$)", "Neutral density $(n_n$)", "Normalized electron temperature $(T_en_e$)"]
 for (i, (label, field_warpx, field_theory, tolerance)) in enumerate(zip( labels, data_warpx, data_theory, tolerances)):

From ac4b357345a7f9504f84c86277f4da6197f4ae9e Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Fri, 20 Dec 2024 03:45:10 +0000
Subject: [PATCH 17/25] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 Examples/Tests/ionization_dsmc/CMakeLists.txt |   1 -
 .../analysis_ionization_dsmc_3d.py            | 125 +++++++++++-------
 .../test_3d_ionization_dsmc.json              |   1 -
 3 files changed, 75 insertions(+), 52 deletions(-)

diff --git a/Examples/Tests/ionization_dsmc/CMakeLists.txt b/Examples/Tests/ionization_dsmc/CMakeLists.txt
index 085902915ff..27514622d62 100644
--- a/Examples/Tests/ionization_dsmc/CMakeLists.txt
+++ b/Examples/Tests/ionization_dsmc/CMakeLists.txt
@@ -10,4 +10,3 @@ add_warpx_test(
     diags/diag1/  # output
     OFF  # dependency
 )
-
diff --git a/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py b/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py
index dc0bbee4c0f..1e7ee12d31b 100755
--- a/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py
+++ b/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py
@@ -1,15 +1,17 @@
 #!/usr/bin/env python3
-# DSMC ionization test script: 
+# DSMC ionization test script:
 #   - compares WarpX simulation results with theoretical model predictions.
 import os
 import sys
+
 import matplotlib.pyplot as plt
 import numpy as np
-from scipy.stats import qmc
-import argparse
-from openpmd_viewer import OpenPMDTimeSeries
 import tqdm
+from openpmd_viewer import OpenPMDTimeSeries
+from scipy.stats import qmc
+
 from pywarpx import picmi
+
 constants = picmi.constants
 
 sys.path.insert(1, "../../../../warpx/Regression/Checksum/")
@@ -36,7 +38,7 @@
 
 L = [0.1] * 3
 
-sigma_iz_file = '../../../../warpx-data/MCC_cross_sections/Xe/ionization.dat'
+sigma_iz_file = "../../../../warpx-data/MCC_cross_sections/Xe/ionization.dat"
 iz_data = np.loadtxt(sigma_iz_file)
 
 energy_eV = iz_data[:, 0]
@@ -46,14 +48,17 @@
 
 def get_Te(ts):
     T_e = []
-    for iteration in tqdm.tqdm(ts.iterations):   
-        ux, uy, uz = ts.get_particle(["ux", "uy", "uz"], species="electrons", iteration=iteration)
+    for iteration in tqdm.tqdm(ts.iterations):
+        ux, uy, uz = ts.get_particle(
+            ["ux", "uy", "uz"], species="electrons", iteration=iteration
+        )
         v_std_x = np.std(ux * clight)
         v_std_y = np.std(uy * clight)
         v_std_z = np.std(uz * clight)
         v_std = (v_std_x + v_std_y + v_std_z) / 3
-        T_e.append( m_e * v_std**2 / q_e)
-    return T_e        
+        T_e.append(m_e * v_std**2 / q_e)
+    return T_e
+
 
 def get_density(ts):
     number_data = np.genfromtxt("./diags/counts.txt")
@@ -65,9 +70,10 @@ def get_density(ts):
     nn = neutral_weight / total_volume
     print(np.size(ne), np.size(Te))
 
-    return [ne, nn, ne*Te]
+    return [ne, nn, ne * Te]
+
 
-def compute_rate_coefficients(temperatures_eV, energy_eV, sigma_m2, num_samples = 1024):
+def compute_rate_coefficients(temperatures_eV, energy_eV, sigma_m2, num_samples=1024):
     # """Integrate cross sections over maxwellian VDF to obtain reaction rate coefficients
     # Given electron energy in eV (`energy_eV`) and reaction cross sections at those energies (`sigma_m2`),
     # this function computes the reaction rate coefficient $k(T_e)$ for maxwellian electrons at
@@ -85,38 +91,39 @@ def compute_rate_coefficients(temperatures_eV, energy_eV, sigma_m2, num_samples
     # """
     thermal_speed_scale = np.sqrt(q_e / m_e)
     k = np.zeros(temperatures_eV.size)
-    
+
     # obtain low-discrepancy samples of normal dist
     dist = qmc.MultivariateNormalQMC(np.zeros(3), np.eye(3))
     v = dist.random(num_samples)
 
-    for (i, T) in enumerate(temperatures_eV):
+    for i, T in enumerate(temperatures_eV):
         # scale velocities to proper temperature
         # compute energies corresponding to each sampled velocity vector
-        speed_squared = (v[:, 0]**2 + v[:, 1]**2 + v[:, 2]**2) * T
+        speed_squared = (v[:, 0] ** 2 + v[:, 1] ** 2 + v[:, 2] ** 2) * T
         e = 0.5 * speed_squared
         speed = np.sqrt(speed_squared) * thermal_speed_scale
         # get cross section by interpolating on table
-        sigma = np.interp(e, energy_eV, sigma_m2, left = 0)
+        sigma = np.interp(e, energy_eV, sigma_m2, left=0)
         k[i] = np.mean(sigma * speed)
-    
+
     return k
 
+
 def rhs(state, params):
     # """ Compute the right-hand side of ODE system that solves the global model described below.
-    # The global model solves for the evolution of plasma density ($n_e$), neutral density ($n_n), 
+    # The global model solves for the evolution of plasma density ($n_e$), neutral density ($n_n),
     # and electron temperature ($T_e$) in the presence of ionization.
     # The model equations consist of a continuity equation for electrons and neutrals,
     # combined with an energy equation for electrons.
-    
+
     # $$
     # \frac{\partial n_e}{\partial t} = \dot{n}
     # \frac{\partial n_n}{\partial t} = -\dot{n}
     # \frac{3}{2}\frac{\partial n_e T_e}{\partial t} = \dot{n} \epsilon_{iz},
     # $$
-    
-    # where 
-    
+
+    # where
+
     # $$
     # \dot{n} = n_n n_e k_{iz}(T_e),
     # $$
@@ -129,48 +136,53 @@ def rhs(state, params):
 
     # compute rhs
     Te = energy / ne / 1.5
-    rate_coeff = np.interp(Te, Te_table, kiz_table, left = 0)
+    rate_coeff = np.interp(Te, Te_table, kiz_table, left=0)
     ndot = ne * nn * rate_coeff
     f = np.empty(3)
-    
+
     # fill in ionization rate
-    f[0] =  ndot   # d(ne)/dt
+    f[0] = ndot  # d(ne)/dt
 
     # at present, the dsmc solver does not deplete electron energy but does deplete neutrals
     # the opposite is true for the MCC solver
-    f[1] = -ndot         # d(nn)/dt
-    f[2] = 0             # d(ne*eps) / dt
+    f[1] = -ndot  # d(nn)/dt
+    f[2] = 0  # d(ne*eps) / dt
 
     return f
-    
+
+
 def solve_theory_model():
     # integrate cross-section table to get rate coefficients
     Te_table = np.linspace(0, 2 * electron_temp, 256)
-    kiz_table = compute_rate_coefficients(Te_table, energy_eV, sigma_m2, num_samples = 4096)
+    kiz_table = compute_rate_coefficients(
+        Te_table, energy_eV, sigma_m2, num_samples=4096
+    )
 
     # set up system
-    num_steps = max_steps+1
+    num_steps = max_steps + 1
     state_vec = np.zeros((num_steps, 3))
-    state_vec[0, :] = np.array([plasma_density, neutral_density, 1.5 * plasma_density * electron_temp])
+    state_vec[0, :] = np.array(
+        [plasma_density, neutral_density, 1.5 * plasma_density * electron_temp]
+    )
     t = np.linspace(0, max_time, num_steps)
     params = (iz_energy, Te_table, kiz_table)
-    
+
     # solve the system (use RK2 integration)
     for i in range(1, num_steps):
-        u = state_vec[i-1, :]
+        u = state_vec[i - 1, :]
         k1 = rhs(u, params)
-        k2 = rhs(u + k1*dt, params)
-        state_vec[i, :] = u + (k1 + k2)*dt/2
+        k2 = rhs(u + k1 * dt, params)
+        state_vec[i, :] = u + (k1 + k2) * dt / 2
 
     # return result
     ne = state_vec[:, 0]
     nn = state_vec[:, 1]
     Te = state_vec[:, 2] / (1.5 * ne)
     print(np.size(ne), np.size(Te))
-    return t, [ne, nn, ne*Te]
+    return t, [ne, nn, ne * Te]
 
 
-t_warpx = np.loadtxt("diags/counts.txt")[:,1]
+t_warpx = np.loadtxt("diags/counts.txt")[:, 1]
 data_warpx = get_density(ts)
 
 t_theory, data_theory = solve_theory_model()
@@ -181,16 +193,18 @@ def solve_theory_model():
 method = "dsmc"
 labels = ["$n_e$ [m$^{-3}$]", "$n_n$ [m${-3}$]", "$n_e T_e$ [eVm$^{-3}$]"]
 titles = ["Plasma density", "Neutral density", "Normalized electron temperature"]
-        
-for (i, (title, label, field_warpx, field_theory)) in enumerate(zip(titles, labels, data_warpx, data_theory)):
+
+for i, (title, label, field_warpx, field_theory) in enumerate(
+    zip(titles, labels, data_warpx, data_theory)
+):
     axs[i].set_ylabel(label)
     axs[i].set_title(title)
-    axs[i].plot(t_warpx, field_warpx, label = "WarpX (" + method + ')' )
-    axs[i].plot(t_theory, field_theory, label = "theory", color = 'red', ls = '--' )
+    axs[i].plot(t_warpx, field_warpx, label="WarpX (" + method + ")")
+    axs[i].plot(t_theory, field_theory, label="theory", color="red", ls="--")
 
     axs[i].legend()
 plt.tight_layout()
-plt.savefig('ionization_dsmc__density_Te.png', dpi=150)    
+plt.savefig("ionization_dsmc__density_Te.png", dpi=150)
 
 
 def check_tolerance(array, tolerance):
@@ -199,29 +213,40 @@ def check_tolerance(array, tolerance):
     ), f"Test did not pass: one or more elements exceed the tolerance of {tolerance}."
     print("All elements of are within the tolerance.")
 
+
 # For now, setting very high tolerances on purpose. This needs to be modified once all features are implemented in the PR.
 tolerances = [1, 1e-6, 1e-2]
+
+
 def check_tolerance(array, tolerance):
     assert np.all(
         array <= tolerance
     ), f"Test did not pass: one or more elements exceed the tolerance of {tolerance}."
     print("All elements of are within the tolerance.")
 
-plt.figure()
-labels = ["Plasma density $(n_e$)", "Neutral density $(n_n$)", "Normalized electron temperature $(T_en_e$)"]
-for (i, (label, field_warpx, field_theory, tolerance)) in enumerate(zip( labels, data_warpx, data_theory, tolerances)):
 
-    relative_error = np.array(abs((data_warpx[i] - data_theory[i][::5])/data_theory[i][::5]))
-    plt.plot(t_warpx, relative_error, label = label)
-    plt.ylabel('Relative error')
-    plt.xlabel('Time [s]')
+plt.figure()
+labels = [
+    "Plasma density $(n_e$)",
+    "Neutral density $(n_n$)",
+    "Normalized electron temperature $(T_en_e$)",
+]
+for i, (label, field_warpx, field_theory, tolerance) in enumerate(
+    zip(labels, data_warpx, data_theory, tolerances)
+):
+    relative_error = np.array(
+        abs((data_warpx[i] - data_theory[i][::5]) / data_theory[i][::5])
+    )
+    plt.plot(t_warpx, relative_error, label=label)
+    plt.ylabel("Relative error")
+    plt.xlabel("Time [s]")
     plt.legend()
     check_tolerance(relative_error, tolerance)
-plt.savefig('./relative_error_density_Te.png', dpi=150) 
+plt.savefig("./relative_error_density_Te.png", dpi=150)
 
 # compare checksums
 evaluate_checksum(
     test_name=os.path.split(os.getcwd())[1],
     output_file=sys.argv[1],
     output_format="openpmd",
-)
\ No newline at end of file
+)
diff --git a/Regression/Checksum/benchmarks_json/test_3d_ionization_dsmc.json b/Regression/Checksum/benchmarks_json/test_3d_ionization_dsmc.json
index 6681b90f1c9..6743c19ecb8 100644
--- a/Regression/Checksum/benchmarks_json/test_3d_ionization_dsmc.json
+++ b/Regression/Checksum/benchmarks_json/test_3d_ionization_dsmc.json
@@ -30,4 +30,3 @@
       "particle_momentum_z": 1.2056811544009856e-18
     }
   }
-  
\ No newline at end of file

From 5f9aae5239f02d308ba6f729360c3a28ada68972 Mon Sep 17 00:00:00 2001
From: Olga Shapoval <oshapoval@lbl.gov>
Date: Thu, 19 Dec 2024 23:10:33 -0800
Subject: [PATCH 18/25] Updated CI: splited checksum analysis from test
 analysis

---
 Examples/Tests/ionization_dsmc/CMakeLists.txt |  5 +-
 .../analysis_ionization_dsmc_3d.py            | 21 +-----
 .../inputs_test_3d_ionization_dsmc            | 12 +++-
 .../test_3d_ionization_dsmc.json              | 65 ++++++++++---------
 4 files changed, 48 insertions(+), 55 deletions(-)

diff --git a/Examples/Tests/ionization_dsmc/CMakeLists.txt b/Examples/Tests/ionization_dsmc/CMakeLists.txt
index 085902915ff..db8d3cb7b2a 100644
--- a/Examples/Tests/ionization_dsmc/CMakeLists.txt
+++ b/Examples/Tests/ionization_dsmc/CMakeLists.txt
@@ -6,8 +6,7 @@ add_warpx_test(
     3  # dims
     2  # nprocs
     inputs_test_3d_ionization_dsmc  # inputs
-    analysis_ionization_dsmc_3d.py  # analysis
-    diags/diag1/  # output
+    "analysis_ionization_dsmc_3d.py"  # analysis
+    "analysis_default_regression.py --path diags/diag1000250"  # checksum
     OFF  # dependency
 )
-
diff --git a/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py b/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py
index dc0bbee4c0f..38a5b5b8515 100755
--- a/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py
+++ b/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py
@@ -12,12 +12,7 @@
 from pywarpx import picmi
 constants = picmi.constants
 
-sys.path.insert(1, "../../../../warpx/Regression/Checksum/")
-from checksumAPI import evaluate_checksum
-
-filename = sys.argv[1]
-testname = os.path.split(os.getcwd())[1]
-ts = OpenPMDTimeSeries(filename)
+ts = OpenPMDTimeSeries("diags/diag2/")
 
 q_e = constants.q_e
 m_p = constants.m_p
@@ -56,15 +51,13 @@ def get_Te(ts):
     return T_e        
 
 def get_density(ts):
-    number_data = np.genfromtxt("./diags/counts.txt")
+    number_data = np.genfromtxt("diags/counts.txt")
     Te = get_Te(ts)
     total_volume = L[0] * L[1] * L[2]
     electron_weight = number_data[:, 8]
     neutral_weight = number_data[:, 9]
     ne = electron_weight / total_volume
     nn = neutral_weight / total_volume
-    print(np.size(ne), np.size(Te))
-
     return [ne, nn, ne*Te]
 
 def compute_rate_coefficients(temperatures_eV, energy_eV, sigma_m2, num_samples = 1024):
@@ -99,7 +92,6 @@ def compute_rate_coefficients(temperatures_eV, energy_eV, sigma_m2, num_samples
         # get cross section by interpolating on table
         sigma = np.interp(e, energy_eV, sigma_m2, left = 0)
         k[i] = np.mean(sigma * speed)
-    
     return k
 
 def rhs(state, params):
@@ -140,7 +132,6 @@ def rhs(state, params):
     # the opposite is true for the MCC solver
     f[1] = -ndot         # d(nn)/dt
     f[2] = 0             # d(ne*eps) / dt
-
     return f
     
 def solve_theory_model():
@@ -166,7 +157,6 @@ def solve_theory_model():
     ne = state_vec[:, 0]
     nn = state_vec[:, 1]
     Te = state_vec[:, 2] / (1.5 * ne)
-    print(np.size(ne), np.size(Te))
     return t, [ne, nn, ne*Te]
 
 
@@ -218,10 +208,3 @@ def check_tolerance(array, tolerance):
     plt.legend()
     check_tolerance(relative_error, tolerance)
 plt.savefig('./relative_error_density_Te.png', dpi=150) 
-
-# compare checksums
-evaluate_checksum(
-    test_name=os.path.split(os.getcwd())[1],
-    output_file=sys.argv[1],
-    output_format="openpmd",
-)
\ No newline at end of file
diff --git a/Examples/Tests/ionization_dsmc/inputs_test_3d_ionization_dsmc b/Examples/Tests/ionization_dsmc/inputs_test_3d_ionization_dsmc
index 8e430635787..e13156d8a90 100644
--- a/Examples/Tests/ionization_dsmc/inputs_test_3d_ionization_dsmc
+++ b/Examples/Tests/ionization_dsmc/inputs_test_3d_ionization_dsmc
@@ -95,15 +95,21 @@ coll_elec.type = dsmc
 #################################
 ############ DIAGNOSTICS #############
 #################################
-diagnostics.diags_names = diag1
+diagnostics.diags_names = diag1 diag2
 warpx.reduced_diags_names = counts
 counts.intervals = 5
 counts.path = diags/
 counts.type = ParticleNumber
 
 # Diagnostics
-diag1.intervals = 5
+diag1.intervals = 250
 diag1.diag_type = Full
 diag1.electrons.variables = ux uy uz
 diag1.neutrals.variables = ux uy uz
-diag1.format = openpmd
+diag1.format = plotfile
+
+diag2.intervals = 5
+diag2.diag_type = Full
+diag2.electrons.variables = ux uy uz
+diag2.neutrals.variables = ux uy uz
+diag2.format = openpmd
diff --git a/Regression/Checksum/benchmarks_json/test_3d_ionization_dsmc.json b/Regression/Checksum/benchmarks_json/test_3d_ionization_dsmc.json
index 6681b90f1c9..15e1ac7d0b0 100644
--- a/Regression/Checksum/benchmarks_json/test_3d_ionization_dsmc.json
+++ b/Regression/Checksum/benchmarks_json/test_3d_ionization_dsmc.json
@@ -1,33 +1,38 @@
 {
   "lev=0": {
-      "Bx": 0.0,
-      "By": 0.0,
-      "Bz": 0.0,
-      "Ex": 33418.976671097655,
-      "Ey": 28556.630099887625,
-      "Ez": 33908.990756787396,
-      "jx": 351.709319628182,
-      "jy": 389.97631892921305,
-      "jz": 366.54101928747855
-    },
-    "electrons": {
-      "particle_momentum_x": 4.0749126130424784e-20,
-      "particle_momentum_y": 4.086095034273236e-20,
-      "particle_momentum_z": 4.0905976427530184e-20
-    },
-    "ions": {
-      "particle_position_x": 2503.961424209735,
-      "particle_position_y": 2503.383309677753,
-      "particle_position_z": 2505.603131072416,
-      "particle_momentum_x": 1.2028342698649854e-18,
-      "particle_momentum_y": 1.2008761224148308e-18,
-      "particle_momentum_z": 1.2009924309769205e-18,
-      "particle_weight": 152685546875.00003
-    },
-    "neutrals": {
-      "particle_momentum_x": 1.204336106475819e-18,
-      "particle_momentum_y": 1.2126021881554375e-18,
-      "particle_momentum_z": 1.2056811544009856e-18
-    }
+    "Bx": 0.0,
+    "By": 0.0,
+    "Bz": 0.0,
+    "Ex": 33418.97667109626,
+    "Ey": 28556.630099887385,
+    "Ez": 33908.99075678701,
+    "jx": 351.70931962817417,
+    "jy": 389.9763189292032,
+    "jz": 366.5410192874648
+  },
+  "electrons": {
+    "particle_momentum_x": 4.0749126130424784e-20,
+    "particle_momentum_y": 4.086095034273236e-20,
+    "particle_momentum_z": 4.0905976427530184e-20,
+    "particle_position_x": 2505.1614598048145,
+    "particle_position_y": 2509.0808137882786,
+    "particle_position_z": 2500.717012325057
+  },
+  "ions": {
+    "particle_momentum_x": 1.2028342698649854e-18,
+    "particle_momentum_y": 1.2008761224148308e-18,
+    "particle_momentum_z": 1.2009924309769205e-18,
+    "particle_position_x": 2503.961424209735,
+    "particle_position_y": 2503.383309677753,
+    "particle_position_z": 2505.603131072416,
+    "particle_weight": 152685546875.00003
+  },
+  "neutrals": {
+    "particle_momentum_x": 1.2043361064758194e-18,
+    "particle_momentum_y": 1.2126021881554377e-18,
+    "particle_momentum_z": 1.2056811544009856e-18,
+    "particle_position_x": 2503.9591741369045,
+    "particle_position_y": 2503.390353098188,
+    "particle_position_z": 2505.598048390032
   }
-  
\ No newline at end of file
+}

From a9124457a45d7546a128591e9754af24c6521a8d Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Fri, 20 Dec 2024 08:24:05 +0000
Subject: [PATCH 19/25] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 .../ionization_dsmc/analysis_ionization_dsmc_3d.py   | 12 +++++-------
 1 file changed, 5 insertions(+), 7 deletions(-)

diff --git a/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py b/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py
index e176b06461d..ffa8cbd7d17 100755
--- a/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py
+++ b/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py
@@ -1,8 +1,6 @@
 #!/usr/bin/env python3
 # DSMC ionization test script:
 #   - compares WarpX simulation results with theoretical model predictions.
-import os
-import sys
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -63,7 +61,7 @@ def get_density(ts):
     neutral_weight = number_data[:, 9]
     ne = electron_weight / total_volume
     nn = neutral_weight / total_volume
-    return [ne, nn, ne*Te]
+    return [ne, nn, ne * Te]
 
 
 def compute_rate_coefficients(temperatures_eV, energy_eV, sigma_m2, num_samples=1024):
@@ -137,8 +135,8 @@ def rhs(state, params):
 
     # at present, the dsmc solver does not deplete electron energy but does deplete neutrals
     # the opposite is true for the MCC solver
-    f[1] = -ndot         # d(nn)/dt
-    f[2] = 0             # d(ne*eps) / dt
+    f[1] = -ndot  # d(nn)/dt
+    f[2] = 0  # d(ne*eps) / dt
     return f
 
 
@@ -169,7 +167,7 @@ def solve_theory_model():
     ne = state_vec[:, 0]
     nn = state_vec[:, 1]
     Te = state_vec[:, 2] / (1.5 * ne)
-    return t, [ne, nn, ne*Te]
+    return t, [ne, nn, ne * Te]
 
 
 t_warpx = np.loadtxt("diags/counts.txt")[:, 1]
@@ -232,4 +230,4 @@ def check_tolerance(array, tolerance):
     plt.xlabel("Time [s]")
     plt.legend()
     check_tolerance(relative_error, tolerance)
-plt.savefig('./relative_error_density_Te.png', dpi=150) 
+plt.savefig("./relative_error_density_Te.png", dpi=150)

From 22bc06e19395f2fc14f76f930525bba3b9fe8898 Mon Sep 17 00:00:00 2001
From: Olga Shapoval <oshapoval@lbl.gov>
Date: Fri, 20 Dec 2024 09:44:35 -0800
Subject: [PATCH 20/25] Add link to default regression analysis script

---
 Examples/Tests/ionization_dsmc/analysis_default_regression.py | 1 +
 1 file changed, 1 insertion(+)
 create mode 120000 Examples/Tests/ionization_dsmc/analysis_default_regression.py

diff --git a/Examples/Tests/ionization_dsmc/analysis_default_regression.py b/Examples/Tests/ionization_dsmc/analysis_default_regression.py
new file mode 120000
index 00000000000..d8ce3fca419
--- /dev/null
+++ b/Examples/Tests/ionization_dsmc/analysis_default_regression.py
@@ -0,0 +1 @@
+../../analysis_default_regression.py
\ No newline at end of file

From d5c1b778d20caeed3db455fa6912acbebb3f580b Mon Sep 17 00:00:00 2001
From: Olga Shapoval <oshapoval@lbl.gov>
Date: Fri, 20 Dec 2024 09:53:50 -0800
Subject: [PATCH 21/25] Remove duplicate function from analysis script

---
 .../Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py   | 7 -------
 1 file changed, 7 deletions(-)

diff --git a/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py b/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py
index ffa8cbd7d17..600982910ed 100755
--- a/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py
+++ b/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py
@@ -195,13 +195,6 @@ def solve_theory_model():
 plt.savefig("ionization_dsmc__density_Te.png", dpi=150)
 
 
-def check_tolerance(array, tolerance):
-    assert np.all(
-        array <= tolerance
-    ), f"Test did not pass: one or more elements exceed the tolerance of {tolerance}."
-    print("All elements of are within the tolerance.")
-
-
 # For now, setting very high tolerances on purpose. This needs to be modified once all features are implemented in the PR.
 tolerances = [1, 1e-6, 1e-2]
 

From e21b1c566873aa0c820ce004b287ee3966fed847 Mon Sep 17 00:00:00 2001
From: Olga Shapoval <oshapoval@lbl.gov>
Date: Fri, 20 Dec 2024 11:03:36 -0800
Subject: [PATCH 22/25] Updated benchmarks for ionization_dsmc CI test

---
 .../test_3d_ionization_dsmc.json              | 54 +++++++++----------
 1 file changed, 27 insertions(+), 27 deletions(-)

diff --git a/Regression/Checksum/benchmarks_json/test_3d_ionization_dsmc.json b/Regression/Checksum/benchmarks_json/test_3d_ionization_dsmc.json
index 15e1ac7d0b0..110e706df93 100644
--- a/Regression/Checksum/benchmarks_json/test_3d_ionization_dsmc.json
+++ b/Regression/Checksum/benchmarks_json/test_3d_ionization_dsmc.json
@@ -3,36 +3,36 @@
     "Bx": 0.0,
     "By": 0.0,
     "Bz": 0.0,
-    "Ex": 33418.97667109626,
-    "Ey": 28556.630099887385,
-    "Ez": 33908.99075678701,
-    "jx": 351.70931962817417,
-    "jy": 389.9763189292032,
-    "jz": 366.5410192874648
-  },
-  "electrons": {
-    "particle_momentum_x": 4.0749126130424784e-20,
-    "particle_momentum_y": 4.086095034273236e-20,
-    "particle_momentum_z": 4.0905976427530184e-20,
-    "particle_position_x": 2505.1614598048145,
-    "particle_position_y": 2509.0808137882786,
-    "particle_position_z": 2500.717012325057
+    "Ex": 27962.117841006104,
+    "Ey": 31492.65714998622,
+    "Ez": 32685.446517469143,
+    "jx": 371.57329790056076,
+    "jy": 384.5221496512611,
+    "jz": 357.20969486000763
   },
   "ions": {
-    "particle_momentum_x": 1.2028342698649854e-18,
-    "particle_momentum_y": 1.2008761224148308e-18,
-    "particle_momentum_z": 1.2009924309769205e-18,
-    "particle_position_x": 2503.961424209735,
-    "particle_position_y": 2503.383309677753,
-    "particle_position_z": 2505.603131072416,
-    "particle_weight": 152685546875.00003
+    "particle_momentum_x": 1.2038443279955311e-18,
+    "particle_momentum_y": 1.2018846347141862e-18,
+    "particle_momentum_z": 1.200245213220433e-18,
+    "particle_position_x": 2507.9306710321707,
+    "particle_position_y": 2499.1981759893633,
+    "particle_position_z": 2500.8509915859736,
+    "particle_weight": 152996826171.87503
   },
   "neutrals": {
-    "particle_momentum_x": 1.2043361064758194e-18,
-    "particle_momentum_y": 1.2126021881554377e-18,
-    "particle_momentum_z": 1.2056811544009856e-18,
-    "particle_position_x": 2503.9591741369045,
-    "particle_position_y": 2503.390353098188,
-    "particle_position_z": 2505.598048390032
+    "particle_momentum_x": 1.1979391197017039e-18,
+    "particle_momentum_y": 1.2001180842947048e-18,
+    "particle_momentum_z": 1.1996051570661201e-18,
+    "particle_position_x": 2507.928340994915,
+    "particle_position_y": 2499.2207342985084,
+    "particle_position_z": 2500.8614737989697
+  },
+  "electrons": {
+    "particle_momentum_x": 4.103730302973143e-20,
+    "particle_momentum_y": 4.1075862160708524e-20,
+    "particle_momentum_z": 4.107345023993974e-20,
+    "particle_position_x": 2505.02545711238,
+    "particle_position_y": 2501.887490274123,
+    "particle_position_z": 2508.68572175914
   }
 }

From 106a22fa8e054d227730c1d806fa616d1002cba6 Mon Sep 17 00:00:00 2001
From: Remi Lehe <remi.lehe@normalesup.org>
Date: Fri, 20 Dec 2024 14:14:33 -0800
Subject: [PATCH 23/25] Revert CI test

---
 .../inputs_base_1d_picmi.py                   | 59 ++++++++-----------
 1 file changed, 26 insertions(+), 33 deletions(-)

diff --git a/Examples/Physics_applications/capacitive_discharge/inputs_base_1d_picmi.py b/Examples/Physics_applications/capacitive_discharge/inputs_base_1d_picmi.py
index d5e0db7213c..9f2591d8b2b 100644
--- a/Examples/Physics_applications/capacitive_discharge/inputs_base_1d_picmi.py
+++ b/Examples/Physics_applications/capacitive_discharge/inputs_base_1d_picmi.py
@@ -268,40 +268,33 @@ def setup_run(self):
         #######################################################################
 
         cross_sec_direc = "../../../../warpx-data/MCC_cross_sections/He/"
-        electron_scattering_processes = {
-            "elastic": {"cross_section": cross_sec_direc + "electron_scattering.dat"},
-            "excitation1": {
-                "cross_section": cross_sec_direc + "excitation_1.dat",
-                "energy": 19.82,
-            },
-            "excitation2": {
-                "cross_section": cross_sec_direc + "excitation_2.dat",
-                "energy": 20.61,
-            },
-            "ionization": {
-                "cross_section": cross_sec_direc + "ionization.dat",
-                "energy": 24.55,
-                "species": self.ions,
-            },
-        }
 
-        if self.dsmc:
-            electron_colls = picmi.DSMCCollisions(
-                name="coll_elec",
-                species=[self.electrons, self.neutrals],
-                ndt=5,
-                scattering_processes=electron_scattering_processes,
-            )
-        else:
-            electron_colls = picmi.MCCCollisions(
-                name="coll_elec",
-                species=self.electrons,
-                background_density=self.gas_density,
-                background_temperature=self.gas_temp,
-                background_mass=self.ions.mass,
-                ndt=self.mcc_subcycling_steps,
-                scattering_processes=electron_scattering_processes,
-            )
+        electron_colls = picmi.MCCCollisions(
+            name="coll_elec",
+            species=self.electrons,
+            background_density=self.gas_density,
+            background_temperature=self.gas_temp,
+            background_mass=self.ions.mass,
+            ndt=self.mcc_subcycling_steps,
+            scattering_processes={
+                "elastic": {
+                    "cross_section": cross_sec_direc + "electron_scattering.dat"
+                },
+                "excitation1": {
+                    "cross_section": cross_sec_direc + "excitation_1.dat",
+                    "energy": 19.82,
+                },
+                "excitation2": {
+                    "cross_section": cross_sec_direc + "excitation_2.dat",
+                    "energy": 20.61,
+                },
+                "ionization": {
+                    "cross_section": cross_sec_direc + "ionization.dat",
+                    "energy": 24.55,
+                    "species": self.ions,
+                },
+            },
+        )
 
         ion_scattering_processes = {
             "elastic": {"cross_section": cross_sec_direc + "ion_scattering.dat"},

From c51fc58055bd6888e614618e7bc21b6ef7252098 Mon Sep 17 00:00:00 2001
From: Remi Lehe <remi.lehe@normalesup.org>
Date: Fri, 20 Dec 2024 14:34:00 -0800
Subject: [PATCH 24/25] Fix clang issues

---
 .../Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H    | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
index d46b7b56bac..292f5c9c4fa 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
@@ -239,14 +239,14 @@ public:
                 } else if (mask[i] == int(ScatteringProcessType::IONIZATION)) {
                     // calculate kinetic energy of the incident electron
                     double E_coll;
-                    double energy_cost = 0; // TODO: update this: get from scattering process
+                    double const energy_cost = 0; // TODO: update this: get from scattering process
                     const amrex::ParticleReal u_coll2 = ux1*ux1 + uy1*uy1 + uz1*uz1;
                     ParticleUtils::getEnergy(u_coll2, m1, E_coll);
                     // each of the two electron (the scattered e) gets half the energy (could change this later)
-                    const auto E_out = static_cast<amrex::ParticleReal>((E_coll - energy_cost) / 2.0_prt * PhysConst::q_e);
+                    const amrex::ParticleReal E_out = static_cast<amrex::ParticleReal>((E_coll - energy_cost) / 2.0_prt * PhysConst::q_e);
                     constexpr auto c2 = PhysConst::c * PhysConst::c;
                     const auto mc2 = m1*c2;
-                    const amrex::ParticleReal up = sqrt(E_out * (E_out + 2.0_prt*mc2) / c2) / m1;
+                    const amrex::ParticleReal up = static_cast<amrex::ParticleReal>(std::sqrt(E_out * (E_out + 2.0_prt*mc2) / c2) / m1);
 
                     // isotropically scatter electrons
                     // - The incident electron

From 857c99f241c5c9ec86c7cad62e8b093ca9a5278d Mon Sep 17 00:00:00 2001
From: Remi Lehe <remi.lehe@normalesup.org>
Date: Fri, 20 Dec 2024 14:49:38 -0800
Subject: [PATCH 25/25] Remove the created neutral particle

---
 .../Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H   | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
index 292f5c9c4fa..6ade310ee69 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
@@ -145,7 +145,6 @@ public:
                 // Set the weight of the new particles to p_pair_reaction_weight[i]
                 soa_products_data[0].m_rdata[PIdx::w][product1_index] = p_pair_reaction_weight[i];
 
-                // TODO: when the reaction is ionization, this still creates a second particle of the target...
                 const auto product2_index = products_np_data[1] +
                                            (p_offsets[i]*p_num_products_device[1] + 0);
                 // Make a copy of the particle from species 2
@@ -171,6 +170,12 @@ public:
                                         static_cast<int>(product3_index), engine);
                         // Set the weight of the new particles to p_pair_reaction_weight[i]
                         soa_products_data[2].m_rdata[PIdx::w][product3_index] = p_pair_reaction_weight[i];
+
+                        // The ion created in product 3 replaces the scattered neutral particle in product 2
+                        // (which would have been created by other, elastic scattering processes)
+                        // Therefore, we set the ID of the scattered neutral particle to invalid,
+                        // so that it gets removed in the next call to `ReDistribute` or `RemoveInvalidParticles`
+                        soa_products_data[1].idcpu(product2_index) = amrex::ParticleIdCpus::Invalid;
                     } else {
                         // Set IDs as invalid for slots that have been allocated for potential
                         // ionization events, so that the corresponding particles get removed in the next