From 183cf2a124274e6689c2512630df8a051ad65b5c Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 12:08:03 +0100
Subject: [PATCH 01/37] Add method to calculate equilibration time in plasma

---
 process/models/physics/physics.py | 55 +++++++++++++++++++++++++++++++
 1 file changed, 55 insertions(+)

diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 413e545d4c..d5679edad7 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -5937,6 +5937,61 @@ def calculate_spitzer_resistivity(
             / (temp_plasma_electron_kev * constants.KILOELECTRON_VOLT) ** 1.5
         )
 
+    @staticmethod
+    @nb.njit(cache=True)
+    def calculate_equilibriation_time(
+        temp_plasma_electron_kev: float | np.ndarray,
+        nd_plasma_ions: float | np.ndarray,
+        plasma_coulomb_log_electron_ion: float | np.ndarray,
+        m_ion: float,
+        n_charge_ion: float = 1.0,
+    ):
+        """
+        Calculate the equilibration time (τ_eq) between electrons and ions in a plasma.
+
+        Parameters
+        ----------
+        temp_plasma_electron_kev : float | np.ndarray
+            Electron temperature in keV.
+        nd_plasma_ions : float | np.ndarray
+            Ion density (m^-3).
+        plasma_coulomb_log_electron_ion : float | np.ndarray
+            Coulomb logarithm for electron-ion collisions.
+        m_ion : float
+            Ion mass (kg).
+        n_charge_ion : float, optional
+            Charge number (Z) of the ion. Default is 1.0.
+
+        Returns
+        -------
+        float | np.ndarray
+            Equilibration time (s).
+
+        Notes
+        -----
+        - The equilibration time is the characteristic time for energy exchange between
+        electrons and ions in a plasma, leading to thermal equilibrium.
+        - It is the characteristic timescale required for two different particle species
+        to share heat and reach thermal equilibrium.
+        """
+        return (
+            3
+            * (2 * np.pi) ** 1.5
+            * constants.EPSILON0**2
+            * m_ion
+            * (constants.ELECTRON_MASS)
+            * (
+                (temp_plasma_electron_kev * constants.KILOELECTRON_VOLT)
+                / constants.ELECTRON_MASS
+            )
+            ** 1.5
+        ) / (
+            nd_plasma_ions
+            * n_charge_ion**2
+            * plasma_coulomb_log_electron_ion
+            * constants.ELECTRON_CHARGE**4
+        )
+
     def output_detailed_physics(self):
         """Outputs detailed physics variables to file."""
         po.oheadr(self.outfile, "Detailed Plasma")

From 64027f7831437a0e10fd0c5ee7c9bc743e0e61fc Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 12:11:53 +0100
Subject: [PATCH 02/37] Add electron-deuteron and electron-triton equilibration
 time variables to PhysicsData

---
 process/data_structure/physics_variables.py | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 28e37dee79..6cbd91e779 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1581,6 +1581,22 @@ class PhysicsData:
     t_plasma_electron_alpha_thermal_collision_vol_avg: float = 0.0
     """Volume averaged electron-alpha collision time in plasma (s)"""
 
+    t_plasma_electron_deuteron_equilbriation_vol_avg: float = 0.0
+    """Volume averaged electron-deuteron equilibration time in plasma (s)"""
+
+    t_plasma_electron_deuteron_equilbriation_profile: list[float] = field(
+        default_factory=list
+    )
+    """Profile of electron-deuteron equilibration time in plasma (s)"""
+
+    t_plasma_electron_triton_equilbriation_vol_avg: float = 0.0
+    """Volume averaged electron-triton equilibration time in plasma (s)"""
+
+    t_plasma_electron_triton_equilbriation_profile: list[float] = field(
+        default_factory=list
+    )
+    """Profile of electron-triton equilibration time in plasma (s)"""
+
     freq_plasma_electron_electron_collision_profile: list[float] = field(
         default_factory=list
     )

From 2feee91122ce49ee5e3f07bbfdeefde1a88e9412 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 12:58:26 +0100
Subject: [PATCH 03/37] Add plasma electron-deuteron equilibration time
 calculations and plotting

---
 process/core/io/plot/summary.py   | 45 +++++++++++++++++++++++++++
 process/models/physics/physics.py | 51 +++++++++++++++++++++++++++++++
 2 files changed, 96 insertions(+)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index e8b2fea071..c57afd9092 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -14090,6 +14090,33 @@ def plot_resistivity_profile(axis: plt.Axes, mfile_data: MFile, scan: int) -> No
     axis.legend()
 
 
+def plot_plasma_equilibration_time_profile(
+    axis: plt.Axes, mfile_data: MFile, scan: int
+) -> None:
+    """Plot the plasma equilibration time on the given axis."""
+    t_plasma_electron_deuteron_equilbriation_profile = [
+        mfile_data.data[f"t_plasma_electron_deuteron_equilbriation_profile{i}"].get_scan(
+            scan
+        )
+        for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)))
+    ]
+
+    axis.plot(
+        np.linspace(0, 1, len(t_plasma_electron_deuteron_equilbriation_profile)),
+        t_plasma_electron_deuteron_equilbriation_profile,
+        color="red",
+        linestyle="-",
+        label=r"$\tau_{e-D,eq}$",
+    )
+
+    axis.set_yscale("log")
+    axis.set_ylabel("Equilibration Time [s]")
+    axis.set_xlabel("$\\rho \\ [r/a]$")
+    axis.grid(True, which="both", linestyle="--", alpha=0.5)
+    axis.minorticks_on()
+    axis.legend()
+
+
 def plot_equality_constraint_equations(axis: plt.Axes, m_file_data: MFile, scan: int):
     """Plot the equality constraints for a solution and their normalised residuals
 
@@ -14785,6 +14812,11 @@ def plot_detailed_plasma_parameters(axis: plt.Axes, fig, mfile: MFile, scan: int
         f"$\\langle\\eta_{{Spitzer}}\\rangle$: {mfile.get('res_plasma_fuel_spitzer_vol_avg', scan=scan):.4e} $\\Omega\\mathrm{{m}}$"
     )
 
+    textstr_equilibriation = (
+        f"$\\mathbf{{Equilibration \\ Times:}}$\n\n"
+        f"$\\langle\\tau_{{e-D,Equil}}\\rangle$: {mfile.get('t_plasma_electron_deuteron_equilibration_vol_avg', scan=scan):.4e} s\n"
+    )
+
     light_yellow_box = {
         "boxstyle": "round",
         "facecolor": "lightyellow",
@@ -14825,6 +14857,17 @@ def plot_detailed_plasma_parameters(axis: plt.Axes, fig, mfile: MFile, scan: int
         bbox=light_yellow_box,
     )
 
+    axis.text(
+        0.65,
+        0.45,
+        textstr_equilibriation,
+        fontsize=9,
+        verticalalignment="top",
+        horizontalalignment="left",
+        transform=fig.transFigure,
+        bbox=light_yellow_box,
+    )
+
     light_cyan_box = {
         "boxstyle": "round",
         "facecolor": "lightcyan",
@@ -16114,6 +16157,8 @@ def _add_page(name: str | None = None):
 
     plot_resistivity_profile(pages["detailed_params"].add_subplot(232), m_file, scan)
 
+    plot_plasma_equilibration_time_profile(pages["detailed_params"].add_subplot(233), m_file, scan)
+
     plot_detailed_plasma_parameters(
         pages["detailed_params"].add_subplot(233),
         fig=pages["detailed_params"],
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index d5679edad7..529e10afec 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -5495,6 +5495,38 @@ def run(self):
             temp_plasma_electron_kev=self.plasma_profile.teprofile.profile_y,
         )
 
+        # ================================
+        # Ion-electron equilibration times
+        # ================================
+
+        self.data.physics.t_plasma_electron_deuteron_equilbriation_vol_avg = self.calculate_equilibriation_time(
+            temp_plasma_electron_kev=self.data.physics.temp_plasma_electron_vol_avg_kev,
+            nd_plasma_ions=self.data.physics.nd_plasma_electrons_vol_avg
+            * self.data.physics.f_plasma_fuel_deuterium
+            * (
+                self.data.physics.nd_plasma_fuel_ions_vol_avg
+                / self.data.physics.nd_plasma_electrons_vol_avg
+            ),
+            plasma_coulomb_log_electron_ion=self.data.physics.plasma_coulomb_log_electron_deuteron_vol_avg,
+            m_ion=constants.DEUTERON_MASS,
+            n_charge_ion=1,
+        )
+
+        self.data.physics.t_plasma_electron_deuteron_equilbriation_profile = self.calculate_equilibriation_time(
+            temp_plasma_electron_kev=self.plasma_profile.teprofile.profile_y,
+            nd_plasma_ions=(
+                self.plasma_profile.neprofile.profile_y
+                * self.data.physics.f_plasma_fuel_deuterium
+                * (
+                    self.data.physics.nd_plasma_fuel_ions_vol_avg
+                    / self.data.physics.nd_plasma_electrons_vol_avg
+                )
+            ),
+            plasma_coulomb_log_electron_ion=self.data.physics.plasma_coulomb_log_electron_deuteron_profile,
+            m_ion=constants.DEUTERON_MASS,
+            n_charge_ion=1,
+        )
+
     @staticmethod
     @nb.njit(cache=True)
     def calculate_debye_length(
@@ -6550,3 +6582,22 @@ def output_detailed_physics(self):
                 f"(res_plasma_fuel_spitzer_profile{i})",
                 self.data.physics.res_plasma_fuel_spitzer_profile[i],
             )
+
+        po.osubhd(self.outfile, "Equilibration Times:")
+
+        po.ovarre(
+            self.outfile,
+            "Volume averaged electron-deuteron equilibration time (τ_eq) (s)",
+            "(t_plasma_electron_deuteron_equilbriation_vol_avg)",
+            self.data.physics.t_plasma_electron_deuteron_equilbriation_vol_avg,
+        )
+
+        for i in range(
+            len(self.data.physics.t_plasma_electron_deuteron_equilbriation_profile)
+        ):
+            po.ovarre(
+                self.mfile,
+                f"Electron-deuteron equilibration time at point {i}",
+                f"(t_plasma_electron_deuteron_equilbriation_profile{i})",
+                self.data.physics.t_plasma_electron_deuteron_equilbriation_profile[i],
+            )

From 95cebec542ac60ba5b86bd9e7558c9ef1850b355 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 13:07:32 +0100
Subject: [PATCH 04/37] Fix naming of equilibration error

---
 process/core/io/plot/summary.py             |  8 ++++----
 process/data_structure/physics_variables.py |  8 ++++----
 process/models/physics/physics.py           | 14 +++++++-------
 3 files changed, 15 insertions(+), 15 deletions(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index c57afd9092..a2ca178f03 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -14094,16 +14094,16 @@ def plot_plasma_equilibration_time_profile(
     axis: plt.Axes, mfile_data: MFile, scan: int
 ) -> None:
     """Plot the plasma equilibration time on the given axis."""
-    t_plasma_electron_deuteron_equilbriation_profile = [
-        mfile_data.data[f"t_plasma_electron_deuteron_equilbriation_profile{i}"].get_scan(
+    t_plasma_electron_deuteron_equilibration_profile = [
+        mfile_data.data[f"t_plasma_electron_deuteron_equilibration_profile{i}"].get_scan(
             scan
         )
         for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)))
     ]
 
     axis.plot(
-        np.linspace(0, 1, len(t_plasma_electron_deuteron_equilbriation_profile)),
-        t_plasma_electron_deuteron_equilbriation_profile,
+        np.linspace(0, 1, len(t_plasma_electron_deuteron_equilibration_profile)),
+        t_plasma_electron_deuteron_equilibration_profile,
         color="red",
         linestyle="-",
         label=r"$\tau_{e-D,eq}$",
diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 6cbd91e779..6187559e8b 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1581,18 +1581,18 @@ class PhysicsData:
     t_plasma_electron_alpha_thermal_collision_vol_avg: float = 0.0
     """Volume averaged electron-alpha collision time in plasma (s)"""
 
-    t_plasma_electron_deuteron_equilbriation_vol_avg: float = 0.0
+    t_plasma_electron_deuteron_equilibration_vol_avg: float = 0.0
     """Volume averaged electron-deuteron equilibration time in plasma (s)"""
 
-    t_plasma_electron_deuteron_equilbriation_profile: list[float] = field(
+    t_plasma_electron_deuteron_equilibration_profile: list[float] = field(
         default_factory=list
     )
     """Profile of electron-deuteron equilibration time in plasma (s)"""
 
-    t_plasma_electron_triton_equilbriation_vol_avg: float = 0.0
+    t_plasma_electron_triton_equilibration_vol_avg: float = 0.0
     """Volume averaged electron-triton equilibration time in plasma (s)"""
 
-    t_plasma_electron_triton_equilbriation_profile: list[float] = field(
+    t_plasma_electron_triton_equilibration_profile: list[float] = field(
         default_factory=list
     )
     """Profile of electron-triton equilibration time in plasma (s)"""
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 529e10afec..8cb114df14 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -5499,7 +5499,7 @@ def run(self):
         # Ion-electron equilibration times
         # ================================
 
-        self.data.physics.t_plasma_electron_deuteron_equilbriation_vol_avg = self.calculate_equilibriation_time(
+        self.data.physics.t_plasma_electron_deuteron_equilibration_vol_avg = self.calculate_equilibriation_time(
             temp_plasma_electron_kev=self.data.physics.temp_plasma_electron_vol_avg_kev,
             nd_plasma_ions=self.data.physics.nd_plasma_electrons_vol_avg
             * self.data.physics.f_plasma_fuel_deuterium
@@ -5512,7 +5512,7 @@ def run(self):
             n_charge_ion=1,
         )
 
-        self.data.physics.t_plasma_electron_deuteron_equilbriation_profile = self.calculate_equilibriation_time(
+        self.data.physics.t_plasma_electron_deuteron_equilibration_profile = self.calculate_equilibriation_time(
             temp_plasma_electron_kev=self.plasma_profile.teprofile.profile_y,
             nd_plasma_ions=(
                 self.plasma_profile.neprofile.profile_y
@@ -6588,16 +6588,16 @@ def output_detailed_physics(self):
         po.ovarre(
             self.outfile,
             "Volume averaged electron-deuteron equilibration time (τ_eq) (s)",
-            "(t_plasma_electron_deuteron_equilbriation_vol_avg)",
-            self.data.physics.t_plasma_electron_deuteron_equilbriation_vol_avg,
+            "(t_plasma_electron_deuteron_equilibration_vol_avg)",
+            self.data.physics.t_plasma_electron_deuteron_equilibration_vol_avg,
         )
 
         for i in range(
-            len(self.data.physics.t_plasma_electron_deuteron_equilbriation_profile)
+            len(self.data.physics.t_plasma_electron_deuteron_equilibration_profile)
         ):
             po.ovarre(
                 self.mfile,
                 f"Electron-deuteron equilibration time at point {i}",
-                f"(t_plasma_electron_deuteron_equilbriation_profile{i})",
-                self.data.physics.t_plasma_electron_deuteron_equilbriation_profile[i],
+                f"(t_plasma_electron_deuteron_equilibration_profile{i})",
+                self.data.physics.t_plasma_electron_deuteron_equilibration_profile[i],
             )

From 3ddf89348f53fcbecbedb43f6fea5a145f519fcc Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 14:12:35 +0100
Subject: [PATCH 05/37] Add electron-triton equilibration time calculations and
 plotting

---
 process/core/io/plot/summary.py   | 18 ++++++++++++-
 process/models/physics/physics.py | 45 +++++++++++++++++++++++++++++++
 2 files changed, 62 insertions(+), 1 deletion(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index a2ca178f03..57ff833f84 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -14101,14 +14101,29 @@ def plot_plasma_equilibration_time_profile(
         for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)))
     ]
 
+    t_plasma_electron_triton_equilibration_profile = [
+        mfile_data.data[f"t_plasma_electron_triton_equilibration_profile{i}"].get_scan(
+            scan
+        )
+        for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)))
+    ]
+
     axis.plot(
         np.linspace(0, 1, len(t_plasma_electron_deuteron_equilibration_profile)),
         t_plasma_electron_deuteron_equilibration_profile,
-        color="red",
+        color="pink",
         linestyle="-",
         label=r"$\tau_{e-D,eq}$",
     )
 
+    axis.plot(
+        np.linspace(0, 1, len(t_plasma_electron_triton_equilibration_profile)),
+        t_plasma_electron_triton_equilibration_profile,
+        color="green",
+        linestyle="-",
+        label=r"$\tau_{e-T,eq}$",
+    )
+
     axis.set_yscale("log")
     axis.set_ylabel("Equilibration Time [s]")
     axis.set_xlabel("$\\rho \\ [r/a]$")
@@ -14815,6 +14830,7 @@ def plot_detailed_plasma_parameters(axis: plt.Axes, fig, mfile: MFile, scan: int
     textstr_equilibriation = (
         f"$\\mathbf{{Equilibration \\ Times:}}$\n\n"
         f"$\\langle\\tau_{{e-D,Equil}}\\rangle$: {mfile.get('t_plasma_electron_deuteron_equilibration_vol_avg', scan=scan):.4e} s\n"
+        f"$\\langle\\tau_{{e-T,Equil}}\\rangle$: {mfile.get('t_plasma_electron_triton_equilibration_vol_avg', scan=scan):.4e} s\n"
     )
 
     light_yellow_box = {
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 8cb114df14..60baeb66c7 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -5527,6 +5527,34 @@ def run(self):
             n_charge_ion=1,
         )
 
+        self.data.physics.t_plasma_electron_triton_equilibration_vol_avg = self.calculate_equilibriation_time(
+            temp_plasma_electron_kev=self.data.physics.temp_plasma_electron_vol_avg_kev,
+            nd_plasma_ions=self.data.physics.nd_plasma_electrons_vol_avg
+            * self.data.physics.f_plasma_fuel_tritium
+            * (
+                self.data.physics.nd_plasma_fuel_ions_vol_avg
+                / self.data.physics.nd_plasma_electrons_vol_avg
+            ),
+            plasma_coulomb_log_electron_ion=self.data.physics.plasma_coulomb_log_electron_triton_vol_avg,
+            m_ion=constants.TRITON_MASS,
+            n_charge_ion=1,
+        )
+
+        self.data.physics.t_plasma_electron_triton_equilibration_profile = self.calculate_equilibriation_time(
+            temp_plasma_electron_kev=self.plasma_profile.teprofile.profile_y,
+            nd_plasma_ions=(
+                self.plasma_profile.neprofile.profile_y
+                * self.data.physics.f_plasma_fuel_tritium
+                * (
+                    self.data.physics.nd_plasma_fuel_ions_vol_avg
+                    / self.data.physics.nd_plasma_electrons_vol_avg
+                )
+            ),
+            plasma_coulomb_log_electron_ion=self.data.physics.plasma_coulomb_log_electron_triton_profile,
+            m_ion=constants.TRITON_MASS,
+            n_charge_ion=1,
+        )
+
     @staticmethod
     @nb.njit(cache=True)
     def calculate_debye_length(
@@ -6601,3 +6629,20 @@ def output_detailed_physics(self):
                 f"(t_plasma_electron_deuteron_equilibration_profile{i})",
                 self.data.physics.t_plasma_electron_deuteron_equilibration_profile[i],
             )
+
+        po.ovarre(
+            self.outfile,
+            "Volume averaged electron-triton equilibration time (τ_eq) (s)",
+            "(t_plasma_electron_triton_equilibration_vol_avg)",
+            self.data.physics.t_plasma_electron_triton_equilibration_vol_avg,
+        )
+
+        for i in range(
+            len(self.data.physics.t_plasma_electron_triton_equilibration_profile)
+        ):
+            po.ovarre(
+                self.mfile,
+                f"Electron-triton equilibration time at point {i}",
+                f"(t_plasma_electron_triton_equilibration_profile{i})",
+                self.data.physics.t_plasma_electron_triton_equilibration_profile[i],
+            )

From 3f443522f99eef215ef42de0b84b81acdaf9ef03 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 14:15:43 +0100
Subject: [PATCH 06/37] Rename rether to
 calculate_ion_electron_equilibration_power and update references

---
 process/models/physics/physics.py         | 20 +++++++++++---------
 process/models/stellarator/stellarator.py | 23 ++++++++++++++---------
 tests/unit/models/physics/test_physics.py |  6 +++---
 3 files changed, 28 insertions(+), 21 deletions(-)

diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 60baeb66c7..1ce5a1fdd2 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -98,7 +98,7 @@ def calculate_cylindrical_safety_factor(
 
 
 @nb.jit(nopython=True, cache=True)
-def rether(
+def calculate_ion_electron_equilibration_power(
     alphan,
     alphat,
     nd_plasma_electrons_vol_avg,
@@ -707,14 +707,16 @@ def run(self):
 
         # Calculate ion/electron equilibration power
 
-        self.data.physics.pden_ion_electron_equilibration_mw = rether(
-            self.data.physics.alphan,
-            self.data.physics.alphat,
-            self.data.physics.nd_plasma_electrons_vol_avg,
-            self.data.physics.dlamie,
-            self.data.physics.temp_plasma_electron_vol_avg_kev,
-            self.data.physics.temp_plasma_ion_vol_avg_kev,
-            self.data.physics.n_charge_plasma_effective_mass_weighted_vol_avg,
+        self.data.physics.pden_ion_electron_equilibration_mw = (
+            calculate_ion_electron_equilibration_power(
+                self.data.physics.alphan,
+                self.data.physics.alphat,
+                self.data.physics.nd_plasma_electrons_vol_avg,
+                self.data.physics.dlamie,
+                self.data.physics.temp_plasma_electron_vol_avg_kev,
+                self.data.physics.temp_plasma_ion_vol_avg_kev,
+                self.data.physics.n_charge_plasma_effective_mass_weighted_vol_avg,
+            )
         )
 
         # Calculate radiation power
diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py
index 4d23487128..c5e699c126 100644
--- a/process/models/stellarator/stellarator.py
+++ b/process/models/stellarator/stellarator.py
@@ -15,7 +15,10 @@
 from process.core.model import Model
 from process.data_structure.physics_variables import PlasmaIgnitionModel
 from process.models.engineering.pumping import CoolantType
-from process.models.physics.physics import Physics, rether
+from process.models.physics.physics import (
+    Physics,
+    calculate_ion_electron_equilibration_power,
+)
 from process.models.power import PumpingPowerModelTypes
 from process.models.stellarator.build import st_build
 from process.models.stellarator.coils.calculate import st_coil
@@ -2104,14 +2107,16 @@ def st_phys(self, output):
 
         # Calculate ion/electron equilibration power
 
-        self.data.physics.pden_ion_electron_equilibration_mw = rether(
-            self.data.physics.alphan,
-            self.data.physics.alphat,
-            self.data.physics.nd_plasma_electrons_vol_avg,
-            self.data.physics.dlamie,
-            self.data.physics.temp_plasma_electron_vol_avg_kev,
-            self.data.physics.temp_plasma_ion_vol_avg_kev,
-            self.data.physics.n_charge_plasma_effective_mass_weighted_vol_avg,
+        self.data.physics.pden_ion_electron_equilibration_mw = (
+            calculate_ion_electron_equilibration_power(
+                self.data.physics.alphan,
+                self.data.physics.alphat,
+                self.data.physics.nd_plasma_electrons_vol_avg,
+                self.data.physics.dlamie,
+                self.data.physics.temp_plasma_electron_vol_avg_kev,
+                self.data.physics.temp_plasma_ion_vol_avg_kev,
+                self.data.physics.n_charge_plasma_effective_mass_weighted_vol_avg,
+            )
         )
 
         # Calculate radiation power
diff --git a/tests/unit/models/physics/test_physics.py b/tests/unit/models/physics/test_physics.py
index ccef3e81db..053514389f 100644
--- a/tests/unit/models/physics/test_physics.py
+++ b/tests/unit/models/physics/test_physics.py
@@ -14,9 +14,9 @@
     PlasmaBeta,
     PlasmaInductance,
     calculate_cylindrical_safety_factor,
+    calculate_ion_electron_equilibration_power,
     ps_fraction_scene,
     res_diff_time,
-    rether,
 )
 
 
@@ -1970,8 +1970,8 @@ def test_phyaux(phyauxparam, monkeypatch, physics):
     assert t_alpha_confinement == pytest.approx(phyauxparam.expected_t_alpha_confinement)
 
 
-def test_rether():
-    assert rether(
+def test_calculate_ion_electron_equilibration_power():
+    assert calculate_ion_electron_equilibration_power(
         1.0, 1.45, 7.5e19, 17.81065204, 12.0, 13.0, 0.43258985
     ) == pytest.approx(0.028360489673597476)
 

From edc3b07a5c036f4406bc2224f753498f66fd2a23 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 14:33:47 +0100
Subject: [PATCH 07/37] Update electron-deuteron and electron-triton
 equilibration time variables and adjust calculation function

---
 process/data_structure/physics_variables.py |  4 +-
 process/models/physics/physics.py           | 82 +++++++++------------
 tests/unit/models/physics/test_physics.py   |  7 +-
 3 files changed, 41 insertions(+), 52 deletions(-)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 6187559e8b..007b00520d 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1581,7 +1581,7 @@ class PhysicsData:
     t_plasma_electron_alpha_thermal_collision_vol_avg: float = 0.0
     """Volume averaged electron-alpha collision time in plasma (s)"""
 
-    t_plasma_electron_deuteron_equilibration_vol_avg: float = 0.0
+    t_plasma_electron_deuteron_equilibration_vol_avg: float = 3.0
     """Volume averaged electron-deuteron equilibration time in plasma (s)"""
 
     t_plasma_electron_deuteron_equilibration_profile: list[float] = field(
@@ -1589,7 +1589,7 @@ class PhysicsData:
     )
     """Profile of electron-deuteron equilibration time in plasma (s)"""
 
-    t_plasma_electron_triton_equilibration_vol_avg: float = 0.0
+    t_plasma_electron_triton_equilibration_vol_avg: float = 4.5
     """Volume averaged electron-triton equilibration time in plasma (s)"""
 
     t_plasma_electron_triton_equilibration_profile: list[float] = field(
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 1ce5a1fdd2..248e3c7f94 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -99,56 +99,41 @@ def calculate_cylindrical_safety_factor(
 
 @nb.jit(nopython=True, cache=True)
 def calculate_ion_electron_equilibration_power(
-    alphan,
-    alphat,
-    nd_plasma_electrons_vol_avg,
-    dlamie,
-    te,
-    temp_plasma_ion_vol_avg_kev,
-    n_charge_plasma_effective_mass_weighted_vol_avg,
-):
+    nd_electrons: float | np.ndarray,
+    temp_plasma_electron_kev: float | np.ndarray,
+    temp_plasma_ion_kev: float | np.ndarray,
+    t_plasma_ion_electron_equilibration: float | np.ndarray,
+) -> float | np.ndarray:
     """Routine to find the equilibration power between the
-    ions and electrons
-    This routine calculates the equilibration power between the
-    ions and electrons.
-    Unknown origin
+    ions and electrons (Pₑᵢ) in the plasma.
 
     Parameters
     ----------
-    alphan :
-        density profile index
-    alphat :
-        temperature profile index
-    nd_plasma_electrons_vol_avg :
-        electron density (/m3)
-    dlamie :
-        ion-electron coulomb logarithm
-    te :
-        electron temperature (keV)
-    temp_plasma_ion_vol_avg_kev :
-        ion temperature (keV)
-    n_charge_plasma_effective_mass_weighted_vol_avg :
-        mass weighted plasma effective charge
+    nd_electrons : float | np.ndarray
+        Electron density [/m³]
+    temp_plasma_electron_kev : float | np.ndarray
+        Electron temperature [keV]
+    temp_plasma_ion_kev : float | np.ndarray
+        Ion temperature [keV]
+    t_plasma_ion_electron_equilibration : float | np.ndarray
+        Ion-electron equilibration time [s]
 
     Returns
     -------
-    pden_ion_electron_equilibration_mw  :
-        ion/electron equilibration power (MW/m3)
+    pden_ion_electron_equilibration  :
+        ion/electron equilibration power (Pₑᵢ) [W/m³]
 
     """
-    profie = (1.0 + alphan) ** 2 / (
-        (2.0 * alphan - 0.5 * alphat + 1.0) * np.sqrt(1.0 + alphat)
-    )
-    conie = (
-        2.42165e-41
-        * dlamie
-        * nd_plasma_electrons_vol_avg**2
-        * n_charge_plasma_effective_mass_weighted_vol_avg
-        * profie
+    return (
+        (3 / 2)
+        * nd_electrons
+        * constants.KILOELECTRON_VOLT
+        * (
+            (temp_plasma_electron_kev - temp_plasma_ion_kev)
+            / t_plasma_ion_electron_equilibration
+        )
     )
 
-    return conie * (temp_plasma_ion_vol_avg_kev - te) / (te**1.5)
-
 
 # -----------------------------------------------------
 # Diamagnetic and Pfirsch-Schlüter Current Calculations
@@ -707,16 +692,17 @@ def run(self):
 
         # Calculate ion/electron equilibration power
 
+        pden_ion_electron_equilibration = calculate_ion_electron_equilibration_power(
+            nd_electrons=self.data.physics.nd_plasma_electrons_vol_avg,
+            temp_plasma_electron_kev=self.data.physics.temp_plasma_electron_vol_avg_kev,
+            temp_plasma_ion_kev=self.data.physics.temp_plasma_ion_vol_avg_kev,
+            t_plasma_ion_electron_equilibration=np.average([
+                self.data.physics.t_plasma_electron_deuteron_equilibration_vol_avg,
+                self.data.physics.t_plasma_electron_triton_equilibration_vol_avg,
+            ]),
+        )
         self.data.physics.pden_ion_electron_equilibration_mw = (
-            calculate_ion_electron_equilibration_power(
-                self.data.physics.alphan,
-                self.data.physics.alphat,
-                self.data.physics.nd_plasma_electrons_vol_avg,
-                self.data.physics.dlamie,
-                self.data.physics.temp_plasma_electron_vol_avg_kev,
-                self.data.physics.temp_plasma_ion_vol_avg_kev,
-                self.data.physics.n_charge_plasma_effective_mass_weighted_vol_avg,
-            )
+            pden_ion_electron_equilibration / 1.0e6
         )
 
         # Calculate radiation power
diff --git a/tests/unit/models/physics/test_physics.py b/tests/unit/models/physics/test_physics.py
index 053514389f..5aa70f4689 100644
--- a/tests/unit/models/physics/test_physics.py
+++ b/tests/unit/models/physics/test_physics.py
@@ -1972,8 +1972,11 @@ def test_phyaux(phyauxparam, monkeypatch, physics):
 
 def test_calculate_ion_electron_equilibration_power():
     assert calculate_ion_electron_equilibration_power(
-        1.0, 1.45, 7.5e19, 17.81065204, 12.0, 13.0, 0.43258985
-    ) == pytest.approx(0.028360489673597476)
+        nd_electrons=7.5e19,
+        temp_plasma_electron_kev=17.81065204,
+        temp_plasma_ion_kev=12.0,
+        t_plasma_ion_electron_equilibration=0.43258985,
+    ) == pytest.approx(242109.29342520377)
 
 
 class PohmParam(NamedTuple):

From d1c882f8f81208ba4d9f95f2f0ccb549a818a762 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 14:44:35 +0100
Subject: [PATCH 08/37] Add f_temp_plasma_ion_electron iteration variable for
 physics calculations

---
 process/core/solver/iteration_variables.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py
index ac79dc2789..61c07979c3 100644
--- a/process/core/solver/iteration_variables.py
+++ b/process/core/solver/iteration_variables.py
@@ -235,6 +235,7 @@ class IterationVariable:
     175: IterationVariable("kappa", "physics", 0.00, 10.00),
     176: IterationVariable("f_st_coil_aspect", "stellarator", 0.70, 1.30),
     177: IterationVariable("f_a_tf_turn_cable_space_extra_void", "tfcoil", 0.01, 1.0),
+    178: IterationVariable("f_temp_plasma_ion_electron", "physics", 0.0, 2.0),
 }
 
 

From 23d19ad6ee04e17863ec708d975eaea6f46ac4ea Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 15:36:15 +0100
Subject: [PATCH 09/37] Refactor plasma power balance documentation for clarity
 and completeness

---
 .../physics-models/plasma_power_balance.md    | 102 ++++++++++++++----
 1 file changed, 82 insertions(+), 20 deletions(-)

diff --git a/documentation/source/physics-models/plasma_power_balance.md b/documentation/source/physics-models/plasma_power_balance.md
index cdab37b71c..791640435a 100644
--- a/documentation/source/physics-models/plasma_power_balance.md
+++ b/documentation/source/physics-models/plasma_power_balance.md
@@ -1,20 +1,82 @@
-# Plasma Core Power Balance
-
-The figure below shows the flow of power as calculated by the code.
-
-<figure markdown>
-![Power balance](../images/Overall-power-flow.png){ width="100%"}
-<figcaption>Figure 1: Machine build for D-shaped major components</figcaption>
-</figure>
-
-The primary sources of power are the fusion reactions themselves, ohmic power
-due to resistive heating within the plasma, and any auxiliary power provided
-for heating and current drive. The power carried by the fusion-generated 
-neutrons is lost from the plasma, but is deposited in the surrounding material. 
-A fraction `f_p_alpha_plasma_deposited` of the alpha particle power is assumed to stay within the
-plasma core to contribute to the plasma power balance. The sum of this core
-alpha power, any power carried by non-alpha charged particles, the ohmic power
-and any injected power, is converted into charged particle transport power
-($P_{\mbox{loss}}$) plus core radiation power, as shown in the Figure. The core 
-power balance calculation is turned on using constraint equation no. 2 (which 
-should therefore always be used).
\ No newline at end of file
+# Plasma Power Balance
+
+There are several equations that need to be satisfied in order to ensure that the plasma is in 
+power equilibrium. This includes the total energy and power leaving the plasma as a whole and then also the rate of energy transfer between the ion and electron species themselves.
+
+!!! warning "Treatment of particle energy confinement"
+
+    Currently the electrons and ions in `PROCESS` are assumed to have the same energy confinement time and which is equal to the calculated global energy confinement time, $\tau_{e} = \tau_{i} = \tau_{E}$.
+
+## Global plasma power balance
+
+This constraint can be activated by stating `icc = 2` in the input file.
+
+This constraint ensures self consistency between the the transport loss power used for the confinement scalings and the calculated confinement time in relation to the plasmas total thermal energy:
+
+$$
+P_{\text{L}} = \frac{W}{\tau_{\text{E}}}
+$$
+
+$$
+\underbrace{\frac{3}{2}\frac{n_{\text{i}} \langle T_{\text{i}} \rangle_{\text{n}}}{\tau_{\text{E}}} + \frac{3}{2}\frac{n_{\text{e}} \langle T_{\text{e}} \rangle_{\text{n}}}{\tau_{\text{E}}}}_{\frac{W}{\tau_{\text{E}}}}  = \underbrace{\frac{f_{\alpha}P_{\alpha} + P_{\text{c}} + P_{\text{OH}} + P_{\text{HCD}}}{V_{\text{P}}} - \frac{P_{\text{rad}}}{V_{\text{p}}}}_{P_{\text{L}}}
+$$
+
+The $\frac{3}{2}n_{\text{i}} \langle T_{\text{i}} \rangle_{\text{n}}$ value is simply the volume averaged ion thermal energy density where $\langle T_{\text{i}} \rangle_{\text{n}}$ is the density weighted temperature. The same goes for the $\frac{3}{2}n_{\text{e}} \langle T_{\text{e}} \rangle_{\text{e}}$ electron thermal energy density term. $\tau_{\text{E}}$ is the confinement time calculated from the chosen confinement scaling via `i_confinement_time`. 
+
+The constraint uses the loss power and thermal densities hence the inclusion of the $V_{\text{p}}$ plasma volume term. The constraint is adapted depending on the condition of `i_rad_loss` which governs the radiation contribution to the loss power definition, see the [radiation and energy confinement section](#effect-of-radiation-on-energy-confinement) for more info. The injected heating and current drive contribution $P_{\text{HCD}}$ is also included or excluded depending if the plasma is deemed to be ignited with the `i_plasma_ignited` switch.
+
+**It is highly recommended to always have this constraint on as it is a global consistency checker**
+
+----------
+
+## Plasma Species Power Balance
+
+While the plasma may be in global power equilibrium, the electrons and ions are continually transferring energy to each other via Coulomb collisions and certain heating and current drive systems only heat particular species. Therefore in order to be in true equilibrium the net power transfer between all species in the plasma must also be equal to zero.
+
+!!! tip "Solving for particle power balance"
+
+    It is highly recommended to have `f_temp_plasma_ion_electron` which represent the ion to electron temperature ratio $\frac{T_i}{T_e}$, as an iteration variable. This allows the solver to have another degree of freedom upon which to change the fusion power and the equilibration power.
+
+-------------------
+
+## Electron Power Balance
+
+This constraint can be activated by stating `icc = 3` in the input file.
+
+For the electrons this internal power blance must be satisifed:
+
+$$
+P_{\Omega} + P_{\text{HCD,e}}+ (f_{e}\times(f_{\alpha}P_{\alpha})) - P_{\text{rad}} - P_{\text{ei}} - \frac{3}{2}\frac{n_{\text{e}} \langle T_{\text{e}} \rangle_{\text{n}}}{\tau_{\text{E}}} = 0
+$$
+
+where $P_{\Omega}$ is the plasma ohmic heating power, $P_{\text{HCD,e}}$ is the external heating and current drive power that only goes to electrons,$(f_{e}\times(f_{\alpha}P_{\alpha}))$ is the fraction of the coupled alpha particle power that goes to the electrons, $P_{\text{rad}}$ is the total radiation power given off by the plasma, $P_{\text{ei}}$ is the [electron-ion equilibration power](#ion-electron-equilibration-power-density--calculate_ion_electron_equilibration_power).
+
+**It is highly recommended to always have this constraint on as it is a plasma equilibrium checker**
+
+------------
+
+## Ion Power Balance
+
+This constraint can be activated by stating `icc = 4` in the input file.
+
+For the ions this internal power blance must be satisifed:
+
+$$
+P_{\text{HCD,i}} + (f_{i}\times(f_{\alpha}P_{\alpha})) + P_{\text{ei}} - \frac{3}{2}\frac{n_{\text{i}} \langle T_{\text{i}} \rangle_{\text{n}}}{\tau_{\text{E}}} = 0
+$$
+
+where $P_{\text{HCD,i}}$ is the external heating and current drive power that only goes to ions,$(f_{e}\times(f_{\alpha}P_{\alpha}))$ is the fraction of the coupled alpha particle power that goes to the ions, $P_{\text{ei}}$ is the [electron-ion equilibration power](#ion-electron-equilibration-power-density--calculate_ion_electron_equilibration_power).
+
+**It is highly recommended to always have this constraint on as it is a plasma equilibrium checker**
+
+-----------
+
+## Ion-electron equilibration power density | `calculate_ion_electron_equilibration_power()`
+
+The equilibration power $P_{\text{ei}}$ refers to the rate of energy transfer between plasma species—typically ions and electrons—due to Coulomb collisions. It is driven by the temperature difference between the species.
+
+$$
+P_{\text{ei}} = \frac{3}{2} \frac{n_e \text{e} (T_i - T_e)}{\tau_{\text{eq}}}
+$$
+
+where $n_e$ is the electron density, $T_i$ and $T_e$ is the ion and electron temperatures in $[\text{eV}]$ and $\tau_{\text{eq}}$ is the ion-electron equilibration time.
\ No newline at end of file

From 987a808d103de3007d1919849aef7a87dfc9f7a0 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 15:43:13 +0100
Subject: [PATCH 10/37] Post rebase fixes

---
 process/core/init.py            | 2 +-
 process/core/io/plot/summary.py | 4 +++-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/process/core/init.py b/process/core/init.py
index 5aba77d345..483719dd9e 100644
--- a/process/core/init.py
+++ b/process/core/init.py
@@ -294,7 +294,7 @@ def check_process(inputs, data):  # noqa: ARG001
         )
 
     # Deprecate constraints
-    for depcrecated_constraint in [3, 4, 10, 74, 42]:
+    for depcrecated_constraint in [10, 74, 42]:
         if (
             data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns]
             == depcrecated_constraint
diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index 57ff833f84..fefe9487cc 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -16173,7 +16173,9 @@ def _add_page(name: str | None = None):
 
     plot_resistivity_profile(pages["detailed_params"].add_subplot(232), m_file, scan)
 
-    plot_plasma_equilibration_time_profile(pages["detailed_params"].add_subplot(233), m_file, scan)
+    plot_plasma_equilibration_time_profile(
+        pages["detailed_params"].add_subplot(233), m_file, scan
+    )
 
     plot_detailed_plasma_parameters(
         pages["detailed_params"].add_subplot(233),

From 5f3e8aa4f7c5b64d7c915070da5394bfb739050d Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 16:19:58 +0100
Subject: [PATCH 11/37] Rename f_alpha_electron to f_p_alpha_total_electron and
 update references in physics models and tests

---
 process/data_structure/physics_variables.py        |  4 ++--
 process/models/physics/physics.py                  | 10 +++++-----
 process/models/stellarator/stellarator.py          |  2 +-
 tests/unit/models/physics/test_fusion_reactions.py |  8 ++++----
 tests/unit/models/physics/test_physics.py          | 10 +++++-----
 5 files changed, 17 insertions(+), 17 deletions(-)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 007b00520d..2ee8ea99a7 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -684,8 +684,8 @@ class PhysicsData:
     f_c_plasma_inductive: float = 0.0
     """fraction of plasma current produced inductively"""
 
-    f_alpha_electron: float = 0.0
-    """fraction of alpha energy to electrons"""
+    f_p_alpha_total_electron: float = 0.0
+    """Fraction of total alpha power used to heat the electrons"""
 
     f_p_alpha_plasma_deposited: float = 0.95
     """Fraction of alpha power deposited in plasma. Default of 0.95 taken from https://doi.org/10.1088/0029-5515/39/12/305."""
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 248e3c7f94..d400e43fed 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -666,7 +666,7 @@ def run(self):
             self.data.physics.p_charged_particle_mw,
             self.data.physics.p_fusion_total_mw,
         ) = reactions.set_fusion_powers(
-            self.data.physics.f_alpha_electron,
+            self.data.physics.f_p_alpha_total_electron,
             self.data.physics.f_alpha_ion,
             self.data.physics.p_beam_alpha_mw,
             self.data.physics.pden_non_alpha_charged_mw,
@@ -1347,10 +1347,10 @@ def plasma_composition(self):
         else:
             pc = self.data.physics.f_temp_plasma_electron_density_vol_avg
 
-        self.data.physics.f_alpha_electron = 0.88155 * np.exp(
+        self.data.physics.f_p_alpha_total_electron = 0.88155 * np.exp(
             -self.data.physics.temp_plasma_electron_vol_avg_kev * pc / 67.4036
         )
-        self.data.physics.f_alpha_ion = 1.0 - self.data.physics.f_alpha_electron
+        self.data.physics.f_alpha_ion = 1.0 - self.data.physics.f_p_alpha_total_electron
 
         # ======================================================================
 
@@ -2283,8 +2283,8 @@ def outplas(self):
         po.ovarre(
             self.outfile,
             "Fraction of alpha power to electrons",
-            "(f_alpha_electron)",
-            self.data.physics.f_alpha_electron,
+            "(f_p_alpha_total_electron)",
+            self.data.physics.f_p_alpha_total_electron,
             "OP ",
         )
         po.ovarre(
diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py
index c5e699c126..89844df9bc 100644
--- a/process/models/stellarator/stellarator.py
+++ b/process/models/stellarator/stellarator.py
@@ -2055,7 +2055,7 @@ def st_phys(self, output):
             self.data.physics.p_charged_particle_mw,
             self.data.physics.p_fusion_total_mw,
         ) = reactions.set_fusion_powers(
-            self.data.physics.f_alpha_electron,
+            self.data.physics.f_p_alpha_total_electron,
             self.data.physics.f_alpha_ion,
             self.data.physics.p_beam_alpha_mw,
             self.data.physics.pden_non_alpha_charged_mw,
diff --git a/tests/unit/models/physics/test_fusion_reactions.py b/tests/unit/models/physics/test_fusion_reactions.py
index b64821e975..5eafb79abc 100644
--- a/tests/unit/models/physics/test_fusion_reactions.py
+++ b/tests/unit/models/physics/test_fusion_reactions.py
@@ -12,7 +12,7 @@
 class SetFusionPowersParam(NamedTuple):
     f_p_alpha_plasma_deposited: Any = None
 
-    f_alpha_electron: Any = None
+    f_p_alpha_total_electron: Any = None
 
     f_alpha_ion: Any = None
 
@@ -50,7 +50,7 @@ class SetFusionPowersParam(NamedTuple):
     [
         SetFusionPowersParam(
             f_p_alpha_plasma_deposited=0.95,
-            f_alpha_electron=0.68,
+            f_p_alpha_total_electron=0.68,
             f_alpha_ion=0.32,
             p_beam_alpha_mw=0,
             pden_non_alpha_charged_mw=0.00066,
@@ -69,7 +69,7 @@ class SetFusionPowersParam(NamedTuple):
         ),
         SetFusionPowersParam(
             f_p_alpha_plasma_deposited=0.95,
-            f_alpha_electron=0.68,
+            f_p_alpha_total_electron=0.68,
             f_alpha_ion=0.32,
             p_beam_alpha_mw=100.5,
             pden_non_alpha_charged_mw=0.00066,
@@ -115,7 +115,7 @@ def test_set_fusion_powers(setfusionpowersparam):
         p_charged_particle_mw,
         p_fusion_total_mw,
     ) = reactions.set_fusion_powers(
-        f_alpha_electron=setfusionpowersparam.f_alpha_electron,
+        f_p_alpha_total_electron=setfusionpowersparam.f_p_alpha_total_electron,
         f_alpha_ion=setfusionpowersparam.f_alpha_ion,
         p_beam_alpha_mw=setfusionpowersparam.p_beam_alpha_mw,
         pden_non_alpha_charged_mw=setfusionpowersparam.pden_non_alpha_charged_mw,
diff --git a/tests/unit/models/physics/test_physics.py b/tests/unit/models/physics/test_physics.py
index 5aa70f4689..13cf59a052 100644
--- a/tests/unit/models/physics/test_physics.py
+++ b/tests/unit/models/physics/test_physics.py
@@ -1235,7 +1235,7 @@ class PlasmaCompositionParam(NamedTuple):
 
     i_plasma_ignited: Any = None
 
-    f_alpha_electron: Any = None
+    f_p_alpha_total_electron: Any = None
 
     m_fuel_amu: Any = None
 
@@ -1379,7 +1379,7 @@ class PlasmaCompositionParam(NamedTuple):
             ]),
             alphat=1.45,
             i_plasma_ignited=PlasmaIgnitionModel.NON_IGNITED,
-            f_alpha_electron=0,
+            f_p_alpha_total_electron=0,
             m_fuel_amu=0,
             f_plasma_fuel_tritium=0.5,
             nd_plasma_fuel_ions_vol_avg=0,
@@ -1485,7 +1485,7 @@ class PlasmaCompositionParam(NamedTuple):
             ).transpose(),
             alphat=1.45,
             i_plasma_ignited=PlasmaIgnitionModel.NON_IGNITED,
-            f_alpha_electron=0.6845930883190634,
+            f_p_alpha_total_electron=0.6845930883190634,
             m_fuel_amu=2.5,
             f_plasma_fuel_tritium=0.5,
             nd_plasma_fuel_ions_vol_avg=5.8589175702454272e19,
@@ -1581,7 +1581,7 @@ def test_plasma_composition(plasmacompositionparam, monkeypatch, physics):
     for field in [
         "alphat",
         "i_plasma_ignited",
-        "f_alpha_electron",
+        "f_p_alpha_total_electron",
         "m_fuel_amu",
         "f_plasma_fuel_tritium",
         "nd_plasma_fuel_ions_vol_avg",
@@ -1621,7 +1621,7 @@ def test_plasma_composition(plasmacompositionparam, monkeypatch, physics):
         plasmacompositionparam.expected_impurity_arr_frac
     )
 
-    assert physics.data.physics.f_alpha_electron == pytest.approx(
+    assert physics.data.physics.f_p_alpha_total_electron == pytest.approx(
         plasmacompositionparam.expected_f_alpha_electron
     )
 

From db558c1bacb5e5a2928e0c8b1394cb2544356dfb Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 16:22:09 +0100
Subject: [PATCH 12/37] Rename f_alpha_ion to f_p_alpha_total_ions and update
 references in fusion reactions and tests

---
 process/data_structure/physics_variables.py        |  4 ++--
 process/models/physics/fusion_reactions.py         |  8 +++++---
 process/models/physics/physics.py                  | 10 ++++++----
 process/models/stellarator/stellarator.py          |  2 +-
 tests/unit/models/physics/test_fusion_reactions.py |  8 ++++----
 tests/unit/models/physics/test_physics.py          | 10 +++++-----
 6 files changed, 23 insertions(+), 19 deletions(-)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 2ee8ea99a7..1db51d8f2b 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -690,8 +690,8 @@ class PhysicsData:
     f_p_alpha_plasma_deposited: float = 0.95
     """Fraction of alpha power deposited in plasma. Default of 0.95 taken from https://doi.org/10.1088/0029-5515/39/12/305."""
 
-    f_alpha_ion: float = 0.0
-    """fraction of alpha power to ions"""
+    f_p_alpha_total_ions: float = 0.0
+    """Fraction of total alpha power used to heat the ions"""
 
     f_plasma_fuel_deuterium: float = 0.5
     """Plasma deuterium fuel fraction"""
diff --git a/process/models/physics/fusion_reactions.py b/process/models/physics/fusion_reactions.py
index 53d44266ff..5201708eec 100644
--- a/process/models/physics/fusion_reactions.py
+++ b/process/models/physics/fusion_reactions.py
@@ -796,7 +796,7 @@ def bosch_hale_reactivity(
 
 def set_fusion_powers(
     f_alpha_electron: float,
-    f_alpha_ion: float,
+    f_p_alpha_total_ions: float,
     p_beam_alpha_mw: float,
     pden_non_alpha_charged_mw: float,
     pden_plasma_neutron_mw: float,
@@ -811,7 +811,7 @@ def set_fusion_powers(
     ----------
     f_alpha_electron :
         float
-    f_alpha_ion :
+    f_p_alpha_total_ions :
         float
     p_beam_alpha_mw :
         float
@@ -896,7 +896,9 @@ def set_fusion_powers(
     # Alpha power to electrons and ions (used with electron
     # and ion power balance equations only)
     # No consideration of pden_non_alpha_charged_mw here.
-    f_pden_alpha_ions_mw = f_p_alpha_plasma_deposited * pden_alpha_total_mw * f_alpha_ion
+    f_pden_alpha_ions_mw = (
+        f_p_alpha_plasma_deposited * pden_alpha_total_mw * f_p_alpha_total_ions
+    )
     f_pden_alpha_electron_mw = (
         f_p_alpha_plasma_deposited * pden_alpha_total_mw * f_alpha_electron
     )
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index d400e43fed..1d1f4729f0 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -667,7 +667,7 @@ def run(self):
             self.data.physics.p_fusion_total_mw,
         ) = reactions.set_fusion_powers(
             self.data.physics.f_p_alpha_total_electron,
-            self.data.physics.f_alpha_ion,
+            self.data.physics.f_p_alpha_total_ions,
             self.data.physics.p_beam_alpha_mw,
             self.data.physics.pden_non_alpha_charged_mw,
             self.data.physics.pden_plasma_neutron_mw,
@@ -1350,7 +1350,9 @@ def plasma_composition(self):
         self.data.physics.f_p_alpha_total_electron = 0.88155 * np.exp(
             -self.data.physics.temp_plasma_electron_vol_avg_kev * pc / 67.4036
         )
-        self.data.physics.f_alpha_ion = 1.0 - self.data.physics.f_p_alpha_total_electron
+        self.data.physics.f_p_alpha_total_ions = (
+            1.0 - self.data.physics.f_p_alpha_total_electron
+        )
 
         # ======================================================================
 
@@ -2290,8 +2292,8 @@ def outplas(self):
         po.ovarre(
             self.outfile,
             "Fraction of alpha power to ions",
-            "(f_alpha_ion)",
-            self.data.physics.f_alpha_ion,
+            "(f_p_alpha_total_ions)",
+            self.data.physics.f_p_alpha_total_ions,
             "OP ",
         )
         po.ovarre(
diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py
index 89844df9bc..4c46e3b8d4 100644
--- a/process/models/stellarator/stellarator.py
+++ b/process/models/stellarator/stellarator.py
@@ -2056,7 +2056,7 @@ def st_phys(self, output):
             self.data.physics.p_fusion_total_mw,
         ) = reactions.set_fusion_powers(
             self.data.physics.f_p_alpha_total_electron,
-            self.data.physics.f_alpha_ion,
+            self.data.physics.f_p_alpha_total_ions,
             self.data.physics.p_beam_alpha_mw,
             self.data.physics.pden_non_alpha_charged_mw,
             self.data.physics.pden_plasma_neutron_mw,
diff --git a/tests/unit/models/physics/test_fusion_reactions.py b/tests/unit/models/physics/test_fusion_reactions.py
index 5eafb79abc..a9d4654314 100644
--- a/tests/unit/models/physics/test_fusion_reactions.py
+++ b/tests/unit/models/physics/test_fusion_reactions.py
@@ -14,7 +14,7 @@ class SetFusionPowersParam(NamedTuple):
 
     f_p_alpha_total_electron: Any = None
 
-    f_alpha_ion: Any = None
+    f_p_alpha_total_ions: Any = None
 
     p_beam_alpha_mw: Any = None
 
@@ -51,7 +51,7 @@ class SetFusionPowersParam(NamedTuple):
         SetFusionPowersParam(
             f_p_alpha_plasma_deposited=0.95,
             f_p_alpha_total_electron=0.68,
-            f_alpha_ion=0.32,
+            f_p_alpha_total_ions=0.32,
             p_beam_alpha_mw=0,
             pden_non_alpha_charged_mw=0.00066,
             vol_plasma=2426.25,
@@ -70,7 +70,7 @@ class SetFusionPowersParam(NamedTuple):
         SetFusionPowersParam(
             f_p_alpha_plasma_deposited=0.95,
             f_p_alpha_total_electron=0.68,
-            f_alpha_ion=0.32,
+            f_p_alpha_total_ions=0.32,
             p_beam_alpha_mw=100.5,
             pden_non_alpha_charged_mw=0.00066,
             vol_plasma=2426.25,
@@ -116,7 +116,7 @@ def test_set_fusion_powers(setfusionpowersparam):
         p_fusion_total_mw,
     ) = reactions.set_fusion_powers(
         f_p_alpha_total_electron=setfusionpowersparam.f_p_alpha_total_electron,
-        f_alpha_ion=setfusionpowersparam.f_alpha_ion,
+        f_p_alpha_total_ions=setfusionpowersparam.f_p_alpha_total_ions,
         p_beam_alpha_mw=setfusionpowersparam.p_beam_alpha_mw,
         pden_non_alpha_charged_mw=setfusionpowersparam.pden_non_alpha_charged_mw,
         pden_plasma_neutron_mw=setfusionpowersparam.pden_plasma_neutron_mw,
diff --git a/tests/unit/models/physics/test_physics.py b/tests/unit/models/physics/test_physics.py
index 13cf59a052..a241a81b56 100644
--- a/tests/unit/models/physics/test_physics.py
+++ b/tests/unit/models/physics/test_physics.py
@@ -1255,7 +1255,7 @@ class PlasmaCompositionParam(NamedTuple):
 
     f_nd_plasma_oxygen_electron: Any = None
 
-    f_alpha_ion: Any = None
+    f_p_alpha_total_ions: Any = None
 
     f_nd_alpha_thermal_electron: Any = None
 
@@ -1389,7 +1389,7 @@ class PlasmaCompositionParam(NamedTuple):
             n_charge_plasma_effective_mass_weighted_vol_avg=0,
             f_nd_plasma_carbon_electron=0,
             f_nd_plasma_oxygen_electron=0,
-            f_alpha_ion=0,
+            f_p_alpha_total_ions=0,
             f_nd_alpha_thermal_electron=0.10000000000000001,
             f_nd_beam_electron=0,
             n_charge_plasma_effective_vol_avg=0,
@@ -1495,7 +1495,7 @@ class PlasmaCompositionParam(NamedTuple):
             n_charge_plasma_effective_mass_weighted_vol_avg=0.43046641789338563,
             f_nd_plasma_carbon_electron=0,
             f_nd_plasma_oxygen_electron=0,
-            f_alpha_ion=0.3154069116809366,
+            f_p_alpha_total_ions=0.3154069116809366,
             f_nd_alpha_thermal_electron=0.10000000000000001,
             f_nd_beam_electron=0,
             n_charge_plasma_effective_vol_avg=2.0909945616489103,
@@ -1591,7 +1591,7 @@ def test_plasma_composition(plasmacompositionparam, monkeypatch, physics):
         "n_charge_plasma_effective_mass_weighted_vol_avg",
         "f_nd_plasma_carbon_electron",
         "f_nd_plasma_oxygen_electron",
-        "f_alpha_ion",
+        "f_p_alpha_total_ions",
         "f_nd_alpha_thermal_electron",
         "f_nd_beam_electron",
         "n_charge_plasma_effective_vol_avg",
@@ -1646,7 +1646,7 @@ def test_plasma_composition(plasmacompositionparam, monkeypatch, physics):
         == pytest.approx(plasmacompositionparam.expected_zeffai)
     )
 
-    assert physics.data.physics.f_alpha_ion == pytest.approx(
+    assert physics.data.physics.f_p_alpha_total_ions == pytest.approx(
         plasmacompositionparam.expected_f_alpha_ion
     )
 

From 65c27564a7b9197e6b2e032b0dbe2bcaee4c9653 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 16:28:52 +0100
Subject: [PATCH 13/37] Rename f_pden_alpha_electron_mw to
 pden_alpha_heating_electrons_mw and update references across the codebase

---
 process/core/io/plot/summary.py                    |  2 +-
 process/core/solver/constraints.py                 |  6 +++---
 process/data_structure/physics_variables.py        | 10 +++++-----
 process/models/physics/fusion_reactions.py         |  6 +++---
 process/models/physics/physics.py                  |  6 +++---
 process/models/stellarator/stellarator.py          |  2 +-
 tests/unit/models/physics/test_fusion_reactions.py |  4 ++--
 7 files changed, 18 insertions(+), 18 deletions(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index fefe9487cc..e6ad1dc488 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -11811,7 +11811,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int):
         f"Rate density, plasma: {mfile.get('fusden_plasma_alpha', scan=scan):.4e} particles/m3/sec\n\n"
         f"Total power density: {mfile.get('pden_alpha_total_mw', scan=scan):.4e} MW/m3\n"
         f"Plasma power density: {mfile.get('pden_plasma_alpha_mw', scan=scan):.4e} MW/m3\n\n"
-        f"Power per unit volume transferred to electrons: {mfile.get('f_pden_alpha_electron_mw', scan=scan):.4e} MW/m3\n"
+        f"Power per unit volume transferred to electrons: {mfile.get('pden_alpha_heating_electrons_mw', scan=scan):.4e} MW/m3\n"
         f"Power per unit volume transferred to ions: {mfile.get('f_pden_alpha_ions_mw', scan=scan):.4e} MW/m3\n\n"
     )
 
diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py
index 3c2406507e..34daf3d09b 100644
--- a/process/core/solver/constraints.py
+++ b/process/core/solver/constraints.py
@@ -385,7 +385,7 @@ def constraint_equation_4(constraint_registration, data):
         total core radiation power per volume (MW/m3)
     f_p_alpha_plasma_deposited:
         fraction of alpha power deposited in plasma
-    f_pden_alpha_electron_mw:
+    pden_alpha_heating_electrons_mw:
         alpha power per volume to electrons (MW/m3)
     pden_ion_electron_equilibration_mw:
         ion/electron equilibration power per volume (MW/m3)
@@ -411,7 +411,7 @@ def constraint_equation_4(constraint_registration, data):
     ):
         pdenom = (
             data.physics.f_p_alpha_plasma_deposited
-            * data.physics.f_pden_alpha_electron_mw
+            * data.physics.pden_alpha_heating_electrons_mw
             + data.physics.pden_ion_electron_equilibration_mw
             + data.current_drive.p_hcd_injected_electrons_mw / data.physics.vol_plasma
         )
@@ -419,7 +419,7 @@ def constraint_equation_4(constraint_registration, data):
         # if plasma ignited
         pdenom = (
             data.physics.f_p_alpha_plasma_deposited
-            * data.physics.f_pden_alpha_electron_mw
+            * data.physics.pden_alpha_heating_electrons_mw
             + data.physics.pden_ion_electron_equilibration_mw
         )
 
diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 1db51d8f2b..e949343b29 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1020,19 +1020,19 @@ class PhysicsData:
     """branching ratio for DD -> T"""
 
     pden_plasma_alpha_mw: float = 0.0
-    """Alpha power per volume just from plasma [MW/m3]"""
+    """Alpha power per volume just from plasma [MW/m³]"""
 
     pden_alpha_total_mw: float = 0.0
-    """Alpha power per volume from plasma and beams [MW/m3]"""
+    """Alpha power per volume from plasma and beams [MW/m³]"""
 
-    f_pden_alpha_electron_mw: float = 0.0
-    """Alpha power per volume to electrons [MW/m3]"""
+    pden_alpha_heating_electrons_mw: float = 0.0
+    """Alpha power heating per unit volume to electrons [MW/m³]"""
 
     p_fw_alpha_mw: float = 0.0
     """alpha power escaping plasma and reaching first wall (MW)"""
 
     f_pden_alpha_ions_mw: float = 0.0
-    """alpha power per volume to ions (MW/m3)"""
+    """alpha power per volume to ions [MW/m³]"""
 
     p_plasma_alpha_mw: float = 0.0
     """Alpha power from only the plasma (MW)"""
diff --git a/process/models/physics/fusion_reactions.py b/process/models/physics/fusion_reactions.py
index 5201708eec..5987fa5c59 100644
--- a/process/models/physics/fusion_reactions.py
+++ b/process/models/physics/fusion_reactions.py
@@ -840,7 +840,7 @@ def set_fusion_powers(
         - p_non_alpha_charged_mw (float): Other total charged particle fusion power [MW].
         - pden_alpha_total_mw (float): Alpha power per unit volume, from beams and
           plasma [MW/m³].
-        - f_pden_alpha_electron_mw (float): Alpha power per unit volume to
+        - pden_alpha_heating_electrons_mw (float): Alpha power per unit volume to
           electrons [MW/m³].
         - f_pden_alpha_ions_mw (float): Alpha power per unit volume to ions [MW/m³].
         - p_charged_particle_mw (float): Charged particle fusion power [MW].
@@ -899,7 +899,7 @@ def set_fusion_powers(
     f_pden_alpha_ions_mw = (
         f_p_alpha_plasma_deposited * pden_alpha_total_mw * f_p_alpha_total_ions
     )
-    f_pden_alpha_electron_mw = (
+    pden_alpha_heating_electrons_mw = (
         f_p_alpha_plasma_deposited * pden_alpha_total_mw * f_alpha_electron
     )
 
@@ -911,7 +911,7 @@ def set_fusion_powers(
         p_neutron_total_mw,
         p_non_alpha_charged_mw,
         pden_alpha_total_mw,
-        f_pden_alpha_electron_mw,
+        pden_alpha_heating_electrons_mw,
         f_pden_alpha_ions_mw,
         p_charged_particle_mw,
         p_fusion_total_mw,
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 1d1f4729f0..1a7986acfe 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -661,7 +661,7 @@ def run(self):
             self.data.physics.p_neutron_total_mw,
             self.data.physics.p_non_alpha_charged_mw,
             self.data.physics.pden_alpha_total_mw,
-            self.data.physics.f_pden_alpha_electron_mw,
+            self.data.physics.pden_alpha_heating_electrons_mw,
             self.data.physics.f_pden_alpha_ions_mw,
             self.data.physics.p_charged_particle_mw,
             self.data.physics.p_fusion_total_mw,
@@ -1994,8 +1994,8 @@ def outplas(self):
         po.ovarre(
             self.outfile,
             "Alpha power per unit volume transferred to electrons (MW/m³)",
-            "(f_pden_alpha_electron_mw)",
-            self.data.physics.f_pden_alpha_electron_mw,
+            "(pden_alpha_heating_electrons_mw)",
+            self.data.physics.pden_alpha_heating_electrons_mw,
             "OP ",
         )
         po.ovarre(
diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py
index 4c46e3b8d4..c865426546 100644
--- a/process/models/stellarator/stellarator.py
+++ b/process/models/stellarator/stellarator.py
@@ -2050,7 +2050,7 @@ def st_phys(self, output):
             self.data.physics.p_neutron_total_mw,
             self.data.physics.p_non_alpha_charged_mw,
             self.data.physics.pden_alpha_total_mw,
-            self.data.physics.f_pden_alpha_electron_mw,
+            self.data.physics.pden_alpha_heating_electrons_mw,
             self.data.physics.f_pden_alpha_ions_mw,
             self.data.physics.p_charged_particle_mw,
             self.data.physics.p_fusion_total_mw,
diff --git a/tests/unit/models/physics/test_fusion_reactions.py b/tests/unit/models/physics/test_fusion_reactions.py
index a9d4654314..8b318ad819 100644
--- a/tests/unit/models/physics/test_fusion_reactions.py
+++ b/tests/unit/models/physics/test_fusion_reactions.py
@@ -110,7 +110,7 @@ def test_set_fusion_powers(setfusionpowersparam):
         p_neutron_total_mw,
         p_non_alpha_charged_mw,
         pden_alpha_total_mw,
-        f_pden_alpha_electron_mw,
+        pden_alpha_heating_electrons_mw,
         f_pden_alpha_ions_mw,
         p_charged_particle_mw,
         p_fusion_total_mw,
@@ -140,7 +140,7 @@ def test_set_fusion_powers(setfusionpowersparam):
     assert p_non_alpha_charged_mw == pytest.approx(
         setfusionpowersparam.expected_non_alpha_charged_power
     )
-    assert f_pden_alpha_electron_mw == pytest.approx(
+    assert pden_alpha_heating_electrons_mw == pytest.approx(
         setfusionpowersparam.expected_alpha_power_electron_density
     )
     assert f_pden_alpha_ions_mw == pytest.approx(

From aab5e395a3907009a9b93fa998f78d383e0a951d Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 16:30:05 +0100
Subject: [PATCH 14/37] Rename f_pden_alpha_ions_mw to
 pden_alpha_heating_ions_mw and update references across the codebase

---
 process/core/io/plot/summary.py                    | 2 +-
 process/core/solver/constraints.py                 | 9 ++++++---
 process/data_structure/physics_variables.py        | 4 ++--
 process/models/physics/fusion_reactions.py         | 6 +++---
 process/models/physics/physics.py                  | 6 +++---
 process/models/stellarator/stellarator.py          | 2 +-
 tests/unit/models/physics/test_fusion_reactions.py | 4 ++--
 7 files changed, 18 insertions(+), 15 deletions(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index e6ad1dc488..9be2d668e6 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -11812,7 +11812,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int):
         f"Total power density: {mfile.get('pden_alpha_total_mw', scan=scan):.4e} MW/m3\n"
         f"Plasma power density: {mfile.get('pden_plasma_alpha_mw', scan=scan):.4e} MW/m3\n\n"
         f"Power per unit volume transferred to electrons: {mfile.get('pden_alpha_heating_electrons_mw', scan=scan):.4e} MW/m3\n"
-        f"Power per unit volume transferred to ions: {mfile.get('f_pden_alpha_ions_mw', scan=scan):.4e} MW/m3\n\n"
+        f"Power per unit volume transferred to ions: {mfile.get('pden_alpha_heating_ions_mw', scan=scan):.4e} MW/m3\n\n"
     )
 
     axis.text(
diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py
index 34daf3d09b..da2f13dea3 100644
--- a/process/core/solver/constraints.py
+++ b/process/core/solver/constraints.py
@@ -326,7 +326,7 @@ def constraint_equation_3(constraint_registration, data):
         ion/electron equilibration power per volume (MW/m3)
     f_p_alpha_plasma_deposited:
         fraction of alpha power deposited in plasma
-    f_pden_alpha_ions_mw:
+    pden_alpha_heating_ions_mw:
         alpha power per volume to ions (MW/m3)
     p_hcd_injected_ions_mw:
         auxiliary injected power to ions (MW)
@@ -345,7 +345,7 @@ def constraint_equation_3(constraint_registration, data):
             ),
             (
                 data.physics.f_p_alpha_plasma_deposited
-                * data.physics.f_pden_alpha_ions_mw
+                * data.physics.pden_alpha_heating_ions_mw
                 + data.current_drive.p_hcd_injected_ions_mw / data.physics.vol_plasma
             ),
             constraint_registration,
@@ -357,7 +357,10 @@ def constraint_equation_3(constraint_registration, data):
             data.physics.pden_ion_transport_loss_mw
             + data.physics.pden_ion_electron_equilibration_mw
         ),
-        (data.physics.f_p_alpha_plasma_deposited * data.physics.f_pden_alpha_ions_mw),
+        (
+            data.physics.f_p_alpha_plasma_deposited
+            * data.physics.pden_alpha_heating_ions_mw
+        ),
         constraint_registration,
     )
 
diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index e949343b29..7cc02d5518 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1031,8 +1031,8 @@ class PhysicsData:
     p_fw_alpha_mw: float = 0.0
     """alpha power escaping plasma and reaching first wall (MW)"""
 
-    f_pden_alpha_ions_mw: float = 0.0
-    """alpha power per volume to ions [MW/m³]"""
+    pden_alpha_heating_ions_mw: float = 0.0
+    """Alpha power heating per unit volume to ions [MW/m³]"""
 
     p_plasma_alpha_mw: float = 0.0
     """Alpha power from only the plasma (MW)"""
diff --git a/process/models/physics/fusion_reactions.py b/process/models/physics/fusion_reactions.py
index 5987fa5c59..995ec44a59 100644
--- a/process/models/physics/fusion_reactions.py
+++ b/process/models/physics/fusion_reactions.py
@@ -842,7 +842,7 @@ def set_fusion_powers(
           plasma [MW/m³].
         - pden_alpha_heating_electrons_mw (float): Alpha power per unit volume to
           electrons [MW/m³].
-        - f_pden_alpha_ions_mw (float): Alpha power per unit volume to ions [MW/m³].
+        - pden_alpha_heating_ions_mw (float): Alpha power per unit volume to ions [MW/m³].
         - p_charged_particle_mw (float): Charged particle fusion power [MW].
         - p_fusion_total_mw (float): Total fusion power [MW].
 
@@ -896,7 +896,7 @@ def set_fusion_powers(
     # Alpha power to electrons and ions (used with electron
     # and ion power balance equations only)
     # No consideration of pden_non_alpha_charged_mw here.
-    f_pden_alpha_ions_mw = (
+    pden_alpha_heating_ions_mw = (
         f_p_alpha_plasma_deposited * pden_alpha_total_mw * f_p_alpha_total_ions
     )
     pden_alpha_heating_electrons_mw = (
@@ -912,7 +912,7 @@ def set_fusion_powers(
         p_non_alpha_charged_mw,
         pden_alpha_total_mw,
         pden_alpha_heating_electrons_mw,
-        f_pden_alpha_ions_mw,
+        pden_alpha_heating_ions_mw,
         p_charged_particle_mw,
         p_fusion_total_mw,
     )
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 1a7986acfe..c1aa4b36a7 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -662,7 +662,7 @@ def run(self):
             self.data.physics.p_non_alpha_charged_mw,
             self.data.physics.pden_alpha_total_mw,
             self.data.physics.pden_alpha_heating_electrons_mw,
-            self.data.physics.f_pden_alpha_ions_mw,
+            self.data.physics.pden_alpha_heating_ions_mw,
             self.data.physics.p_charged_particle_mw,
             self.data.physics.p_fusion_total_mw,
         ) = reactions.set_fusion_powers(
@@ -2001,8 +2001,8 @@ def outplas(self):
         po.ovarre(
             self.outfile,
             "Alpha power per unit volume transferred to ions (MW/m³)",
-            "(f_pden_alpha_ions_mw)",
-            self.data.physics.f_pden_alpha_ions_mw,
+            "(pden_alpha_heating_ions_mw)",
+            self.data.physics.pden_alpha_heating_ions_mw,
             "OP ",
         )
 
diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py
index c865426546..15d674f88a 100644
--- a/process/models/stellarator/stellarator.py
+++ b/process/models/stellarator/stellarator.py
@@ -2051,7 +2051,7 @@ def st_phys(self, output):
             self.data.physics.p_non_alpha_charged_mw,
             self.data.physics.pden_alpha_total_mw,
             self.data.physics.pden_alpha_heating_electrons_mw,
-            self.data.physics.f_pden_alpha_ions_mw,
+            self.data.physics.pden_alpha_heating_ions_mw,
             self.data.physics.p_charged_particle_mw,
             self.data.physics.p_fusion_total_mw,
         ) = reactions.set_fusion_powers(
diff --git a/tests/unit/models/physics/test_fusion_reactions.py b/tests/unit/models/physics/test_fusion_reactions.py
index 8b318ad819..9bbc1f1a3d 100644
--- a/tests/unit/models/physics/test_fusion_reactions.py
+++ b/tests/unit/models/physics/test_fusion_reactions.py
@@ -111,7 +111,7 @@ def test_set_fusion_powers(setfusionpowersparam):
         p_non_alpha_charged_mw,
         pden_alpha_total_mw,
         pden_alpha_heating_electrons_mw,
-        f_pden_alpha_ions_mw,
+        pden_alpha_heating_ions_mw,
         p_charged_particle_mw,
         p_fusion_total_mw,
     ) = reactions.set_fusion_powers(
@@ -143,7 +143,7 @@ def test_set_fusion_powers(setfusionpowersparam):
     assert pden_alpha_heating_electrons_mw == pytest.approx(
         setfusionpowersparam.expected_alpha_power_electron_density
     )
-    assert f_pden_alpha_ions_mw == pytest.approx(
+    assert pden_alpha_heating_ions_mw == pytest.approx(
         setfusionpowersparam.expected_alpha_power_ion_density
     )
     assert p_charged_particle_mw == pytest.approx(

From 90d44493076f7a639378f860320868fd579caf78 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 16:42:33 +0100
Subject: [PATCH 15/37] Refactor constraint_equation_3 to simplify ion heating
 and loss calculations

---
 process/core/solver/constraints.py | 34 +++++++++++++++---------------
 1 file changed, 17 insertions(+), 17 deletions(-)

diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py
index da2f13dea3..d2bc14c588 100644
--- a/process/core/solver/constraints.py
+++ b/process/core/solver/constraints.py
@@ -338,29 +338,29 @@ def constraint_equation_3(constraint_registration, data):
         PlasmaIgnitionModel(data.physics.i_plasma_ignited)
         == PlasmaIgnitionModel.NON_IGNITED
     ):
+        p_ion_heating = (
+            data.physics.pden_alpha_heating_ions_mw
+            + (data.current_drive.p_hcd_injected_ions_mw / data.physics.vol_plasma)
+            + data.physics.pden_ion_electron_equilibration_mw
+        )
+
+        p_ion_loss = data.physics.pden_ion_transport_loss_mw
         return eq(
-            (
-                data.physics.pden_ion_transport_loss_mw
-                + data.physics.pden_ion_electron_equilibration_mw
-            ),
-            (
-                data.physics.f_p_alpha_plasma_deposited
-                * data.physics.pden_alpha_heating_ions_mw
-                + data.current_drive.p_hcd_injected_ions_mw / data.physics.vol_plasma
-            ),
+            (p_ion_loss),
+            (p_ion_heating),
             constraint_registration,
         )
 
     # Plasma ignited
+    p_ion_heating = (
+        data.physics.pden_alpha_heating_ions_mw
+        + data.physics.pden_ion_electron_equilibration_mw
+    )
+
+    p_ion_loss = data.physics.pden_ion_transport_loss_mw
     return eq(
-        (
-            data.physics.pden_ion_transport_loss_mw
-            + data.physics.pden_ion_electron_equilibration_mw
-        ),
-        (
-            data.physics.f_p_alpha_plasma_deposited
-            * data.physics.pden_alpha_heating_ions_mw
-        ),
+        (p_ion_loss),
+        (p_ion_heating),
         constraint_registration,
     )
 

From 8467b1deffc11f1003fecf83d4130e0f88076f22 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 17:23:32 +0100
Subject: [PATCH 16/37] Update constraint_equation_3 to include non-alpha
 charged particle contributions to heating and loss calculations

---
 process/core/solver/constraints.py | 31 +++++++++++++++---------------
 1 file changed, 16 insertions(+), 15 deletions(-)

diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py
index d2bc14c588..ea6f2e710e 100644
--- a/process/core/solver/constraints.py
+++ b/process/core/solver/constraints.py
@@ -342,6 +342,7 @@ def constraint_equation_3(constraint_registration, data):
             data.physics.pden_alpha_heating_ions_mw
             + (data.current_drive.p_hcd_injected_ions_mw / data.physics.vol_plasma)
             + data.physics.pden_ion_electron_equilibration_mw
+            + data.physics.pden_non_alpha_charged_mw
         )
 
         p_ion_loss = data.physics.pden_ion_transport_loss_mw
@@ -355,6 +356,7 @@ def constraint_equation_3(constraint_registration, data):
     p_ion_heating = (
         data.physics.pden_alpha_heating_ions_mw
         + data.physics.pden_ion_electron_equilibration_mw
+        + data.physics.pden_non_alpha_charged_mw
     )
 
     p_ion_loss = data.physics.pden_ion_transport_loss_mw
@@ -398,35 +400,34 @@ def constraint_equation_4(constraint_registration, data):
     """
     # pscaling: total transport power per volume (MW/m3)
 
-    pscaling = data.physics.pden_electron_transport_loss_mw
+    p_electron_loss = (
+        data.physics.pden_electron_transport_loss_mw
+        + data.physics.pden_ion_electron_equilibration_mw
+    )
     # Total power lost is scaling power plus radiation:
     if data.physics.i_rad_loss == 0:
-        pnumerator = pscaling + data.physics.pden_plasma_rad_mw
+        p_electron_loss += data.physics.pden_plasma_rad_mw
     elif data.physics.i_rad_loss == 1:
-        pnumerator = pscaling + data.physics.pden_plasma_core_rad_mw
-    else:
-        pnumerator = pscaling
+        p_electron_loss += data.physics.pden_plasma_core_rad_mw
 
     # if plasma not ignited include injected power
     if (
         PlasmaIgnitionModel(data.physics.i_plasma_ignited)
         == PlasmaIgnitionModel.NON_IGNITED
     ):
-        pdenom = (
-            data.physics.f_p_alpha_plasma_deposited
-            * data.physics.pden_alpha_heating_electrons_mw
-            + data.physics.pden_ion_electron_equilibration_mw
-            + data.current_drive.p_hcd_injected_electrons_mw / data.physics.vol_plasma
+        p_electron_heating = (
+            data.physics.pden_alpha_heating_electrons_mw
+            + (data.current_drive.p_hcd_injected_electrons_mw / data.physics.vol_plasma)
+            + data.physics.pden_plasma_ohmic_mw
         )
     else:
         # if plasma ignited
-        pdenom = (
-            data.physics.f_p_alpha_plasma_deposited
-            * data.physics.pden_alpha_heating_electrons_mw
-            + data.physics.pden_ion_electron_equilibration_mw
+        p_electron_heating = (
+            data.physics.pden_alpha_heating_electrons_mw
+            + data.physics.pden_plasma_ohmic_mw
         )
 
-    return eq(pnumerator, pdenom, constraint_registration)
+    return eq(p_electron_loss, p_electron_heating, constraint_registration)
 
 
 @ConstraintManager.register_constraint(5, "/m³", "<=")

From 95a62f60e5b9190906e6b18749e82137e208dd2d Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 17:25:29 +0100
Subject: [PATCH 17/37] Rename pden_ion_electron_equilibration_mw to
 pden_ion_electron_equilibration_vol_avg_mw and update references across the
 codebase

---
 process/core/solver/constraints.py          | 10 +++++-----
 process/data_structure/physics_variables.py |  4 ++--
 process/models/physics/physics.py           |  6 +++---
 process/models/stellarator/stellarator.py   |  2 +-
 4 files changed, 11 insertions(+), 11 deletions(-)

diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py
index ea6f2e710e..4e8f50393c 100644
--- a/process/core/solver/constraints.py
+++ b/process/core/solver/constraints.py
@@ -322,7 +322,7 @@ def constraint_equation_3(constraint_registration, data):
 
     pden_ion_transport_loss_mw:
         ion transport power per volume (MW/m3)
-    pden_ion_electron_equilibration_mw:
+    pden_ion_electron_equilibration_vol_avg_mw:
         ion/electron equilibration power per volume (MW/m3)
     f_p_alpha_plasma_deposited:
         fraction of alpha power deposited in plasma
@@ -341,7 +341,7 @@ def constraint_equation_3(constraint_registration, data):
         p_ion_heating = (
             data.physics.pden_alpha_heating_ions_mw
             + (data.current_drive.p_hcd_injected_ions_mw / data.physics.vol_plasma)
-            + data.physics.pden_ion_electron_equilibration_mw
+            + data.physics.pden_ion_electron_equilibration_vol_avg_mw
             + data.physics.pden_non_alpha_charged_mw
         )
 
@@ -355,7 +355,7 @@ def constraint_equation_3(constraint_registration, data):
     # Plasma ignited
     p_ion_heating = (
         data.physics.pden_alpha_heating_ions_mw
-        + data.physics.pden_ion_electron_equilibration_mw
+        + data.physics.pden_ion_electron_equilibration_vol_avg_mw
         + data.physics.pden_non_alpha_charged_mw
     )
 
@@ -392,7 +392,7 @@ def constraint_equation_4(constraint_registration, data):
         fraction of alpha power deposited in plasma
     pden_alpha_heating_electrons_mw:
         alpha power per volume to electrons (MW/m3)
-    pden_ion_electron_equilibration_mw:
+    pden_ion_electron_equilibration_vol_avg_mw:
         ion/electron equilibration power per volume (MW/m3)
     p_hcd_injected_electrons_mw:
         auxiliary injected power to electrons (MW)
@@ -402,7 +402,7 @@ def constraint_equation_4(constraint_registration, data):
 
     p_electron_loss = (
         data.physics.pden_electron_transport_loss_mw
-        + data.physics.pden_ion_electron_equilibration_mw
+        + data.physics.pden_ion_electron_equilibration_vol_avg_mw
     )
     # Total power lost is scaling power plus radiation:
     if data.physics.i_rad_loss == 0:
diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 7cc02d5518..933b02de94 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1109,8 +1109,8 @@ class PhysicsData:
     pflux_fw_rad_mw: float = 0.0
     """Nominal mean radiation load on inside surface of reactor (MW/m2)"""
 
-    pden_ion_electron_equilibration_mw: float = 0.0
-    """ion/electron equilibration power per volume (MW/m3)"""
+    pden_ion_electron_equilibration_vol_avg_mw: float = 0.0
+    """Ion/electron volume averaged equilibration power per unit volume [MW/m³]"""
 
     plasma_current: float = 0.0
     """Plasma current (Iₚ) [A]"""
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index c1aa4b36a7..8bf960b541 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -692,7 +692,7 @@ def run(self):
 
         # Calculate ion/electron equilibration power
 
-        pden_ion_electron_equilibration = calculate_ion_electron_equilibration_power(
+        pden_ion_electron_equilibration_vol_avg = calculate_ion_electron_equilibration_power(
             nd_electrons=self.data.physics.nd_plasma_electrons_vol_avg,
             temp_plasma_electron_kev=self.data.physics.temp_plasma_electron_vol_avg_kev,
             temp_plasma_ion_kev=self.data.physics.temp_plasma_ion_vol_avg_kev,
@@ -701,8 +701,8 @@ def run(self):
                 self.data.physics.t_plasma_electron_triton_equilibration_vol_avg,
             ]),
         )
-        self.data.physics.pden_ion_electron_equilibration_mw = (
-            pden_ion_electron_equilibration / 1.0e6
+        self.data.physics.pden_ion_electron_equilibration_vol_avg_mw = (
+            pden_ion_electron_equilibration_vol_avg / 1.0e6
         )
 
         # Calculate radiation power
diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py
index 15d674f88a..b99325ed68 100644
--- a/process/models/stellarator/stellarator.py
+++ b/process/models/stellarator/stellarator.py
@@ -2107,7 +2107,7 @@ def st_phys(self, output):
 
         # Calculate ion/electron equilibration power
 
-        self.data.physics.pden_ion_electron_equilibration_mw = (
+        self.data.physics.pden_ion_electron_equilibration_vol_avg_mw = (
             calculate_ion_electron_equilibration_power(
                 self.data.physics.alphan,
                 self.data.physics.alphat,

From 34e2ee984cdd4c1b6fa9a988c4060aab3db00eaf Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 18:36:04 +0100
Subject: [PATCH 18/37] Refactor plasma power balance output by consolidating
 heating and transport power calculations into a new method

---
 process/models/physics/physics.py | 163 ++++++++++++++++++------------
 1 file changed, 96 insertions(+), 67 deletions(-)

diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 8bf960b541..3ef5ddcb11 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -2261,73 +2261,10 @@ def outplas(self):
                 )
 
         po.oblnkl(self.outfile)
-        po.ovarre(
-            self.outfile,
-            "Total heating power given to the plasma (Pₕₑₐₜ) [MW]",
-            "(p_plasma_heating_total_mw)",
-            self.data.physics.p_plasma_heating_total_mw,
-            "OP ",
-        )
-        po.ovarre(
-            self.outfile,
-            "Ohmic heating power (MW)",
-            "(p_plasma_ohmic_mw)",
-            self.data.physics.p_plasma_ohmic_mw,
-            "OP ",
-        )
-        po.ovarre(
-            self.outfile,
-            "Fraction of alpha power deposited in plasma",
-            "(f_p_alpha_plasma_deposited)",
-            self.data.physics.f_p_alpha_plasma_deposited,
-            "IP",
-        )
-        po.ovarre(
-            self.outfile,
-            "Fraction of alpha power to electrons",
-            "(f_p_alpha_total_electron)",
-            self.data.physics.f_p_alpha_total_electron,
-            "OP ",
-        )
-        po.ovarre(
-            self.outfile,
-            "Fraction of alpha power to ions",
-            "(f_p_alpha_total_ions)",
-            self.data.physics.f_p_alpha_total_ions,
-            "OP ",
-        )
-        po.ovarre(
-            self.outfile,
-            "Ion transport (MW)",
-            "(p_ion_transport_loss_mw)",
-            self.data.physics.p_ion_transport_loss_mw,
-            "OP ",
-        )
-        po.ovarre(
-            self.outfile,
-            "Electron transport (MW)",
-            "(p_electron_transport_loss_mw)",
-            self.data.physics.p_electron_transport_loss_mw,
-            "OP ",
-        )
-        po.ovarre(
-            self.outfile,
-            "Injection power to ions (MW)",
-            "(p_hcd_injected_ions_mw)",
-            self.data.current_drive.p_hcd_injected_ions_mw,
-            "OP ",
-        )
-        po.ovarre(
-            self.outfile,
-            "Injection power to electrons (MW)",
-            "(p_hcd_injected_electrons_mw)",
-            self.data.current_drive.p_hcd_injected_electrons_mw,
-            "OP ",
-        )
-        if (
-            PlasmaIgnitionModel(self.data.physics.i_plasma_ignited)
-            == PlasmaIgnitionModel.IGNITED
-        ):
+
+        self.output_plasma_power_balance()
+
+        if self.data.physics.i_plasma_ignited == 1:
             po.ocmmnt(self.outfile, "  (Injected power only used for start-up phase)")
 
         self.exhaust.output()
@@ -2431,6 +2368,98 @@ def outplas(self):
                 self.data.reinke.fzactual,
             )
 
+    def output_plasma_power_balance(self) -> None:
+        """Output information about plasma power balance."""
+        po.oheadr(self.outfile, "Plasma Power Balance")
+
+        po.osubhd(self.outfile, "Global power balance :")
+
+        po.ovarre(
+            self.outfile,
+            "Ohmic heating power (MW)",
+            "(p_plasma_ohmic_mw)",
+            self.data.physics.p_plasma_ohmic_mw,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "Total alpha heating power (MW)",
+            "(p_alpha_total_mw)",
+            self.data.physics.p_alpha_total_mw,
+            "OP ",
+        )
+        po.ovarrf(
+            self.outfile,
+            "Fraction of alpha power deposited in plasma",
+            "(f_p_alpha_plasma_deposited)",
+            self.data.physics.f_p_alpha_plasma_deposited,
+            "IP",
+        )
+        po.ovarrf(
+            self.outfile,
+            "Total power injected into the plasma (MW)",
+            "(p_hcd_injected_total_mw)",
+            self.data.current_drive.p_hcd_injected_total_mw,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "Non-alpha heating power (MW)",
+            "(p_non_alpha_charged_mw)",
+            self.data.physics.p_non_alpha_charged_mw,
+            "OP ",
+        )
+        po.oblnkl(self.outfile)
+        po.ovarre(
+            self.outfile,
+            "Total heating power given to the plasma (Pₕₑₐₜ) [MW]",
+            "(p_plasma_heating_total_mw)",
+            self.data.physics.p_plasma_heating_total_mw,
+            "OP ",
+        )
+        po.oblnkl(self.outfile)
+
+        po.osubhd(self.outfile, "Electron power balance :")
+
+        po.ovarre(
+            self.outfile,
+            "Electron transport (MW)",
+            "(p_electron_transport_loss_mw)",
+            self.data.physics.p_electron_transport_loss_mw,
+            "OP ",
+        )
+
+        po.ovarre(
+            self.outfile,
+            "Injection power to electrons (MW)",
+            "(p_hcd_injected_electrons_mw)",
+            self.data.current_drive.p_hcd_injected_electrons_mw,
+            "OP ",
+        )
+        po.osubhd(self.outfile, "Ions power balance :")
+
+        po.ovarrf(
+            self.outfile,
+            "Fraction of alpha power to ions",
+            "(f_p_alpha_total_ions)",
+            self.data.physics.f_p_alpha_total_ions,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "Ion transport (MW)",
+            "(p_ion_transport_loss_mw)",
+            self.data.physics.p_ion_transport_loss_mw,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "Injection power to ions (MW)",
+            "(p_hcd_injected_ions_mw)",
+            self.data.current_drive.p_hcd_injected_ions_mw,
+            "OP ",
+        )
+
     def output_temperature_density_profile_info(self) -> None:
         """Output information about plasma temperature and density profiles."""
         po.oheadr(self.outfile, "Plasma Density and Temperature Profiles")

From d273f479dbfd63fbe9d7930a6335a6d882bfa920 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 19:04:19 +0100
Subject: [PATCH 19/37] Update ion-electron equilibration temperature
 references to use density-weighted averages

---
 process/models/physics/physics.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 3ef5ddcb11..6b4a545f90 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -694,8 +694,8 @@ def run(self):
 
         pden_ion_electron_equilibration_vol_avg = calculate_ion_electron_equilibration_power(
             nd_electrons=self.data.physics.nd_plasma_electrons_vol_avg,
-            temp_plasma_electron_kev=self.data.physics.temp_plasma_electron_vol_avg_kev,
-            temp_plasma_ion_kev=self.data.physics.temp_plasma_ion_vol_avg_kev,
+            temp_plasma_electron_kev=self.data.physics.temp_plasma_electron_density_weighted_kev,
+            temp_plasma_ion_kev=self.data.physics.temp_plasma_ion_density_weighted_kev,
             t_plasma_ion_electron_equilibration=np.average([
                 self.data.physics.t_plasma_electron_deuteron_equilibration_vol_avg,
                 self.data.physics.t_plasma_electron_triton_equilibration_vol_avg,

From 2055861c8ac4a29d1931651330271ba3234d5151 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 19:56:48 +0100
Subject: [PATCH 20/37] Add ion/electron volume averaged equilibration power
 variable and update calculations

---
 process/data_structure/physics_variables.py |  3 +++
 process/models/physics/physics.py           | 11 +++++++++++
 2 files changed, 14 insertions(+)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 933b02de94..3ff7da0b17 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1111,6 +1111,9 @@ class PhysicsData:
 
     pden_ion_electron_equilibration_vol_avg_mw: float = 0.0
     """Ion/electron volume averaged equilibration power per unit volume [MW/m³]"""
+    
+    p_ion_electron_equilibration_vol_avg_mw: float = 0.0
+    """Ion/electron volume averaged equilibration power [MW]"""
 
     plasma_current: float = 0.0
     """Plasma current (Iₚ) [A]"""
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 6b4a545f90..e5681502ac 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -704,6 +704,10 @@ def run(self):
         self.data.physics.pden_ion_electron_equilibration_vol_avg_mw = (
             pden_ion_electron_equilibration_vol_avg / 1.0e6
         )
+        
+        self.data.physics.p_ion_electron_equilibration_vol_avg_mw = (
+            self.data.physics.pden_ion_electron_equilibration_vol_avg_mw * self.data.physics.vol_plasma
+        )
 
         # Calculate radiation power
 
@@ -2436,6 +2440,13 @@ def output_plasma_power_balance(self) -> None:
             self.data.current_drive.p_hcd_injected_electrons_mw,
             "OP ",
         )
+        po.ovarre(
+            self.outfile,
+            "Ion-electron equilibration power (MW)",
+            "(p_ion_electron_equilibration_vol_avg_mw)",
+            self.data.physics.p_ion_electron_equilibration_vol_avg_mw,
+            "OP ",
+        )
         po.osubhd(self.outfile, "Ions power balance :")
 
         po.ovarrf(

From 6ef28deee2743cefbe5a218112c710f6c8affc5c Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 21:07:52 +0100
Subject: [PATCH 21/37] Refactor ion/electron equilibration variables and
 update calculations for clarity

---
 process/data_structure/physics_variables.py |  2 +-
 process/models/physics/confinement_time.py  |  9 ++++---
 process/models/physics/physics.py           | 29 ++++++++++++++++-----
 3 files changed, 28 insertions(+), 12 deletions(-)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 3ff7da0b17..ff964b4cbf 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1111,7 +1111,7 @@ class PhysicsData:
 
     pden_ion_electron_equilibration_vol_avg_mw: float = 0.0
     """Ion/electron volume averaged equilibration power per unit volume [MW/m³]"""
-    
+
     p_ion_electron_equilibration_vol_avg_mw: float = 0.0
     """Ion/electron volume averaged equilibration power [MW]"""
 
diff --git a/process/models/physics/confinement_time.py b/process/models/physics/confinement_time.py
index 2e9dd2ae8a..1fb83d080c 100644
--- a/process/models/physics/confinement_time.py
+++ b/process/models/physics/confinement_time.py
@@ -987,13 +987,14 @@ def calculate_confinement_time(
             eden_plasma_electrons_thermal_vol_avg / t_electron_energy_confinement
         ) / 1e6  # Convert from W/m³ to MW/m³
 
-        ratio = eden_plasma_ions_thermal_vol_avg / eden_plasma_electrons_thermal_vol_avg
+        # ratio = (nd_plasma_ions_total_vol_avg / nd_plasma_electrons_vol_avg) * (
+        #     temp_plasma_ion_density_weighted_kev
+        #     / temp_plasma_electron_density_weighted_kev
+        # )
 
         # Global energy confinement time
 
-        t_energy_confinement = (ratio + 1.0e0) / (
-            ratio / t_ion_energy_confinement + 1.0e0 / t_electron_energy_confinement
-        )
+        t_energy_confinement = t_ion_energy_confinement
 
         # For comparison directly calculate the confinement time from the stored energy
         # calculated from the total plasma beta and the loss power used above.
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index e5681502ac..fa9002bd08 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -704,9 +704,10 @@ def run(self):
         self.data.physics.pden_ion_electron_equilibration_vol_avg_mw = (
             pden_ion_electron_equilibration_vol_avg / 1.0e6
         )
-        
+
         self.data.physics.p_ion_electron_equilibration_vol_avg_mw = (
-            self.data.physics.pden_ion_electron_equilibration_vol_avg_mw * self.data.physics.vol_plasma
+            self.data.physics.pden_ion_electron_equilibration_vol_avg_mw
+            * self.data.physics.vol_plasma
         )
 
         # Calculate radiation power
@@ -5531,6 +5532,8 @@ def run(self):
 
         self.data.physics.t_plasma_electron_deuteron_equilibration_vol_avg = self.calculate_equilibriation_time(
             temp_plasma_electron_kev=self.data.physics.temp_plasma_electron_vol_avg_kev,
+            temp_plasma_ion_kev=self.data.physics.temp_plasma_electron_vol_avg_kev
+            * self.data.physics.f_temp_plasma_ion_electron,
             nd_plasma_ions=self.data.physics.nd_plasma_electrons_vol_avg
             * self.data.physics.f_plasma_fuel_deuterium
             * (
@@ -5544,6 +5547,8 @@ def run(self):
 
         self.data.physics.t_plasma_electron_deuteron_equilibration_profile = self.calculate_equilibriation_time(
             temp_plasma_electron_kev=self.plasma_profile.teprofile.profile_y,
+            temp_plasma_ion_kev=self.plasma_profile.teprofile.profile_y
+            * self.data.physics.f_temp_plasma_ion_electron,
             nd_plasma_ions=(
                 self.plasma_profile.neprofile.profile_y
                 * self.data.physics.f_plasma_fuel_deuterium
@@ -5559,6 +5564,8 @@ def run(self):
 
         self.data.physics.t_plasma_electron_triton_equilibration_vol_avg = self.calculate_equilibriation_time(
             temp_plasma_electron_kev=self.data.physics.temp_plasma_electron_vol_avg_kev,
+            temp_plasma_ion_kev=self.data.physics.temp_plasma_electron_vol_avg_kev
+            * self.data.physics.f_temp_plasma_ion_electron,
             nd_plasma_ions=self.data.physics.nd_plasma_electrons_vol_avg
             * self.data.physics.f_plasma_fuel_tritium
             * (
@@ -5572,6 +5579,8 @@ def run(self):
 
         self.data.physics.t_plasma_electron_triton_equilibration_profile = self.calculate_equilibriation_time(
             temp_plasma_electron_kev=self.plasma_profile.teprofile.profile_y,
+            temp_plasma_ion_kev=self.plasma_profile.teprofile.profile_y
+            * self.data.physics.f_temp_plasma_ion_electron,
             nd_plasma_ions=(
                 self.plasma_profile.neprofile.profile_y
                 * self.data.physics.f_plasma_fuel_tritium
@@ -6031,6 +6040,7 @@ def calculate_spitzer_resistivity(
     @nb.njit(cache=True)
     def calculate_equilibriation_time(
         temp_plasma_electron_kev: float | np.ndarray,
+        temp_plasma_ion_kev: float | np.ndarray,
         nd_plasma_ions: float | np.ndarray,
         plasma_coulomb_log_electron_ion: float | np.ndarray,
         m_ion: float,
@@ -6043,6 +6053,8 @@ def calculate_equilibriation_time(
         ----------
         temp_plasma_electron_kev : float | np.ndarray
             Electron temperature in keV.
+        temp_plasma_ion_kev : float | np.ndarray
+            Ion temperature in keV.
         nd_plasma_ions : float | np.ndarray
             Ion density (m^-3).
         plasma_coulomb_log_electron_ion : float | np.ndarray
@@ -6066,17 +6078,20 @@ def calculate_equilibriation_time(
         """
         return (
             3
-            * (2 * np.pi) ** 1.5
-            * constants.EPSILON0**2
+            * (4 * np.pi * constants.EPSILON0) ** 2.0
             * m_ion
             * (constants.ELECTRON_MASS)
             * (
-                (temp_plasma_electron_kev * constants.KILOELECTRON_VOLT)
-                / constants.ELECTRON_MASS
+                (
+                    (temp_plasma_electron_kev * constants.KILOELECTRON_VOLT)
+                    / constants.ELECTRON_MASS
+                )
+                + ((temp_plasma_ion_kev * constants.KILOELECTRON_VOLT) / m_ion)
             )
             ** 1.5
         ) / (
-            nd_plasma_ions
+            (8 * np.sqrt(2 * np.pi))
+            * nd_plasma_ions
             * n_charge_ion**2
             * plasma_coulomb_log_electron_ion
             * constants.ELECTRON_CHARGE**4

From dca1cf91022731e75ad1136bd8faa4f011ad947b Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 22:53:28 +0100
Subject: [PATCH 22/37] Add electron-alpha equilibration time calculations and
 profiles

---
 process/core/io/plot/summary.py             | 16 ++++++++
 process/data_structure/physics_variables.py |  8 ++++
 process/models/physics/physics.py           | 43 +++++++++++++++++++++
 3 files changed, 67 insertions(+)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index 9be2d668e6..de399365a6 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -14108,6 +14108,13 @@ def plot_plasma_equilibration_time_profile(
         for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)))
     ]
 
+    t_plasma_electron_alpha_equilibration_profile = [
+        mfile_data.data[f"t_plasma_electron_alpha_equilibration_profile{i}"].get_scan(
+            scan
+        )
+        for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)))
+    ]
+
     axis.plot(
         np.linspace(0, 1, len(t_plasma_electron_deuteron_equilibration_profile)),
         t_plasma_electron_deuteron_equilibration_profile,
@@ -14124,6 +14131,14 @@ def plot_plasma_equilibration_time_profile(
         label=r"$\tau_{e-T,eq}$",
     )
 
+    axis.plot(
+        np.linspace(0, 1, len(t_plasma_electron_alpha_equilibration_profile)),
+        t_plasma_electron_alpha_equilibration_profile,
+        color="red",
+        linestyle="-",
+        label=r"$\tau_{e-\alpha,eq}$",
+    )
+
     axis.set_yscale("log")
     axis.set_ylabel("Equilibration Time [s]")
     axis.set_xlabel("$\\rho \\ [r/a]$")
@@ -14831,6 +14846,7 @@ def plot_detailed_plasma_parameters(axis: plt.Axes, fig, mfile: MFile, scan: int
         f"$\\mathbf{{Equilibration \\ Times:}}$\n\n"
         f"$\\langle\\tau_{{e-D,Equil}}\\rangle$: {mfile.get('t_plasma_electron_deuteron_equilibration_vol_avg', scan=scan):.4e} s\n"
         f"$\\langle\\tau_{{e-T,Equil}}\\rangle$: {mfile.get('t_plasma_electron_triton_equilibration_vol_avg', scan=scan):.4e} s\n"
+        f"$\\langle\\tau_{{e-\\alpha,Equil}}\\rangle$: {mfile.get('t_plasma_electron_alpha_equilibration_vol_avg', scan=scan):.4e} s"
     )
 
     light_yellow_box = {
diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index ff964b4cbf..327a706748 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1600,6 +1600,14 @@ class PhysicsData:
     )
     """Profile of electron-triton equilibration time in plasma (s)"""
 
+    t_plasma_electron_alpha_equilibration_vol_avg: float = 4.5
+    """Volume averaged electron-alpha equilibration time in plasma (s)"""
+
+    t_plasma_electron_alpha_equilibration_profile: list[float] = field(
+        default_factory=list
+    )
+    """Profile of electron-alpha equilibration time in plasma (s)"""
+
     freq_plasma_electron_electron_collision_profile: list[float] = field(
         default_factory=list
     )
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index fa9002bd08..3a9f2a9911 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -5594,6 +5594,32 @@ def run(self):
             n_charge_ion=1,
         )
 
+        self.data.physics.t_plasma_electron_alpha_equilibration_vol_avg = self.calculate_equilibriation_time(
+            temp_plasma_electron_kev=self.data.physics.temp_plasma_electron_vol_avg_kev,
+            temp_plasma_ion_kev=self.data.physics.temp_plasma_electron_vol_avg_kev
+            * self.data.physics.f_temp_plasma_ion_electron,
+            nd_plasma_ions=self.data.physics.nd_plasma_alphas_vol_avg,
+            plasma_coulomb_log_electron_ion=self.data.physics.plasma_coulomb_log_electron_alpha_thermal_vol_avg,
+            m_ion=constants.ALPHA_MASS,
+            n_charge_ion=2,
+        )
+
+        self.data.physics.t_plasma_electron_alpha_equilibration_profile = self.calculate_equilibriation_time(
+            temp_plasma_electron_kev=self.plasma_profile.teprofile.profile_y,
+            temp_plasma_ion_kev=self.plasma_profile.teprofile.profile_y
+            * self.data.physics.f_temp_plasma_ion_electron,
+            nd_plasma_ions=(
+                self.plasma_profile.neprofile.profile_y
+                * (
+                    self.data.physics.nd_plasma_alphas_vol_avg
+                    / self.data.physics.nd_plasma_electrons_vol_avg
+                )
+            ),
+            plasma_coulomb_log_electron_ion=self.data.physics.plasma_coulomb_log_electron_alpha_thermal_profile,
+            m_ion=constants.ALPHA_MASS,
+            n_charge_ion=2,
+        )
+
     @staticmethod
     @nb.njit(cache=True)
     def calculate_debye_length(
@@ -6691,3 +6717,20 @@ def output_detailed_physics(self):
                 f"(t_plasma_electron_triton_equilibration_profile{i})",
                 self.data.physics.t_plasma_electron_triton_equilibration_profile[i],
             )
+
+        po.ovarre(
+            self.outfile,
+            "Volume averaged electron-alpha equilibration time (τ_eq) (s)",
+            "(t_plasma_electron_alpha_equilibration_vol_avg)",
+            self.data.physics.t_plasma_electron_alpha_equilibration_vol_avg,
+        )
+
+        for i in range(
+            len(self.data.physics.t_plasma_electron_alpha_equilibration_profile)
+        ):
+            po.ovarre(
+                self.mfile,
+                f"Electron-alpha equilibration time at point {i}",
+                f"(t_plasma_electron_alpha_equilibration_profile{i})",
+                self.data.physics.t_plasma_electron_alpha_equilibration_profile[i],
+            )

From b5f2df08d81d4120a5b6a585aae00d42b8c52811 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 24 Jun 2026 22:58:10 +0100
Subject: [PATCH 23/37] Update ion-electron equilibration time calculation to
 include alpha particle contributions

---
 process/models/physics/physics.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 3a9f2a9911..7633b578e6 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -696,10 +696,12 @@ def run(self):
             nd_electrons=self.data.physics.nd_plasma_electrons_vol_avg,
             temp_plasma_electron_kev=self.data.physics.temp_plasma_electron_density_weighted_kev,
             temp_plasma_ion_kev=self.data.physics.temp_plasma_ion_density_weighted_kev,
-            t_plasma_ion_electron_equilibration=np.average([
-                self.data.physics.t_plasma_electron_deuteron_equilibration_vol_avg,
-                self.data.physics.t_plasma_electron_triton_equilibration_vol_avg,
-            ]),
+            t_plasma_ion_electron_equilibration=1.0
+            / (
+                1.0 / self.data.physics.t_plasma_electron_deuteron_equilibration_vol_avg
+                + 1.0 / self.data.physics.t_plasma_electron_triton_equilibration_vol_avg
+                + 1.0 / self.data.physics.t_plasma_electron_alpha_equilibration_vol_avg
+            ),
         )
         self.data.physics.pden_ion_electron_equilibration_vol_avg_mw = (
             pden_ion_electron_equilibration_vol_avg / 1.0e6

From 23f283968c7059f50344ccdd3bbb4f9b1817a60b Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 29 Jun 2026 14:04:59 +0100
Subject: [PATCH 24/37] Add fuel ion to electron energy confinement time ratio
 variable

---
 process/data_structure/physics_variables.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 327a706748..0cb1eb101f 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1351,6 +1351,9 @@ class PhysicsData:
 
     f_t_alpha_energy_confinement: float = 0.0
     """Alpha particle to energy confinement time ratio"""
+    
+    f_t_fuel_ion_electron_energy_confinement: float = 0.0
+    """Fuel ion to electron energy confinement time ratio"""
 
     temp_plasma_electron_vol_avg_kev: float = 12.9
     """Plasma volume averaged electron temperature (⟨Tₑ⟩) [keV] (`iteration variable 4`)"""

From 20f4065b1ddea59ca9118357564390a44c79fd0a Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 29 Jun 2026 14:44:53 +0100
Subject: [PATCH 25/37] Add fuel ion to electron energy confinement time ratio
 variable and update related calculations

---
 process/core/input.py                       |  3 ++
 process/core/solver/iteration_variables.py  |  3 ++
 process/data_structure/physics_variables.py |  8 ++--
 process/models/physics/confinement_time.py  | 51 ++++++++++++++++-----
 process/models/physics/physics.py           |  8 ++++
 5 files changed, 57 insertions(+), 16 deletions(-)

diff --git a/process/core/input.py b/process/core/input.py
index 1d7c132b6a..270dd3ba3a 100644
--- a/process/core/input.py
+++ b/process/core/input.py
@@ -1184,6 +1184,9 @@ def bounds(self) -> tuple[NumberType | None, NumberType | None]:
     "force_vmcon_inequality_tolerance": InputVariable(
         "numerics", float, range=(0.0, 1e10)
     ),
+    "f_t_fuel_ion_electron_energy_confinement": InputVariable(
+        "physics", float, range=(0.1, 5.0)
+    ),
 }
 
 
diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py
index 61c07979c3..5ab1998a7f 100644
--- a/process/core/solver/iteration_variables.py
+++ b/process/core/solver/iteration_variables.py
@@ -236,6 +236,9 @@ class IterationVariable:
     176: IterationVariable("f_st_coil_aspect", "stellarator", 0.70, 1.30),
     177: IterationVariable("f_a_tf_turn_cable_space_extra_void", "tfcoil", 0.01, 1.0),
     178: IterationVariable("f_temp_plasma_ion_electron", "physics", 0.0, 2.0),
+    179: IterationVariable(
+        "f_t_fuel_ion_electron_energy_confinement", "physics", 0.1, 5.0
+    ),
 }
 
 
diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 0cb1eb101f..1ddd777de9 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1109,10 +1109,10 @@ class PhysicsData:
     pflux_fw_rad_mw: float = 0.0
     """Nominal mean radiation load on inside surface of reactor (MW/m2)"""
 
-    pden_ion_electron_equilibration_vol_avg_mw: float = 0.0
+    pden_ion_electron_equilibration_vol_avg_mw: float = 0.5
     """Ion/electron volume averaged equilibration power per unit volume [MW/m³]"""
 
-    p_ion_electron_equilibration_vol_avg_mw: float = 0.0
+    p_ion_electron_equilibration_vol_avg_mw: float = 80.0
     """Ion/electron volume averaged equilibration power [MW]"""
 
     plasma_current: float = 0.0
@@ -1351,8 +1351,8 @@ class PhysicsData:
 
     f_t_alpha_energy_confinement: float = 0.0
     """Alpha particle to energy confinement time ratio"""
-    
-    f_t_fuel_ion_electron_energy_confinement: float = 0.0
+
+    f_t_fuel_ion_electron_energy_confinement: float = 1.0
     """Fuel ion to electron energy confinement time ratio"""
 
     temp_plasma_electron_vol_avg_kev: float = 12.9
diff --git a/process/models/physics/confinement_time.py b/process/models/physics/confinement_time.py
index 1fb83d080c..3255a7aef8 100644
--- a/process/models/physics/confinement_time.py
+++ b/process/models/physics/confinement_time.py
@@ -931,12 +931,6 @@ def calculate_confinement_time(
                 i_confinement_time=i_confinement_time,
             )
 
-        # Apply H-factor correction to chosen scaling
-        t_electron_energy_confinement = hfact * t_electron_confinement
-
-        # Ion energy confinement time
-        t_ion_energy_confinement = t_electron_energy_confinement
-
         # Calculate H* non-radiation corrected H factor
         # Note: we will assume the IPB-98y2 scaling.
         if (
@@ -978,14 +972,47 @@ def calculate_confinement_time(
         # temp_plasma_electron_density_weighted_kev are in keV, and
         # pden_electron_transport_loss_mw and pden_ion_transport_loss_mw are in MW/m3)
 
+        e_ions_stored_thermal = (
+            (3 / 2)
+            * (constants.ELECTRON_CHARGE / 1e3)
+            * nd_plasma_ions_total_vol_avg
+            * temp_plasma_ion_density_weighted_kev
+        )
+
+        e_electrons_stored_thermal = (
+            (3 / 2)
+            * (constants.ELECTRON_CHARGE / 1e3)
+            * nd_plasma_electrons_vol_avg
+            * temp_plasma_electron_density_weighted_kev
+        )
+        
+        t_energy_confinement = t_electron_confinement
+
+        # Apply H-factor correction to chosen scaling
+        t_electron_energy_confinement = (
+            hfact
+            * t_energy_confinement
+            * (
+                (
+                    e_electrons_stored_thermal
+                    + (
+                        e_ions_stored_thermal
+                        / self.data.physics.f_t_fuel_ion_electron_energy_confinement
+                    )
+                )
+                / (e_electrons_stored_thermal + e_ions_stored_thermal)
+            )
+        )
+
+        # Ion energy confinement time
+        t_ion_energy_confinement = t_electron_energy_confinement*self.data.physics.f_t_fuel_ion_electron_energy_confinement
+
         # The transport losses is just the electron and ion thermal energies divided by
         # the confinement time.
-        pden_ion_transport_loss_mw = (
-            eden_plasma_ions_thermal_vol_avg / t_ion_energy_confinement
-        ) / 1e6  # Convert from W/m³ to MW/m³
+        pden_ion_transport_loss_mw = e_ions_stored_thermal / t_ion_energy_confinement
         pden_electron_transport_loss_mw = (
-            eden_plasma_electrons_thermal_vol_avg / t_electron_energy_confinement
-        ) / 1e6  # Convert from W/m³ to MW/m³
+            e_electrons_stored_thermal / t_electron_energy_confinement
+        )
 
         # ratio = (nd_plasma_ions_total_vol_avg / nd_plasma_electrons_vol_avg) * (
         #     temp_plasma_ion_density_weighted_kev
@@ -994,7 +1021,7 @@ def calculate_confinement_time(
 
         # Global energy confinement time
 
-        t_energy_confinement = t_ion_energy_confinement
+        #t_energy_confinement = t_ion_energy_confinement
 
         # For comparison directly calculate the confinement time from the stored energy
         # calculated from the total plasma beta and the loss power used above.
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 7633b578e6..d52061e3f6 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -711,6 +711,7 @@ def run(self):
             self.data.physics.pden_ion_electron_equilibration_vol_avg_mw
             * self.data.physics.vol_plasma
         )
+        print(self.data.physics.f_t_fuel_ion_electron_energy_confinement)
 
         # Calculate radiation power
 
@@ -2450,6 +2451,13 @@ def output_plasma_power_balance(self) -> None:
             self.data.physics.p_ion_electron_equilibration_vol_avg_mw,
             "OP ",
         )
+        po.ovarre(
+            self.outfile,
+            "Ratio of fuel ion to electron energy confinement times",
+            "(f_t_fuel_ion_electron_energy_confinement)",
+            self.data.physics.f_t_fuel_ion_electron_energy_confinement,
+            "OP ",
+        )
         po.osubhd(self.outfile, "Ions power balance :")
 
         po.ovarrf(

From cd819c5a375b09a45920b539bae64ebd0c939136 Mon Sep 17 00:00:00 2001
From: chris-ashe <christopher.ashe@ukaea.uk>
Date: Mon, 29 Jun 2026 22:53:20 +0100
Subject: [PATCH 26/37] Fix bug where hfact wasnt applied to global energy
 confinement time

---
 process/models/physics/confinement_time.py | 11 +++++++----
 process/models/physics/physics.py          |  1 -
 2 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/process/models/physics/confinement_time.py b/process/models/physics/confinement_time.py
index 3255a7aef8..288a2cc941 100644
--- a/process/models/physics/confinement_time.py
+++ b/process/models/physics/confinement_time.py
@@ -985,8 +985,8 @@ def calculate_confinement_time(
             * nd_plasma_electrons_vol_avg
             * temp_plasma_electron_density_weighted_kev
         )
-        
-        t_energy_confinement = t_electron_confinement
+
+        t_energy_confinement = t_electron_confinement * hfact
 
         # Apply H-factor correction to chosen scaling
         t_electron_energy_confinement = (
@@ -1005,7 +1005,10 @@ def calculate_confinement_time(
         )
 
         # Ion energy confinement time
-        t_ion_energy_confinement = t_electron_energy_confinement*self.data.physics.f_t_fuel_ion_electron_energy_confinement
+        t_ion_energy_confinement = (
+            t_electron_energy_confinement
+            * self.data.physics.f_t_fuel_ion_electron_energy_confinement
+        )
 
         # The transport losses is just the electron and ion thermal energies divided by
         # the confinement time.
@@ -1021,7 +1024,7 @@ def calculate_confinement_time(
 
         # Global energy confinement time
 
-        #t_energy_confinement = t_ion_energy_confinement
+        # t_energy_confinement = t_ion_energy_confinement
 
         # For comparison directly calculate the confinement time from the stored energy
         # calculated from the total plasma beta and the loss power used above.
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index d52061e3f6..5c681b30f2 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -711,7 +711,6 @@ def run(self):
             self.data.physics.pden_ion_electron_equilibration_vol_avg_mw
             * self.data.physics.vol_plasma
         )
-        print(self.data.physics.f_t_fuel_ion_electron_energy_confinement)
 
         # Calculate radiation power
 

From c9ed40b43e1155d3b86f7827ee9237c14dc75a06 Mon Sep 17 00:00:00 2001
From: chris-ashe <christopher.ashe@ukaea.uk>
Date: Mon, 29 Jun 2026 23:35:13 +0100
Subject: [PATCH 27/37] Add correct thermal energy calculation

---
 process/core/solver/constraints.py | 18 +++++++++++++++---
 1 file changed, 15 insertions(+), 3 deletions(-)

diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py
index 4e8f50393c..7ba7a9b310 100644
--- a/process/core/solver/constraints.py
+++ b/process/core/solver/constraints.py
@@ -279,10 +279,22 @@ def constraint_equation_2(constraint_registration, data):
     """
     # pscaling: total transport power per volume (MW/m3)
 
-    pscaling = (
-        data.physics.pden_electron_transport_loss_mw
-        + data.physics.pden_ion_transport_loss_mw
+    e_ions_stored_thermal = (
+        (3 / 2)
+        * (constants.ELECTRON_CHARGE / 1e3)
+        * data.physics.nd_plasma_ions_total_vol_avg
+        * data.physics.temp_plasma_ion_density_weighted_kev
     )
+
+    e_electrons_stored_thermal = (
+        (3 / 2)
+        * (constants.ELECTRON_CHARGE / 1e3)
+        * data.physics.nd_plasma_electrons_vol_avg
+        * data.physics.temp_plasma_electron_density_weighted_kev
+    )
+    pscaling = (
+        e_ions_stored_thermal + e_electrons_stored_thermal
+    ) / data.physics.t_energy_confinement
     # Total power lost is scaling power plus radiation:
     if data.physics.i_rad_loss == 0:
         pnumerator = pscaling + data.physics.pden_plasma_rad_mw

From de473c49e30a7c3c71fb4c0d53ddf32e3ca96f48 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Wed, 1 Jul 2026 09:49:54 +0100
Subject: [PATCH 28/37] Update docs to add new equations and fix equation
 formatting

---
 .../source/physics-models/detailed_physics.md          | 10 ++++++++++
 .../source/physics-models/plasma_power_balance.md      | 10 +++++-----
 2 files changed, 15 insertions(+), 5 deletions(-)

diff --git a/documentation/source/physics-models/detailed_physics.md b/documentation/source/physics-models/detailed_physics.md
index bc1844a899..b928feb8f5 100644
--- a/documentation/source/physics-models/detailed_physics.md
+++ b/documentation/source/physics-models/detailed_physics.md
@@ -160,3 +160,13 @@ For $T_\text{e}$ in eV
 $$
 \eta_{\text{spitzer}} = \frac{4 \sqrt{2 \pi}}{3 }\frac{Z_i e^2 \sqrt{m_{\text{e}}} \ln \Lambda_{\text{ei}}}{\left(4 \pi \epsilon_0 \right)^2 T_{\text{e}}^{\frac{3}{2}}}
 $$
+
+------------------
+
+### Ion-electron equilibration time | `calculate_equilibriation_time()`
+
+For species temperatures $T$, measured in $[\text{eV}]$
+
+$$
+\frac{3(4\pi\epsilon_0)^2m_{\text{i}} m_{\text{e}}}{8\sqrt{2\pi}n_{\text{i}} \text{e}^4 Z_{\text{i}}^2 \ln\Lambda_{\text{ei}}}\left(\frac{T_{\text{e}}}{m_{\text{e}}}+\frac{T_{\text{i}}}{m_{\text{i}}}\right)^\frac{3}{2}
+$$
diff --git a/documentation/source/physics-models/plasma_power_balance.md b/documentation/source/physics-models/plasma_power_balance.md
index 791640435a..0725625a3b 100644
--- a/documentation/source/physics-models/plasma_power_balance.md
+++ b/documentation/source/physics-models/plasma_power_balance.md
@@ -18,7 +18,7 @@ P_{\text{L}} = \frac{W}{\tau_{\text{E}}}
 $$
 
 $$
-\underbrace{\frac{3}{2}\frac{n_{\text{i}} \langle T_{\text{i}} \rangle_{\text{n}}}{\tau_{\text{E}}} + \frac{3}{2}\frac{n_{\text{e}} \langle T_{\text{e}} \rangle_{\text{n}}}{\tau_{\text{E}}}}_{\frac{W}{\tau_{\text{E}}}}  = \underbrace{\frac{f_{\alpha}P_{\alpha} + P_{\text{c}} + P_{\text{OH}} + P_{\text{HCD}}}{V_{\text{P}}} - \frac{P_{\text{rad}}}{V_{\text{p}}}}_{P_{\text{L}}}
+\underbrace{\frac{3}{2}\frac{n_{\text{i}} \langle T_{\text{i}} \rangle_{\text{n}}}{\tau_{\text{i}}} + \frac{3}{2}\frac{n_{\text{e}} \langle T_{\text{e}} \rangle_{\text{n}}}{\tau_{\text{e}}}}_{\frac{W}{\tau_{\text{E}}}}  = \underbrace{\frac{f_{\alpha}P_{\alpha} + P_{\text{c}} + P_{\text{OH}} + P_{\text{HCD}}}{V_{\text{P}}} - \frac{P_{\text{rad}}}{V_{\text{p}}}}_{P_{\text{L}}}
 $$
 
 The $\frac{3}{2}n_{\text{i}} \langle T_{\text{i}} \rangle_{\text{n}}$ value is simply the volume averaged ion thermal energy density where $\langle T_{\text{i}} \rangle_{\text{n}}$ is the density weighted temperature. The same goes for the $\frac{3}{2}n_{\text{e}} \langle T_{\text{e}} \rangle_{\text{e}}$ electron thermal energy density term. $\tau_{\text{E}}$ is the confinement time calculated from the chosen confinement scaling via `i_confinement_time`. 
@@ -46,7 +46,7 @@ This constraint can be activated by stating `icc = 3` in the input file.
 For the electrons this internal power blance must be satisifed:
 
 $$
-P_{\Omega} + P_{\text{HCD,e}}+ (f_{e}\times(f_{\alpha}P_{\alpha})) - P_{\text{rad}} - P_{\text{ei}} - \frac{3}{2}\frac{n_{\text{e}} \langle T_{\text{e}} \rangle_{\text{n}}}{\tau_{\text{E}}} = 0
+\underbrace{P_{\Omega} + P_{\text{HCD,e}}+ (f_{e}\times(f_{\alpha}P_{\alpha}))}_{\text{Power gain}} = \underbrace{P_{\text{rad}} - P_{\text{ei}} - \frac{3}{2}\frac{n_{\text{e}} \langle T_{\text{e}} \rangle_{\text{n}}}{\tau_{\text{e}}}}_{\text{Power loss}}
 $$
 
 where $P_{\Omega}$ is the plasma ohmic heating power, $P_{\text{HCD,e}}$ is the external heating and current drive power that only goes to electrons,$(f_{e}\times(f_{\alpha}P_{\alpha}))$ is the fraction of the coupled alpha particle power that goes to the electrons, $P_{\text{rad}}$ is the total radiation power given off by the plasma, $P_{\text{ei}}$ is the [electron-ion equilibration power](#ion-electron-equilibration-power-density--calculate_ion_electron_equilibration_power).
@@ -62,7 +62,7 @@ This constraint can be activated by stating `icc = 4` in the input file.
 For the ions this internal power blance must be satisifed:
 
 $$
-P_{\text{HCD,i}} + (f_{i}\times(f_{\alpha}P_{\alpha})) + P_{\text{ei}} - \frac{3}{2}\frac{n_{\text{i}} \langle T_{\text{i}} \rangle_{\text{n}}}{\tau_{\text{E}}} = 0
+\underbrace{P_{\text{HCD,i}} + (f_{i}\times(f_{\alpha}P_{\alpha})) + P_{\text{ei}}}_{\text{Power gain}} = \underbrace{\frac{3}{2}\frac{n_{\text{i}} \langle T_{\text{i}} \rangle_{\text{n}}}{\tau_{\text{i}}}}_{\text{Power Loss}}
 $$
 
 where $P_{\text{HCD,i}}$ is the external heating and current drive power that only goes to ions,$(f_{e}\times(f_{\alpha}P_{\alpha}))$ is the fraction of the coupled alpha particle power that goes to the ions, $P_{\text{ei}}$ is the [electron-ion equilibration power](#ion-electron-equilibration-power-density--calculate_ion_electron_equilibration_power).
@@ -73,10 +73,10 @@ where $P_{\text{HCD,i}}$ is the external heating and current drive power that on
 
 ## Ion-electron equilibration power density | `calculate_ion_electron_equilibration_power()`
 
-The equilibration power $P_{\text{ei}}$ refers to the rate of energy transfer between plasma species—typically ions and electrons—due to Coulomb collisions. It is driven by the temperature difference between the species.
+The equilibration power $P_{\text{ei}}$ refers to the rate of energy transfer between plasma species—typically ions and electrons—due to elastic Coulomb collisions. It is driven by the temperature difference between the species. The equation below represents the rate of energy transfer from the electrons to the ions due to these elastic Coulomb collisions caused by the temperature difference between the electrons and the ions:
 
 $$
-P_{\text{ei}} = \frac{3}{2} \frac{n_e \text{e} (T_i - T_e)}{\tau_{\text{eq}}}
+P_{\text{ei}} = \frac{3}{2} \frac{n_e \text{e} (T_e - T_i)}{\tau_{\text{eq}}}
 $$
 
 where $n_e$ is the electron density, $T_i$ and $T_e$ is the ion and electron temperatures in $[\text{eV}]$ and $\tau_{\text{eq}}$ is the ion-electron equilibration time.
\ No newline at end of file

From 28d4181671c8e5d52ab22e7c30280190dd86dc22 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Thu, 16 Jul 2026 12:09:29 +0100
Subject: [PATCH 29/37] Post rebase fixes

---
 process/models/physics/confinement_time.py | 31 ++++++++--------------
 process/models/physics/physics.py          |  4 +--
 2 files changed, 13 insertions(+), 22 deletions(-)

diff --git a/process/models/physics/confinement_time.py b/process/models/physics/confinement_time.py
index 288a2cc941..6b031def36 100644
--- a/process/models/physics/confinement_time.py
+++ b/process/models/physics/confinement_time.py
@@ -972,20 +972,6 @@ def calculate_confinement_time(
         # temp_plasma_electron_density_weighted_kev are in keV, and
         # pden_electron_transport_loss_mw and pden_ion_transport_loss_mw are in MW/m3)
 
-        e_ions_stored_thermal = (
-            (3 / 2)
-            * (constants.ELECTRON_CHARGE / 1e3)
-            * nd_plasma_ions_total_vol_avg
-            * temp_plasma_ion_density_weighted_kev
-        )
-
-        e_electrons_stored_thermal = (
-            (3 / 2)
-            * (constants.ELECTRON_CHARGE / 1e3)
-            * nd_plasma_electrons_vol_avg
-            * temp_plasma_electron_density_weighted_kev
-        )
-
         t_energy_confinement = t_electron_confinement * hfact
 
         # Apply H-factor correction to chosen scaling
@@ -994,13 +980,16 @@ def calculate_confinement_time(
             * t_energy_confinement
             * (
                 (
-                    e_electrons_stored_thermal
+                    eden_plasma_electrons_thermal_vol_avg
                     + (
-                        e_ions_stored_thermal
+                        eden_plasma_ions_thermal_vol_avg
                         / self.data.physics.f_t_fuel_ion_electron_energy_confinement
                     )
                 )
-                / (e_electrons_stored_thermal + e_ions_stored_thermal)
+                / (
+                    eden_plasma_electrons_thermal_vol_avg
+                    + eden_plasma_ions_thermal_vol_avg
+                )
             )
         )
 
@@ -1012,10 +1001,12 @@ def calculate_confinement_time(
 
         # The transport losses is just the electron and ion thermal energies divided by
         # the confinement time.
-        pden_ion_transport_loss_mw = e_ions_stored_thermal / t_ion_energy_confinement
+        pden_ion_transport_loss_mw = (
+            eden_plasma_ions_thermal_vol_avg / t_ion_energy_confinement
+        ) / 1e6  # Convert from W/m³ to MW/m³
         pden_electron_transport_loss_mw = (
-            e_electrons_stored_thermal / t_electron_energy_confinement
-        )
+            eden_plasma_electrons_thermal_vol_avg / t_electron_energy_confinement
+        ) / 1e6  # Convert from W/m³ to MW/m³
 
         # ratio = (nd_plasma_ions_total_vol_avg / nd_plasma_electrons_vol_avg) * (
         #     temp_plasma_ion_density_weighted_kev
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 5c681b30f2..59580f05f4 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -5607,7 +5607,7 @@ def run(self):
             temp_plasma_electron_kev=self.data.physics.temp_plasma_electron_vol_avg_kev,
             temp_plasma_ion_kev=self.data.physics.temp_plasma_electron_vol_avg_kev
             * self.data.physics.f_temp_plasma_ion_electron,
-            nd_plasma_ions=self.data.physics.nd_plasma_alphas_vol_avg,
+            nd_plasma_ions=self.data.physics.nd_plasma_alphas_thermal_vol_avg,
             plasma_coulomb_log_electron_ion=self.data.physics.plasma_coulomb_log_electron_alpha_thermal_vol_avg,
             m_ion=constants.ALPHA_MASS,
             n_charge_ion=2,
@@ -5620,7 +5620,7 @@ def run(self):
             nd_plasma_ions=(
                 self.plasma_profile.neprofile.profile_y
                 * (
-                    self.data.physics.nd_plasma_alphas_vol_avg
+                    self.data.physics.nd_plasma_alphas_thermal_vol_avg
                     / self.data.physics.nd_plasma_electrons_vol_avg
                 )
             ),

From ab1295ac13e46daa2c613849f5c69593c32273db Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Thu, 16 Jul 2026 13:57:22 +0100
Subject: [PATCH 30/37] Add function for global and species energy confinement
 time species consitency and add derivation to documentation

---
 .../physics-models/plasma_confinement.md      | 83 ++++++++++++++-----
 .../physics-models/plasma_power_balance.md    |  4 -
 process/models/physics/confinement_time.py    | 77 ++++++++++-------
 3 files changed, 111 insertions(+), 53 deletions(-)

diff --git a/documentation/source/physics-models/plasma_confinement.md b/documentation/source/physics-models/plasma_confinement.md
index 80124f1f89..24277c770d 100644
--- a/documentation/source/physics-models/plasma_confinement.md
+++ b/documentation/source/physics-models/plasma_confinement.md
@@ -647,60 +647,105 @@ $$
 
 -------------------------
 
+### Species Confinement Times
+
+Once the global energy confinement time $\tau_{\text{E}}$ is calculated from the scaling we calculate the species and $H$ factor corrected confinement time.
+
+Using out $H$ factor we find the new corrected value for the global energy confinement time:
+
+$$
+\tau_{\text{E}} = H\tau_{\text{E,scaling}}
+$$
+
+We then find the electron confinement time via the consistency relation from [`calculate_electron_species_consistent_energy_confinement_time()`](#global-confinement-time-consistency--calculate_electron_species_consistent_energy_confinement_time)
+
+$$
+\tau_{\text{e}} = \tau_{\text{E}}\left(\frac{W_{\text{e}} +\frac{W_{\text{i}}}{R}}{W_{\text{e}} +W_{\text{i}}}\right)
+$$
+
+The ion energy confinement time is then simply defined as:
+
+$$
+\tau_{\text{i}} = R\tau_{\text{e}}
+$$
+
+
+----------------
+
 ### Transport Powers
 
-After the confinement time scaling with $H$-factor correction has been calculated, the ion and electron transport power densities are found. `PROCESS` assumes the scaling confinement time to be equal to the ion and electron energy confinement time.
+After the confinement time scaling with $H$-factor correction has been calculated along with the specific ion and electron energy confinement times, the ion and electron transport power densities are found.
 
 This is simply the volume averaged thermal energy of the electron and ions divided by the $H$-factor corrected confinement time from the chosen scaling.
 
 $$
-\mathtt{pden\_ion\_transport\_loss\_mw} = \frac{3}{2}\frac{n_{\text{i}} \langle T_{\text{i}} \rangle_{\text{n}}}{\tau_{\text{E}}}
+\overbrace{\left(\frac{P_{\text{i,loss}}}{V}\right)}^{\texttt{pden_ion_transport_loss_mw}} = \frac{3}{2}\frac{\frac{W_{\text{i}}}{V}}{\tau_{\text{E}}}
 $$
 
 $$
-\mathtt{pden\_electron\_transport\_loss\_mw} = \frac{3}{2}\frac{n_{\text{e}} \langle T_{\text{e}} \rangle_{\text{n}}}{\tau_{\text{E}}}
+\overbrace{\left(\frac{P_{\text{e,loss}}}{V}\right)}^{\texttt{pden_electron_transport_loss_mw}} = \frac{3}{2}\frac{\frac{W_{\text{e}}}{V}}{\tau_{\text{E}}}
 $$
 
-Here $\langle T_{\text{i}} \rangle$ and $\langle T_{\text{e}} \rangle$ are the ion and electron density weighted temperatures respectively.
+Here $W_{\text{e}}$ and $W_{\text{i}}$ are the electron and ion stored thermal energies respectively and $V$ is the plasma volume.
+
+
 
-Calculate the density and density weighted ratio:
+----------
+
+## Global confinement time consistency | `calculate_electron_species_consistent_energy_confinement_time()`
+
+It is possible to set a relative ratio between the electron and ion confinement times $\overbrace{R}^{\texttt{f_t_fuel_ion_electron_energy_confinement}} = \frac{\tau_\text{i}}{\tau_\text{e}}$ that still ensures the global energy confinment time is consistent.
+
+It is derived as:
+
+The global energy confinement time is defined as:
 
 $$
-\frac{n_{\text{i}}}{n_{\text{e}}}\frac{\langle T_{\text{i}} \rangle_{\text{n}}}{\langle T_{\text{e}} \rangle_{\text{n}}}
+\tau_{\text{E}} = \frac{W_{\text{e}} +W_{\text{i}}}{P_{\text{e}} +P_{\text{i}}}
 $$
 
-The density weighted global energy confinement time is then found in terms of this ratio:
+where $W_{\text{e}},W_{\text{i}}$ is the stored thermal energy of the electrons and ions respectively and $P_{\text{e}},P_{\text{i}}$ is the electron and ion energy transport powers.
+
+The electron and ion energy confinement times alone are:
 
 $$
-\tau_{\text{E}} = \frac{\frac{n_{\text{i}}}{n_{\text{e}}}\frac{\langle T_{\text{i}} \rangle_{\text{n}}}{\langle T_{\text{e}} \rangle_{\text{n}}} + 1}{\left(\frac{\frac{n_{\text{i}}}{n_{\text{e}}}\frac{\langle T_{\text{i}} \rangle_{\text{n}}}{\langle T_{\text{e}} \rangle_{\text{n}}}}{\tau_{\text{i}}}+\frac{1}{\tau_{\text{e}}}\right)}
+\tau_{\text{e}} = \frac{W_{\text{e}}}{P_{\text{e}}}
 $$
 
-----------
-
-## Key Constraints
+$$
+\tau_{\text{i}} = \frac{W_{\text{i}}}{P_{\text{i}}}
+$$
 
-### Global plasma power balance
+Substituting the individual ion and electron terms into the global term we get:
 
-This constraint can be activated by stating `icc = 2` in the input file.
+$$
+\tau_{\text{E}} = \frac{W_{\text{e}} +W_{\text{i}}}{\frac{W_{\text{e}}}{\tau_{\text{e}}} +\frac{W_{\text{i}}}{\tau_{\text{i}}}}
+$$
 
-This constraint ensures self consistency between the the transport loss power used for the confinement scalings and the calculated confinement time in relation to the plasmas total thermal energy:
+Substituting in our known ratio $\left(R = \frac{\tau_\text{i}}{\tau_\text{e}}\right)$:
 
 $$
-P_{\text{L}} = \frac{W}{\tau_{\text{E}}}
+\tau_{\text{E}} = \frac{W_{\text{e}} +W_{\text{i}}}{\frac{W_{\text{e}}}{\tau_{\text{e}}} +\frac{W_{\text{e}}}{R\tau_{\text{e}}}}
 $$
 
+Factor out the $\frac{1}{\tau_{\text{e}}}$ from the denominator:
+
 $$
-\underbrace{\frac{3}{2}\frac{n_{\text{i}} \langle T_{\text{i}} \rangle_{\text{n}}}{\tau_{\text{E}}} + \frac{3}{2}\frac{n_{\text{e}} \langle T_{\text{e}} \rangle_{\text{n}}}{\tau_{\text{E}}}}_{\frac{W}{\tau_{\text{E}}}}  = \underbrace{\frac{f_{\alpha}P_{\alpha} + P_{\text{c}} + P_{\text{OH}} + P_{\text{HCD}}}{V_{\text{P}}} - \frac{P_{\text{rad}}}{V_{\text{p}}}}_{P_{\text{L}}}
+\tau_{\text{E}} = \tau_{\text{e}}\left(\frac{W_{\text{e}} +W_{\text{i}}}{W_{\text{e}} +\frac{W_{\text{i}}}{R}}\right)
 $$
 
-The $\frac{3}{2}n_{\text{i}} \langle T_{\text{i}} \rangle_{\text{n}}$ value is simply the volume averaged ion thermal energy density where $\langle T_{\text{i}} \rangle_{\text{n}}$ is the density weighted temperature. The same goes for the $\frac{3}{2}n_{\text{e}} \langle T_{\text{e}} \rangle_{\text{e}}$ electron thermal energy density term. $\tau_{\text{E}}$ is the confinement time calculated from the chosen confinement scaling via `i_confinement_time`. 
+Finally re-arrange the equation to isolate $\tau_{\text{e}}$:
 
-The constraint uses the loss power and thermal densities hence the inclusion of the $V_{\text{p}}$ plasma volume term. The constraint is adapted depending on the condition of `i_rad_loss` which governs the radiation contribution to the loss power definition, see the [radiation and energy confinement section](#effect-of-radiation-on-energy-confinement) for more info. The injected heating and current drive contribution $P_{\text{HCD}}$ is also included or excluded depending if the plasma is deemed to be ignited with the `i_plasma_ignited` switch.
+$$
+\tau_{\text{e}} = \tau_{\text{E}}\left(\frac{W_{\text{e}} +\frac{W_{\text{i}}}{R}}{W_{\text{e}} +W_{\text{i}}}\right)
+$$
 
-**It is highly recommended to always have this constraint on as it is a global consistency checker**
+Using the above we can specify plasmas where $\tau_{\text{e}} \neq \tau_{\text{i}} \neq \tau_{\text{E}}$ and still ensure that our global plasma and species energy transport is self consistent.
 
 ----------
 
+## Key Constraints
+
 ### Lower limit on alpha particle confinement time ratio
 
 This constraint can be activated by stating `icc = 62` in the input file.
diff --git a/documentation/source/physics-models/plasma_power_balance.md b/documentation/source/physics-models/plasma_power_balance.md
index 0725625a3b..af83086e0e 100644
--- a/documentation/source/physics-models/plasma_power_balance.md
+++ b/documentation/source/physics-models/plasma_power_balance.md
@@ -3,10 +3,6 @@
 There are several equations that need to be satisfied in order to ensure that the plasma is in 
 power equilibrium. This includes the total energy and power leaving the plasma as a whole and then also the rate of energy transfer between the ion and electron species themselves.
 
-!!! warning "Treatment of particle energy confinement"
-
-    Currently the electrons and ions in `PROCESS` are assumed to have the same energy confinement time and which is equal to the calculated global energy confinement time, $\tau_{e} = \tau_{i} = \tau_{E}$.
-
 ## Global plasma power balance
 
 This constraint can be activated by stating `icc = 2` in the input file.
diff --git a/process/models/physics/confinement_time.py b/process/models/physics/confinement_time.py
index 6b031def36..25aa86d3e5 100644
--- a/process/models/physics/confinement_time.py
+++ b/process/models/physics/confinement_time.py
@@ -966,33 +966,17 @@ def calculate_confinement_time(
         elif self.data.physics.i_rad_loss == ConfinementRadiationLossModel.NO_RADIATION:
             hstar = hfact
 
-        # Calculation of the transport power loss terms
-        # Transport losses in Watts/m3 are 3/2 * n.e.T / tau , with T in eV
-        # (here, temp_plasma_ion_density_weighted_kev and
-        # temp_plasma_electron_density_weighted_kev are in keV, and
-        # pden_electron_transport_loss_mw and pden_ion_transport_loss_mw are in MW/m3)
-
         t_energy_confinement = t_electron_confinement * hfact
 
         # Apply H-factor correction to chosen scaling
         t_electron_energy_confinement = (
-            hfact
-            * t_energy_confinement
-            * (
-                (
-                    eden_plasma_electrons_thermal_vol_avg
-                    + (
-                        eden_plasma_ions_thermal_vol_avg
-                        / self.data.physics.f_t_fuel_ion_electron_energy_confinement
-                    )
-                )
-                / (
-                    eden_plasma_electrons_thermal_vol_avg
-                    + eden_plasma_ions_thermal_vol_avg
-                )
+            self.calculate_electron_species_consistent_energy_confinement_time(
+                self.data.physics.f_t_fuel_ion_electron_energy_confinement,
+                eden_plasma_ions_thermal_vol_avg,
+                eden_plasma_electrons_thermal_vol_avg,
+                t_energy_confinement,
             )
         )
-
         # Ion energy confinement time
         t_ion_energy_confinement = (
             t_electron_energy_confinement
@@ -1008,15 +992,6 @@ def calculate_confinement_time(
             eden_plasma_electrons_thermal_vol_avg / t_electron_energy_confinement
         ) / 1e6  # Convert from W/m³ to MW/m³
 
-        # ratio = (nd_plasma_ions_total_vol_avg / nd_plasma_electrons_vol_avg) * (
-        #     temp_plasma_ion_density_weighted_kev
-        #     / temp_plasma_electron_density_weighted_kev
-        # )
-
-        # Global energy confinement time
-
-        # t_energy_confinement = t_ion_energy_confinement
-
         # For comparison directly calculate the confinement time from the stored energy
         # calculated from the total plasma beta and the loss power used above.
         self.data.physics.t_energy_confinement_beta = (
@@ -1063,6 +1038,48 @@ def calculate_double_and_triple_product(
 
         return ntau, nTtau
 
+    @staticmethod
+    def calculate_electron_species_consistent_energy_confinement_time(
+        f_t_fuel_ion_electron_energy_confinement: float,
+        eden_plasma_ions_thermal_vol_avg: float,
+        eden_plasma_electrons_thermal_vol_avg: float,
+        t_plasma_global_energy_confinement: float,
+    ) -> float:
+        """Calculate the electron energy confinement time based on the global plasma
+        energy confinement time and the species ratio.
+
+        Parameters
+        ----------
+        f_t_fuel_ion_electron_energy_confinement : float
+            Ratio of fuel ion to electron energy confinement times.
+        eden_plasma_ions_thermal_vol_avg : float
+            Volume averaged thermal energy density of ions [J/m³].
+        eden_plasma_electrons_thermal_vol_avg : float
+            Volume averaged thermal energy density of electrons [J/m³].
+        t_plasma_global_energy_confinement : float
+            Global plasma energy confinement time [s].
+
+        Returns
+        -------
+        float
+            Electron energy confinement time [s].
+
+        Notes
+        -----
+        - The variable `f_t_fuel_ion_electron_energy_confinement` is used to relate
+        the energy confinement times of fuel ions and electrons. R = τ_i / τ_e.
+        """
+        return t_plasma_global_energy_confinement * (
+            (
+                eden_plasma_electrons_thermal_vol_avg
+                + (
+                    eden_plasma_ions_thermal_vol_avg
+                    / f_t_fuel_ion_electron_energy_confinement
+                )
+            )
+            / (eden_plasma_electrons_thermal_vol_avg + eden_plasma_ions_thermal_vol_avg)
+        )
+
     def find_other_h_factors(self, i_confinement_time: int) -> float:
         """Function to find H-factor for the equivalent confinement time in other
         scalings.

From e4e465e4d4da9dcb4057a8c76309f1562de29d6e Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Fri, 17 Jul 2026 16:19:59 +0100
Subject: [PATCH 31/37] Rename internal variable to be global instead of
 electron confinement

---
 process/models/physics/confinement_time.py | 136 +++++++++++----------
 1 file changed, 71 insertions(+), 65 deletions(-)

diff --git a/process/models/physics/confinement_time.py b/process/models/physics/confinement_time.py
index 25aa86d3e5..111e7e18e3 100644
--- a/process/models/physics/confinement_time.py
+++ b/process/models/physics/confinement_time.py
@@ -220,7 +220,7 @@ def calculate_confinement_time(
         if (
             model == ConfinementTimeModel.USER_INPUT
         ):  # t_electron_energy_confinement is an input
-            t_electron_confinement = self.data.physics.tauee_in
+            t_global_confinement_scaling = self.data.physics.tauee_in
 
         # ========================================================================
 
@@ -228,7 +228,7 @@ def calculate_confinement_time(
         if (
             model == ConfinementTimeModel.NEO_ALCATOR
         ):  # t_electron_energy_confinement is an input
-            t_electron_confinement = self.neo_alcator_confinement_time(
+            t_global_confinement_scaling = self.neo_alcator_confinement_time(
                 n20, rminor, rmajor, qstar
             )
 
@@ -236,13 +236,15 @@ def calculate_confinement_time(
 
         # "Mirnov"-like scaling (H-mode)
         elif model == ConfinementTimeModel.MIRNOV:  # Mirnov scaling (H-mode)
-            t_electron_confinement = self.mirnov_confinement_time(rminor, kappa95, pcur)
+            t_global_confinement_scaling = self.mirnov_confinement_time(
+                rminor, kappa95, pcur
+            )
 
         # ========================================================================
 
         # Merezhkin-Mukhovatov (MM) OH/L-mode scaling
         elif model == ConfinementTimeModel.MEREZHKIN_MUHKOVATOV:
-            t_electron_confinement = self.merezhkin_muhkovatov_confinement_time(
+            t_global_confinement_scaling = self.merezhkin_muhkovatov_confinement_time(
                 rmajor,
                 rminor,
                 kappa95,
@@ -256,7 +258,7 @@ def calculate_confinement_time(
 
         # Shimomura (S) optimized H-mode scaling
         elif model == ConfinementTimeModel.SHIMOMURA:
-            t_electron_confinement = self.shimomura_confinement_time(
+            t_global_confinement_scaling = self.shimomura_confinement_time(
                 rmajor, rminor, b_plasma_toroidal_on_axis, kappa95, m_fuel_amu
             )
 
@@ -264,7 +266,7 @@ def calculate_confinement_time(
 
         # Kaye-Goldston scaling (L-mode)
         elif model == ConfinementTimeModel.KAYE_GOLDSTON:
-            t_electron_confinement = self.kaye_goldston_confinement_time(
+            t_global_confinement_scaling = self.kaye_goldston_confinement_time(
                 pcur,
                 rmajor,
                 rminor,
@@ -279,7 +281,7 @@ def calculate_confinement_time(
 
         # ITER Power scaling - ITER 89-P (L-mode)
         elif model == ConfinementTimeModel.ITER_89P:
-            t_electron_confinement = self.iter_89p_confinement_time(
+            t_global_confinement_scaling = self.iter_89p_confinement_time(
                 pcur,
                 rmajor,
                 rminor,
@@ -294,7 +296,7 @@ def calculate_confinement_time(
 
         # ITER Offset linear scaling - ITER 89-O (L-mode)
         elif model == ConfinementTimeModel.ITER_89_0:
-            t_electron_confinement = self.iter_89_0_confinement_time(
+            t_global_confinement_scaling = self.iter_89_0_confinement_time(
                 pcur,
                 rmajor,
                 rminor,
@@ -308,7 +310,7 @@ def calculate_confinement_time(
 
         # Rebut-Lallia offset linear scaling (L-mode)
         elif model == ConfinementTimeModel.REBUT_LALLIA:
-            t_electron_confinement = self.rebut_lallia_confinement_time(
+            t_global_confinement_scaling = self.rebut_lallia_confinement_time(
                 rminor,
                 rmajor,
                 kappa,
@@ -324,7 +326,7 @@ def calculate_confinement_time(
 
         # Goldston scaling (L-mode)
         elif model == ConfinementTimeModel.GOLDSTON:  # Goldston scaling (L-mode)
-            t_electron_confinement = self.goldston_confinement_time(
+            t_global_confinement_scaling = self.goldston_confinement_time(
                 pcur, rmajor, rminor, kappa95, m_fuel_amu, p_plasma_loss_mw
             )
 
@@ -332,7 +334,7 @@ def calculate_confinement_time(
 
         # T-10 scaling (L-mode)
         elif model == ConfinementTimeModel.T_10:
-            t_electron_confinement = self.t10_confinement_time(
+            t_global_confinement_scaling = self.t10_confinement_time(
                 dnla20,
                 rmajor,
                 qstar,
@@ -348,7 +350,7 @@ def calculate_confinement_time(
 
         # JAERI / Odajima-Shimomura L-mode scaling
         elif model == ConfinementTimeModel.JAERI:  # JAERI scaling
-            t_electron_confinement = self.jaeri_confinement_time(
+            t_global_confinement_scaling = self.jaeri_confinement_time(
                 kappa95,
                 rminor,
                 m_fuel_amu,
@@ -365,7 +367,7 @@ def calculate_confinement_time(
 
         # Kaye "big"  L-mode scaling (based only on big tokamak data)
         elif model == ConfinementTimeModel.KAYE_BIG:
-            t_electron_confinement = self.kaye_big_confinement_time(
+            t_global_confinement_scaling = self.kaye_big_confinement_time(
                 rmajor,
                 rminor,
                 b_plasma_toroidal_on_axis,
@@ -380,7 +382,7 @@ def calculate_confinement_time(
 
         # ITER H90-P H-mode scaling
         elif model == ConfinementTimeModel.ITER_H90_P:
-            t_electron_confinement = self.iter_h90_p_confinement_time(
+            t_global_confinement_scaling = self.iter_h90_p_confinement_time(
                 pcur,
                 rmajor,
                 rminor,
@@ -395,7 +397,7 @@ def calculate_confinement_time(
 
         # Minimum of ITER 89-P and ITER 89-O
         elif model == ConfinementTimeModel.MINIMUM_OF_ITER_89P_AND_ITER_89_0:
-            t_electron_confinement = min(
+            t_global_confinement_scaling = min(
                 self.iter_89p_confinement_time(
                     pcur,
                     rmajor,
@@ -422,7 +424,7 @@ def calculate_confinement_time(
 
         # Riedel scaling (L-mode)
         elif model == ConfinementTimeModel.RIEDEL_L:
-            t_electron_confinement = self.riedel_l_confinement_time(
+            t_global_confinement_scaling = self.riedel_l_confinement_time(
                 pcur,
                 rmajor,
                 rminor,
@@ -436,7 +438,7 @@ def calculate_confinement_time(
 
         # Christiansen et al scaling (L-mode)
         elif model == ConfinementTimeModel.CHRISTIANSEN:
-            t_electron_confinement = self.christiansen_confinement_time(
+            t_global_confinement_scaling = self.christiansen_confinement_time(
                 pcur,
                 rmajor,
                 rminor,
@@ -451,7 +453,7 @@ def calculate_confinement_time(
 
         # Lackner-Gottardi scaling (L-mode)
         elif model == ConfinementTimeModel.LACKNER_GOTTARDI:
-            t_electron_confinement = self.lackner_gottardi_confinement_time(
+            t_global_confinement_scaling = self.lackner_gottardi_confinement_time(
                 pcur,
                 rmajor,
                 rminor,
@@ -465,7 +467,7 @@ def calculate_confinement_time(
 
         # Neo-Kaye scaling (L-mode)
         elif model == ConfinementTimeModel.NEO_KAYE:
-            t_electron_confinement = self.neo_kaye_confinement_time(
+            t_global_confinement_scaling = self.neo_kaye_confinement_time(
                 pcur,
                 rmajor,
                 rminor,
@@ -479,7 +481,7 @@ def calculate_confinement_time(
 
         # Riedel scaling (H-mode)
         elif model == ConfinementTimeModel.RIEDEL_H:
-            t_electron_confinement = self.riedel_h_confinement_time(
+            t_global_confinement_scaling = self.riedel_h_confinement_time(
                 pcur,
                 rmajor,
                 rminor,
@@ -494,7 +496,7 @@ def calculate_confinement_time(
 
         # Amended version of ITER H90-P law
         elif model == ConfinementTimeModel.ITER_H90_P_AMENDED:
-            t_electron_confinement = self.iter_h90_p_amended_confinement_time(
+            t_global_confinement_scaling = self.iter_h90_p_amended_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 m_fuel_amu,
@@ -507,7 +509,7 @@ def calculate_confinement_time(
 
         # Sudo et al. scaling (stellarators/heliotron)
         elif model == ConfinementTimeModel.SUDO_ET_AL:
-            t_electron_confinement = self.sudo_et_al_confinement_time(
+            t_global_confinement_scaling = self.sudo_et_al_confinement_time(
                 rmajor,
                 rminor,
                 dnla20,
@@ -519,7 +521,7 @@ def calculate_confinement_time(
 
         # Gyro-reduced Bohm scaling
         elif model == ConfinementTimeModel.GYRO_REDUCED_BOHM:
-            t_electron_confinement = self.gyro_reduced_bohm_confinement_time(
+            t_global_confinement_scaling = self.gyro_reduced_bohm_confinement_time(
                 b_plasma_toroidal_on_axis,
                 dnla20,
                 p_plasma_loss_mw,
@@ -531,20 +533,22 @@ def calculate_confinement_time(
 
         # Lackner-Gottardi stellarator scaling
         elif model == ConfinementTimeModel.LACKNER_GOTTARDI_STELLARATOR:
-            t_electron_confinement = self.lackner_gottardi_stellarator_confinement_time(
-                rmajor,
-                rminor,
-                dnla20,
-                b_plasma_toroidal_on_axis,
-                p_plasma_loss_mw,
-                q95,
+            t_global_confinement_scaling = (
+                self.lackner_gottardi_stellarator_confinement_time(
+                    rmajor,
+                    rminor,
+                    dnla20,
+                    b_plasma_toroidal_on_axis,
+                    p_plasma_loss_mw,
+                    q95,
+                )
             )
 
         # ==========================================================================
 
         # ITER_93 ELM-free H-mode scaling
         elif model == ConfinementTimeModel.ITER_93H:
-            t_electron_confinement = self.iter_93h_confinement_time(
+            t_global_confinement_scaling = self.iter_93h_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 p_plasma_loss_mw,
@@ -563,7 +567,7 @@ def calculate_confinement_time(
 
         # ELM-free: ITERH-97P
         elif model == ConfinementTimeModel.ITER_H97P:
-            t_electron_confinement = self.iter_h97p_confinement_time(
+            t_global_confinement_scaling = self.iter_h97p_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 p_plasma_loss_mw,
@@ -578,7 +582,7 @@ def calculate_confinement_time(
 
         # ELMy: ITERH-97P(y)
         elif model == ConfinementTimeModel.ITER_H97P_ELMY:
-            t_electron_confinement = self.iter_h97p_elmy_confinement_time(
+            t_global_confinement_scaling = self.iter_h97p_elmy_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 p_plasma_loss_mw,
@@ -593,7 +597,7 @@ def calculate_confinement_time(
 
         # ITER-96P (= ITER-97L) L-mode scaling
         elif model == ConfinementTimeModel.ITER_96P:
-            t_electron_confinement = self.iter_96p_confinement_time(
+            t_global_confinement_scaling = self.iter_96p_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 kappa95,
@@ -609,7 +613,7 @@ def calculate_confinement_time(
         # Valovic modified ELMy-H mode scaling
         # WARNING: No reference found for this scaling. This may not be its real name
         elif model == ConfinementTimeModel.VALOVIC_ELMY:
-            t_electron_confinement = self.valovic_elmy_confinement_time(
+            t_global_confinement_scaling = self.valovic_elmy_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 dnla19,
@@ -625,7 +629,7 @@ def calculate_confinement_time(
         # Kaye PPPL Workshop April 1998 L-mode scaling
         # WARNING: No reference found for this scaling. This may not be its real name
         elif model == ConfinementTimeModel.KAYE:
-            t_electron_confinement = self.kaye_confinement_time(
+            t_global_confinement_scaling = self.kaye_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 kappa,
@@ -641,7 +645,7 @@ def calculate_confinement_time(
         # ITERH-PB98P(y), ELMy H-mode scaling
         # WARNING: No reference found for this scaling. This may not be its real name
         elif model == ConfinementTimeModel.ITER_PB98P_Y:
-            t_electron_confinement = self.iter_pb98py_confinement_time(
+            t_global_confinement_scaling = self.iter_pb98py_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 dnla19,
@@ -656,7 +660,7 @@ def calculate_confinement_time(
 
         # IPB98(y), ELMy H-mode scaling
         elif model == ConfinementTimeModel.IPB98_Y:
-            t_electron_confinement = self.iter_ipb98y_confinement_time(
+            t_global_confinement_scaling = self.iter_ipb98y_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 dnla19,
@@ -671,7 +675,7 @@ def calculate_confinement_time(
 
         # IPB98(y,1), ELMy H-mode scaling
         elif model == ConfinementTimeModel.ITER_IPB98Y1:
-            t_electron_confinement = self.iter_ipb98y1_confinement_time(
+            t_global_confinement_scaling = self.iter_ipb98y1_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 dnla19,
@@ -686,7 +690,7 @@ def calculate_confinement_time(
 
         # IPB98(y,2), ELMy H-mode scaling
         elif model == ConfinementTimeModel.ITER_IPB98Y2:
-            t_electron_confinement = self.iter_ipb98y2_confinement_time(
+            t_global_confinement_scaling = self.iter_ipb98y2_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 dnla19,
@@ -701,7 +705,7 @@ def calculate_confinement_time(
 
         # IPB98(y,3), ELMy H-mode scaling
         elif model == ConfinementTimeModel.ITER_IPB98Y3:
-            t_electron_confinement = self.iter_ipb98y3_confinement_time(
+            t_global_confinement_scaling = self.iter_ipb98y3_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 dnla19,
@@ -716,7 +720,7 @@ def calculate_confinement_time(
 
         # IPB98(y,4), ELMy H-mode scaling
         elif model == ConfinementTimeModel.ITER_IPB98Y4:
-            t_electron_confinement = self.iter_ipb98y4_confinement_time(
+            t_global_confinement_scaling = self.iter_ipb98y4_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 dnla19,
@@ -733,7 +737,7 @@ def calculate_confinement_time(
         elif model == ConfinementTimeModel.ISS95_STELLARATOR:
             # dummy argument q95 is actual argument iotabar for stellarators
             iotabar = q95
-            t_electron_confinement = self.iss95_stellarator_confinement_time(
+            t_global_confinement_scaling = self.iss95_stellarator_confinement_time(
                 rminor,
                 rmajor,
                 dnla19,
@@ -748,7 +752,7 @@ def calculate_confinement_time(
         elif model == ConfinementTimeModel.ISS04_STELLARATOR:
             # dummy argument q95 is actual argument iotabar for stellarators
             iotabar = q95
-            t_electron_confinement = self.iss04_stellarator_confinement_time(
+            t_global_confinement_scaling = self.iss04_stellarator_confinement_time(
                 rminor,
                 rmajor,
                 dnla19,
@@ -761,7 +765,7 @@ def calculate_confinement_time(
 
         # DS03 beta-independent H-mode scaling
         elif model == ConfinementTimeModel.DS03:
-            t_electron_confinement = self.ds03_confinement_time(
+            t_global_confinement_scaling = self.ds03_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 dnla19,
@@ -776,7 +780,7 @@ def calculate_confinement_time(
 
         #  Murari "Non-power law" scaling
         elif model == ConfinementTimeModel.MURARI:
-            t_electron_confinement = self.murari_confinement_time(
+            t_global_confinement_scaling = self.murari_confinement_time(
                 pcur,
                 rmajor,
                 self.data.physics.kappa_ipb,
@@ -789,7 +793,7 @@ def calculate_confinement_time(
 
         # Petty08, beta independent dimensionless scaling
         elif model == ConfinementTimeModel.PETTY08:
-            t_electron_confinement = self.petty08_confinement_time(
+            t_global_confinement_scaling = self.petty08_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 dnla19,
@@ -803,7 +807,7 @@ def calculate_confinement_time(
 
         # Lang high density relevant confinement scaling
         elif model == ConfinementTimeModel.LANG_HIGH_DENSITY:
-            t_electron_confinement = self.lang_high_density_confinement_time(
+            t_global_confinement_scaling = self.lang_high_density_confinement_time(
                 plasma_current,
                 b_plasma_toroidal_on_axis,
                 nd_plasma_electron_line,
@@ -821,7 +825,7 @@ def calculate_confinement_time(
 
         # Hubbard 2017 I-mode confinement time scaling - nominal
         elif model == ConfinementTimeModel.HUBBARD_NOMINAL:
-            t_electron_confinement = self.hubbard_nominal_confinement_time(
+            t_global_confinement_scaling = self.hubbard_nominal_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 dnla20,
@@ -832,7 +836,7 @@ def calculate_confinement_time(
 
         # Hubbard 2017 I-mode confinement time scaling - lower
         elif model == ConfinementTimeModel.HUBBARD_LOWER:
-            t_electron_confinement = self.hubbard_lower_confinement_time(
+            t_global_confinement_scaling = self.hubbard_lower_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 dnla20,
@@ -843,7 +847,7 @@ def calculate_confinement_time(
 
         # Hubbard 2017 I-mode confinement time scaling - upper
         elif model == ConfinementTimeModel.HUBBARD_UPPER:
-            t_electron_confinement = self.hubbard_upper_confinement_time(
+            t_global_confinement_scaling = self.hubbard_upper_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 dnla20,
@@ -854,7 +858,7 @@ def calculate_confinement_time(
 
         # Menard NSTX, ELMy H-mode scaling
         elif model == ConfinementTimeModel.MENARD_NSTX:
-            t_electron_confinement = self.menard_nstx_confinement_time(
+            t_global_confinement_scaling = self.menard_nstx_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 dnla19,
@@ -869,22 +873,24 @@ def calculate_confinement_time(
 
         # Menard NSTX-Petty08 Hybrid
         elif model == ConfinementTimeModel.MENARD_NSTX_PETTY08_HYBRID:
-            t_electron_confinement = self.menard_nstx_petty08_hybrid_confinement_time(
-                pcur,
-                b_plasma_toroidal_on_axis,
-                dnla19,
-                p_plasma_loss_mw,
-                rmajor,
-                self.data.physics.kappa_ipb,
-                aspect,
-                m_fuel_amu,
+            t_global_confinement_scaling = (
+                self.menard_nstx_petty08_hybrid_confinement_time(
+                    pcur,
+                    b_plasma_toroidal_on_axis,
+                    dnla19,
+                    p_plasma_loss_mw,
+                    rmajor,
+                    self.data.physics.kappa_ipb,
+                    aspect,
+                    m_fuel_amu,
+                )
             )
 
         # ==========================================================================
 
         # NSTX gyro-Bohm (Buxton)
         elif model == ConfinementTimeModel.NSTX_GYRO_BOHM:
-            t_electron_confinement = self.nstx_gyro_bohm_confinement_time(
+            t_global_confinement_scaling = self.nstx_gyro_bohm_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 p_plasma_loss_mw,
@@ -896,7 +902,7 @@ def calculate_confinement_time(
 
         # ITPA20 H-mode scaling
         elif model == ConfinementTimeModel.ITPA20:
-            t_electron_confinement = self.itpa20_confinement_time(
+            t_global_confinement_scaling = self.itpa20_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 dnla19,
@@ -912,7 +918,7 @@ def calculate_confinement_time(
 
         # ITPA20-IL confinement time scaling
         elif model == ConfinementTimeModel.ITPA20_IL:
-            t_electron_confinement = self.itpa20_il_confinement_time(
+            t_global_confinement_scaling = self.itpa20_il_confinement_time(
                 pcur,
                 b_plasma_toroidal_on_axis,
                 p_plasma_loss_mw,
@@ -966,7 +972,7 @@ def calculate_confinement_time(
         elif self.data.physics.i_rad_loss == ConfinementRadiationLossModel.NO_RADIATION:
             hstar = hfact
 
-        t_energy_confinement = t_electron_confinement * hfact
+        t_energy_confinement = t_global_confinement_scaling * hfact
 
         # Apply H-factor correction to chosen scaling
         t_electron_energy_confinement = (

From bb1097adc6b35e6529f2495cb201c5ade698854b Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 20 Jul 2026 10:32:35 +0100
Subject: [PATCH 32/37] Update confinement summary box to display new variables

---
 process/core/io/plot/summary.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index de399365a6..e49fd5bf97 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -2818,17 +2818,17 @@ def plot_main_plasma_information(
         f"$\\mathbf{{Confinement:}}$\n \n"
         f"Confinement scaling law: {mfile.get('tauelaw', scan=scan)}\n"
         f"Confinement $H$ factor: {mfile.get('hfact', scan=scan):.4f}\n"
-        f"Energy confinement time from scaling: {mfile.get('t_energy_confinement', scan=scan):.4f} s\n"
+        f"$\\tau_{{\\text{{E}}}}$: {mfile.get('t_energy_confinement', scan=scan):.4f} s | $\\tau_{{\\text{{e}}}}$: {mfile.get('t_electron_energy_confinement', scan=scan):.4f} s | $\\tau_{{\\text{{i}}}}$: {mfile.get('t_ion_energy_confinement', scan=scan):.4f} s | $\\tau_{{\\text{{i}}}}/\\tau_{{\\text{{e}}}}$: {mfile.get('f_t_fuel_ion_electron_energy_confinement', scan=scan):.4f} \n"
         f"Fusion double product: {mfile.get('ntau', scan=scan):.4e} s/m³\n"
         f"Lawson Triple product: {mfile.get('nttau', scan=scan):.4e} keV·s/m³\n"
         f"Transport loss power assumed in scaling law: {mfile.get('p_plasma_loss_mw', scan=scan):.4f} MW\n"
-        f"Plasma thermal energy (inc. $\\alpha$), $W$: {mfile.get('e_plasma_beta', scan=scan) / 1e9:.4f} GJ\n"
-        f"Alpha particle confinement time: {mfile.get('t_alpha_confinement', scan=scan):.4f} s | $\\tau_{{\\alpha}}/\\tau_{{e}}$: {mfile.get('f_t_alpha_energy_confinement', scan=scan):.4f}"
+        f"$W_{{\\beta}}$: {mfile.get('e_plasma_beta', scan=scan) / 1e9:.4f} GJ | $W_{{th}}$: {mfile.get('e_plasma_thermal_total', scan=scan) / 1e9:.4f} GJ | $W_{{e}}$: {mfile.get('e_plasma_electrons_thermal', scan=scan) / 1e9:.4f} GJ | $W_{{i}}$: {mfile.get('e_plasma_ions_thermal', scan=scan) / 1e9:.4f} GJ\n"
+        f"Particle confinement: $\\tau_{{\\alpha}}$: {mfile.get('t_alpha_confinement', scan=scan):.4f} s | $\\tau_{{\\alpha}}/\\tau_{{e}}$: {mfile.get('f_t_alpha_energy_confinement', scan=scan):.4f}"
     )
 
     axis.text(
         0.025,
-        0.57,
+        0.575,
         textstr_confinement,
         fontsize=9,
         verticalalignment="top",

From 4edf53069ac9a48562e7fddbb9c45de4de22cc67 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 20 Jul 2026 11:46:59 +0100
Subject: [PATCH 33/37] Refactor constraint_equation_2 to streamline plasma
 loss calculations and improve clarity

---
 process/core/solver/constraints.py | 33 ++++++++++--------------------
 1 file changed, 11 insertions(+), 22 deletions(-)

diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py
index 7ba7a9b310..a200fb8bdc 100644
--- a/process/core/solver/constraints.py
+++ b/process/core/solver/constraints.py
@@ -279,36 +279,25 @@ def constraint_equation_2(constraint_registration, data):
     """
     # pscaling: total transport power per volume (MW/m3)
 
-    e_ions_stored_thermal = (
-        (3 / 2)
-        * (constants.ELECTRON_CHARGE / 1e3)
-        * data.physics.nd_plasma_ions_total_vol_avg
-        * data.physics.temp_plasma_ion_density_weighted_kev
-    )
-
-    e_electrons_stored_thermal = (
-        (3 / 2)
-        * (constants.ELECTRON_CHARGE / 1e3)
-        * data.physics.nd_plasma_electrons_vol_avg
-        * data.physics.temp_plasma_electron_density_weighted_kev
-    )
-    pscaling = (
-        e_ions_stored_thermal + e_electrons_stored_thermal
+    p_plasma_loss = (
+        (
+            data.physics.eden_plasma_electrons_thermal_vol_avg
+            + data.physics.eden_plasma_ions_thermal_vol_avg
+        )
+        / 1e6
     ) / data.physics.t_energy_confinement
     # Total power lost is scaling power plus radiation:
     if data.physics.i_rad_loss == 0:
-        pnumerator = pscaling + data.physics.pden_plasma_rad_mw
+        p_plasma_loss += data.physics.pden_plasma_rad_mw
     elif data.physics.i_rad_loss == 1:
-        pnumerator = pscaling + data.physics.pden_plasma_core_rad_mw
-    else:
-        pnumerator = pscaling
+        p_plasma_loss += data.physics.pden_plasma_core_rad_mw
 
     # if plasma not ignited include injected power
     if (
         PlasmaIgnitionModel(data.physics.i_plasma_ignited)
         == PlasmaIgnitionModel.NON_IGNITED
     ):
-        pdenom = (
+        p_plasma_heating = (
             data.physics.f_p_alpha_plasma_deposited * data.physics.pden_alpha_total_mw
             + data.physics.pden_non_alpha_charged_mw
             + data.physics.pden_plasma_ohmic_mw
@@ -316,13 +305,13 @@ def constraint_equation_2(constraint_registration, data):
         )
     else:
         # if plasma ignited
-        pdenom = (
+        p_plasma_heating = (
             data.physics.f_p_alpha_plasma_deposited * data.physics.pden_alpha_total_mw
             + data.physics.pden_non_alpha_charged_mw
             + data.physics.pden_plasma_ohmic_mw
         )
 
-    return eq(pnumerator, pdenom, constraint_registration)
+    return eq(p_plasma_loss, p_plasma_heating, constraint_registration)
 
 
 @ConstraintManager.register_constraint(3, "MW/m³", "=")

From 59580d3d4ae9d68b76e6dfe89568be860da8de07 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 20 Jul 2026 13:53:06 +0100
Subject: [PATCH 34/37] Add warnings for global and species power balance
 constraints in documentation and validation checks

---
 .../source/physics-models/plasma_power_balance.md     |  6 ++++++
 process/core/init.py                                  | 11 +++++++++++
 2 files changed, 17 insertions(+)

diff --git a/documentation/source/physics-models/plasma_power_balance.md b/documentation/source/physics-models/plasma_power_balance.md
index af83086e0e..6cd5e9f588 100644
--- a/documentation/source/physics-models/plasma_power_balance.md
+++ b/documentation/source/physics-models/plasma_power_balance.md
@@ -3,6 +3,12 @@
 There are several equations that need to be satisfied in order to ensure that the plasma is in 
 power equilibrium. This includes the total energy and power leaving the plasma as a whole and then also the rate of energy transfer between the ion and electron species themselves.
 
+!!! warning "Global and species power balance"
+
+    It is not possible to have `icc = 2` for the global plasma power balance on at the same time as the ion `icc = 3` and electron `icc = 4` power balance. This is due to a solver error where the close to convergence the equality Jacobian loses rank. This is because the solver sees only two independant equality conditions but is asked to enforce three.
+
+    Using the ion and electron power balance constraints is preffered and still enforces the global power balance consistency. 
+
 ## Global plasma power balance
 
 This constraint can be activated by stating `icc = 2` in the input file.
diff --git a/process/core/init.py b/process/core/init.py
index 483719dd9e..4a5d108d6c 100644
--- a/process/core/init.py
+++ b/process/core/init.py
@@ -1160,6 +1160,17 @@ def check_process(inputs, data):  # noqa: ARG001
             "turn off CS temperature margin constraint icc = 60 when using REBCO"
         )
 
+    # Cannot use temperature margin constraint with REBCO CS coils
+    if (
+        (data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns] == 2).any()
+        and (data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns] == 3).any()
+        and (data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns] == 4).any()
+    ):
+        raise ProcessValidationError(
+            "Cannot have global power balance on as well as individual particle power "
+            "balance, this will cause a solver error. Either use icc = 2 alone or icc = 3 and 4 alone"
+        )
+
     # Cold end of the cryocooler should be colder than the TF
     if data.tfcoil.temp_tf_cryo > data.tfcoil.tftmp:
         raise ProcessValidationError("temp_tf_cryo should be lower than tftmp")

From 2cec7ac31e2d618c7f29758996ecd917010c6fb0 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 20 Jul 2026 14:34:54 +0100
Subject: [PATCH 35/37] Apply main structure to the plasma power balance output

---
 process/models/physics/confinement_time.py |   7 ++
 process/models/physics/physics.py          | 106 ++++++++++++++++++---
 2 files changed, 102 insertions(+), 11 deletions(-)

diff --git a/process/models/physics/confinement_time.py b/process/models/physics/confinement_time.py
index 111e7e18e3..d7c8da094d 100644
--- a/process/models/physics/confinement_time.py
+++ b/process/models/physics/confinement_time.py
@@ -1239,6 +1239,13 @@ def output_confinement_time_info(self):
         po.ocmmnt(self.outfile, "----------------------------")
         po.oblnkl(self.outfile)
 
+        po.ovarre(
+            self.outfile,
+            "Ratio of fuel ion to electron energy confinement times",
+            "(f_t_fuel_ion_electron_energy_confinement)",
+            self.data.physics.f_t_fuel_ion_electron_energy_confinement,
+            "OP ",
+        )
         po.ovarre(
             self.outfile,
             "Ion energy confinement time, from scaling (s)",
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 59580f05f4..3bf3e74c7c 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -2379,8 +2379,12 @@ def output_plasma_power_balance(self) -> None:
         """Output information about plasma power balance."""
         po.oheadr(self.outfile, "Plasma Power Balance")
 
-        po.osubhd(self.outfile, "Global power balance :")
+        po.ocmmnt(self.outfile, "Global power balance")
+        po.ocmmnt(self.outfile, "----------------------------")
+        po.oblnkl(self.outfile)
 
+        po.ocmmnt(self.outfile, "Heating powers:")
+        po.oblnkl(self.outfile)
         po.ovarre(
             self.outfile,
             "Ohmic heating power (MW)",
@@ -2425,17 +2429,44 @@ def output_plasma_power_balance(self) -> None:
             "OP ",
         )
         po.oblnkl(self.outfile)
+        po.ocmmnt(self.outfile, "Loss powers:")
+        po.oblnkl(self.outfile)
+
+        po.ovarre(
+            self.outfile,
+            "Total transport loss power (Pₗ) [MW]",
+            "(p_plasma_loss_mw)",
+            self.data.physics.p_plasma_loss_mw,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "Total radiated power from the plasma (Pᵧ) [MW]",
+            "(p_plasma_rad_mw)",
+            self.data.physics.p_plasma_rad_mw,
+            "OP ",
+        )
 
+        po.oblnkl(self.outfile)
+        po.ocmmnt(self.outfile, "----------------------------")
         po.osubhd(self.outfile, "Electron power balance :")
+        po.ocmmnt(self.outfile, "Heating powers:")
+        po.oblnkl(self.outfile)
 
         po.ovarre(
             self.outfile,
-            "Electron transport (MW)",
-            "(p_electron_transport_loss_mw)",
-            self.data.physics.p_electron_transport_loss_mw,
+            "Ohmic heating power to electrons (MW)",
+            "(p_plasma_ohmic_mw)",
+            self.data.physics.p_plasma_ohmic_mw,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "Alpha heating power to electrons (MW)",
+            "(pden_alpha_heating_electrons_mw)",
+            self.data.physics.pden_alpha_heating_electrons_mw,
             "OP ",
         )
-
         po.ovarre(
             self.outfile,
             "Injection power to electrons (MW)",
@@ -2443,6 +2474,18 @@ def output_plasma_power_balance(self) -> None:
             self.data.current_drive.p_hcd_injected_electrons_mw,
             "OP ",
         )
+
+        po.oblnkl(self.outfile)
+        po.ocmmnt(self.outfile, "Loss powers:")
+        po.oblnkl(self.outfile)
+
+        po.ovarre(
+            self.outfile,
+            "Electron transport (MW)",
+            "(p_electron_transport_loss_mw)",
+            self.data.physics.p_electron_transport_loss_mw,
+            "OP ",
+        )
         po.ovarre(
             self.outfile,
             "Ion-electron equilibration power (MW)",
@@ -2450,15 +2493,30 @@ def output_plasma_power_balance(self) -> None:
             self.data.physics.p_ion_electron_equilibration_vol_avg_mw,
             "OP ",
         )
+        po.oblnkl(self.outfile)
+        p_electron_balance = (
+            self.data.physics.p_plasma_ohmic_mw
+            + self.data.physics.pden_alpha_heating_electrons_mw
+            + self.data.current_drive.p_hcd_injected_electrons_mw
+            - self.data.physics.p_electron_transport_loss_mw
+            - self.data.physics.p_ion_electron_equilibration_vol_avg_mw
+        )
+        po.oblnkl(self.outfile)
+        po.ocmmnt(self.outfile, "----------------------------")
+        po.oblnkl(self.outfile)
         po.ovarre(
             self.outfile,
-            "Ratio of fuel ion to electron energy confinement times",
-            "(f_t_fuel_ion_electron_energy_confinement)",
-            self.data.physics.f_t_fuel_ion_electron_energy_confinement,
+            "Electron power balance (MW)",
+            "(p_electron_balance_mw)",
+            p_electron_balance,
             "OP ",
         )
+        po.oblnkl(self.outfile)
+        po.ocmmnt(self.outfile, "----------------------------")
         po.osubhd(self.outfile, "Ions power balance :")
 
+        po.ocmmnt(self.outfile, "Heating powers:")
+        po.oblnkl(self.outfile)
         po.ovarrf(
             self.outfile,
             "Fraction of alpha power to ions",
@@ -2468,11 +2526,12 @@ def output_plasma_power_balance(self) -> None:
         )
         po.ovarre(
             self.outfile,
-            "Ion transport (MW)",
-            "(p_ion_transport_loss_mw)",
-            self.data.physics.p_ion_transport_loss_mw,
+            "Alpha heating power to ions (MW)",
+            "(pden_alpha_heating_ions_mw)",
+            self.data.physics.pden_alpha_heating_ions_mw,
             "OP ",
         )
+
         po.ovarre(
             self.outfile,
             "Injection power to ions (MW)",
@@ -2480,6 +2539,31 @@ def output_plasma_power_balance(self) -> None:
             self.data.current_drive.p_hcd_injected_ions_mw,
             "OP ",
         )
+        po.oblnkl(self.outfile)
+        po.ocmmnt(self.outfile, "Loss powers:")
+        po.oblnkl(self.outfile)
+        po.ovarre(
+            self.outfile,
+            "Ion transport (MW)",
+            "(p_ion_transport_loss_mw)",
+            self.data.physics.p_ion_transport_loss_mw,
+            "OP ",
+        )
+        p_ion_balance = (
+            self.data.physics.pden_alpha_heating_ions_mw
+            + self.data.current_drive.p_hcd_injected_ions_mw
+            - self.data.physics.p_ion_transport_loss_mw
+            - self.data.physics.p_ion_electron_equilibration_vol_avg_mw
+        )
+        po.oblnkl(self.outfile)
+        po.ocmmnt(self.outfile, "----------------------------")
+        po.ovarre(
+            self.outfile,
+            "Ion power balance (MW)",
+            "(p_ion_balance_mw)",
+            p_ion_balance,
+            "OP ",
+        )
 
     def output_temperature_density_profile_info(self) -> None:
         """Output information about plasma temperature and density profiles."""

From e1654d42c878f3247fade382d788672475114d5f Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 20 Jul 2026 15:26:20 +0100
Subject: [PATCH 36/37] Refactor alpha heating power calculations and update
 output variables for electrons and ions

---
 process/core/solver/constraints.py          |  1 -
 process/data_structure/physics_variables.py |  6 +++
 process/models/physics/physics.py           | 44 +++++++++++++++++----
 3 files changed, 43 insertions(+), 8 deletions(-)

diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py
index a200fb8bdc..b40dc8c8b5 100644
--- a/process/core/solver/constraints.py
+++ b/process/core/solver/constraints.py
@@ -427,7 +427,6 @@ def constraint_equation_4(constraint_registration, data):
             data.physics.pden_alpha_heating_electrons_mw
             + data.physics.pden_plasma_ohmic_mw
         )
-
     return eq(p_electron_loss, p_electron_heating, constraint_registration)
 
 
diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 1ddd777de9..35724cca4a 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1028,12 +1028,18 @@ class PhysicsData:
     pden_alpha_heating_electrons_mw: float = 0.0
     """Alpha power heating per unit volume to electrons [MW/m³]"""
 
+    p_alpha_heating_electrons_mw: float = 0.0
+    """Alpha power heating to electrons [MW]"""
+
     p_fw_alpha_mw: float = 0.0
     """alpha power escaping plasma and reaching first wall (MW)"""
 
     pden_alpha_heating_ions_mw: float = 0.0
     """Alpha power heating per unit volume to ions [MW/m³]"""
 
+    p_alpha_heating_ions_mw: float = 0.0
+    """Alpha power heating to ions [MW]"""
+
     p_plasma_alpha_mw: float = 0.0
     """Alpha power from only the plasma (MW)"""
 
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 3bf3e74c7c..83c9861313 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -676,6 +676,15 @@ def run(self):
             self.data.physics.f_p_alpha_plasma_deposited,
         )
 
+        self.data.physics.p_alpha_heating_electrons_mw = (
+            self.data.physics.pden_alpha_heating_electrons_mw
+            * self.data.physics.vol_plasma
+        )
+
+        self.data.physics.p_alpha_heating_ions_mw = (
+            self.data.physics.pden_alpha_heating_ions_mw * self.data.physics.vol_plasma
+        )
+
         self.data.physics.beta_fast_alpha = self.beta.fast_alpha_beta(
             b_plasma_poloidal_average=self.data.physics.b_plasma_surface_poloidal_average,
             b_plasma_toroidal_on_axis=self.data.physics.b_plasma_toroidal_on_axis,
@@ -2005,6 +2014,13 @@ def outplas(self):
             self.data.physics.pden_alpha_heating_electrons_mw,
             "OP ",
         )
+        po.ovarre(
+            self.outfile,
+            "Alpha power transferred to electrons [MW]",
+            "(p_alpha_heating_electrons_mw)",
+            self.data.physics.p_alpha_heating_electrons_mw,
+            "OP ",
+        )
         po.ovarre(
             self.outfile,
             "Alpha power per unit volume transferred to ions (MW/m³)",
@@ -2012,6 +2028,13 @@ def outplas(self):
             self.data.physics.pden_alpha_heating_ions_mw,
             "OP ",
         )
+        po.ovarre(
+            self.outfile,
+            "Alpha power transferred to ions [MW]",
+            "(p_alpha_heating_ions_mw)",
+            self.data.physics.p_alpha_heating_ions_mw,
+            "OP ",
+        )
 
         po.oblnkl(self.outfile)
         po.ocmmnt(self.outfile, "----------------------------")
@@ -2463,8 +2486,8 @@ def output_plasma_power_balance(self) -> None:
         po.ovarre(
             self.outfile,
             "Alpha heating power to electrons (MW)",
-            "(pden_alpha_heating_electrons_mw)",
-            self.data.physics.pden_alpha_heating_electrons_mw,
+            "(p_alpha_heating_electrons_mw)",
+            self.data.physics.p_alpha_heating_electrons_mw,
             "OP ",
         )
         po.ovarre(
@@ -2493,15 +2516,22 @@ def output_plasma_power_balance(self) -> None:
             self.data.physics.p_ion_electron_equilibration_vol_avg_mw,
             "OP ",
         )
+        po.ovarre(
+            self.outfile,
+            "Plasma radiation power (MW)",
+            "(p_plasma_rad_mw)",
+            self.data.physics.p_plasma_rad_mw,
+            "OP ",
+        )
         po.oblnkl(self.outfile)
         p_electron_balance = (
             self.data.physics.p_plasma_ohmic_mw
-            + self.data.physics.pden_alpha_heating_electrons_mw
+            + self.data.physics.p_alpha_heating_electrons_mw
             + self.data.current_drive.p_hcd_injected_electrons_mw
             - self.data.physics.p_electron_transport_loss_mw
             - self.data.physics.p_ion_electron_equilibration_vol_avg_mw
+            - self.data.physics.p_plasma_rad_mw
         )
-        po.oblnkl(self.outfile)
         po.ocmmnt(self.outfile, "----------------------------")
         po.oblnkl(self.outfile)
         po.ovarre(
@@ -2527,8 +2557,8 @@ def output_plasma_power_balance(self) -> None:
         po.ovarre(
             self.outfile,
             "Alpha heating power to ions (MW)",
-            "(pden_alpha_heating_ions_mw)",
-            self.data.physics.pden_alpha_heating_ions_mw,
+            "(p_alpha_heating_ions_mw)",
+            self.data.physics.p_alpha_heating_ions_mw,
             "OP ",
         )
 
@@ -2550,7 +2580,7 @@ def output_plasma_power_balance(self) -> None:
             "OP ",
         )
         p_ion_balance = (
-            self.data.physics.pden_alpha_heating_ions_mw
+            self.data.physics.p_alpha_heating_ions_mw
             + self.data.current_drive.p_hcd_injected_ions_mw
             - self.data.physics.p_ion_transport_loss_mw
             - self.data.physics.p_ion_electron_equilibration_vol_avg_mw

From 1a5cf57d8135d56f1b3c816dfd814249fd805965 Mon Sep 17 00:00:00 2001
From: mn3981 <christopher.ashe@ukaea.uk>
Date: Mon, 20 Jul 2026 17:01:14 +0100
Subject: [PATCH 37/37] Rename electron-alpha equilibration variables to
 reflect thermal context and update related calculations

---
 process/core/io/plot/summary.py             | 16 ++++++++--------
 process/data_structure/physics_variables.py |  8 ++++----
 process/models/physics/physics.py           | 19 +++++++++++--------
 3 files changed, 23 insertions(+), 20 deletions(-)

diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py
index e49fd5bf97..73bb5fe489 100644
--- a/process/core/io/plot/summary.py
+++ b/process/core/io/plot/summary.py
@@ -14108,10 +14108,10 @@ def plot_plasma_equilibration_time_profile(
         for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)))
     ]
 
-    t_plasma_electron_alpha_equilibration_profile = [
-        mfile_data.data[f"t_plasma_electron_alpha_equilibration_profile{i}"].get_scan(
-            scan
-        )
+    t_plasma_electron_alpha_thermal_equilibration_profile = [
+        mfile_data.data[
+            f"t_plasma_electron_alpha_thermal_equilibration_profile{i}"
+        ].get_scan(scan)
         for i in range(int(mfile_data.data["n_plasma_profile_elements"].get_scan(scan)))
     ]
 
@@ -14132,11 +14132,11 @@ def plot_plasma_equilibration_time_profile(
     )
 
     axis.plot(
-        np.linspace(0, 1, len(t_plasma_electron_alpha_equilibration_profile)),
-        t_plasma_electron_alpha_equilibration_profile,
+        np.linspace(0, 1, len(t_plasma_electron_alpha_thermal_equilibration_profile)),
+        t_plasma_electron_alpha_thermal_equilibration_profile,
         color="red",
         linestyle="-",
-        label=r"$\tau_{e-\alpha,eq}$",
+        label=r"$\tau_{e-\alpha,thermal,eq}$",
     )
 
     axis.set_yscale("log")
@@ -14846,7 +14846,7 @@ def plot_detailed_plasma_parameters(axis: plt.Axes, fig, mfile: MFile, scan: int
         f"$\\mathbf{{Equilibration \\ Times:}}$\n\n"
         f"$\\langle\\tau_{{e-D,Equil}}\\rangle$: {mfile.get('t_plasma_electron_deuteron_equilibration_vol_avg', scan=scan):.4e} s\n"
         f"$\\langle\\tau_{{e-T,Equil}}\\rangle$: {mfile.get('t_plasma_electron_triton_equilibration_vol_avg', scan=scan):.4e} s\n"
-        f"$\\langle\\tau_{{e-\\alpha,Equil}}\\rangle$: {mfile.get('t_plasma_electron_alpha_equilibration_vol_avg', scan=scan):.4e} s"
+        f"$\\langle\\tau_{{e-\\alpha,Equil}}\\rangle$: {mfile.get('t_plasma_electron_alpha_thermal_equilibration_vol_avg', scan=scan):.4e} s"
     )
 
     light_yellow_box = {
diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py
index 35724cca4a..596b821dfc 100644
--- a/process/data_structure/physics_variables.py
+++ b/process/data_structure/physics_variables.py
@@ -1609,13 +1609,13 @@ class PhysicsData:
     )
     """Profile of electron-triton equilibration time in plasma (s)"""
 
-    t_plasma_electron_alpha_equilibration_vol_avg: float = 4.5
-    """Volume averaged electron-alpha equilibration time in plasma (s)"""
+    t_plasma_electron_alpha_thermal_equilibration_vol_avg: float = 4.5
+    """Volume averaged electron-thermal alphas equilibration time in plasma (s)"""
 
-    t_plasma_electron_alpha_equilibration_profile: list[float] = field(
+    t_plasma_electron_alpha_thermal_equilibration_profile: list[float] = field(
         default_factory=list
     )
-    """Profile of electron-alpha equilibration time in plasma (s)"""
+    """Profile of electron-thermal alphas equilibration time in plasma (s)"""
 
     freq_plasma_electron_electron_collision_profile: list[float] = field(
         default_factory=list
diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py
index 83c9861313..f6a4046f0f 100644
--- a/process/models/physics/physics.py
+++ b/process/models/physics/physics.py
@@ -709,7 +709,8 @@ def run(self):
             / (
                 1.0 / self.data.physics.t_plasma_electron_deuteron_equilibration_vol_avg
                 + 1.0 / self.data.physics.t_plasma_electron_triton_equilibration_vol_avg
-                + 1.0 / self.data.physics.t_plasma_electron_alpha_equilibration_vol_avg
+                + 1.0
+                / self.data.physics.t_plasma_electron_alpha_thermal_equilibration_vol_avg
             ),
         )
         self.data.physics.pden_ion_electron_equilibration_vol_avg_mw = (
@@ -5717,7 +5718,7 @@ def run(self):
             n_charge_ion=1,
         )
 
-        self.data.physics.t_plasma_electron_alpha_equilibration_vol_avg = self.calculate_equilibriation_time(
+        self.data.physics.t_plasma_electron_alpha_thermal_equilibration_vol_avg = self.calculate_equilibriation_time(
             temp_plasma_electron_kev=self.data.physics.temp_plasma_electron_vol_avg_kev,
             temp_plasma_ion_kev=self.data.physics.temp_plasma_electron_vol_avg_kev
             * self.data.physics.f_temp_plasma_ion_electron,
@@ -5727,7 +5728,7 @@ def run(self):
             n_charge_ion=2,
         )
 
-        self.data.physics.t_plasma_electron_alpha_equilibration_profile = self.calculate_equilibriation_time(
+        self.data.physics.t_plasma_electron_alpha_thermal_equilibration_profile = self.calculate_equilibriation_time(
             temp_plasma_electron_kev=self.plasma_profile.teprofile.profile_y,
             temp_plasma_ion_kev=self.plasma_profile.teprofile.profile_y
             * self.data.physics.f_temp_plasma_ion_electron,
@@ -6844,16 +6845,18 @@ def output_detailed_physics(self):
         po.ovarre(
             self.outfile,
             "Volume averaged electron-alpha equilibration time (τ_eq) (s)",
-            "(t_plasma_electron_alpha_equilibration_vol_avg)",
-            self.data.physics.t_plasma_electron_alpha_equilibration_vol_avg,
+            "(t_plasma_electron_alpha_thermal_equilibration_vol_avg)",
+            self.data.physics.t_plasma_electron_alpha_thermal_equilibration_vol_avg,
         )
 
         for i in range(
-            len(self.data.physics.t_plasma_electron_alpha_equilibration_profile)
+            len(self.data.physics.t_plasma_electron_alpha_thermal_equilibration_profile)
         ):
             po.ovarre(
                 self.mfile,
                 f"Electron-alpha equilibration time at point {i}",
-                f"(t_plasma_electron_alpha_equilibration_profile{i})",
-                self.data.physics.t_plasma_electron_alpha_equilibration_profile[i],
+                f"(t_plasma_electron_alpha_thermal_equilibration_profile{i})",
+                self.data.physics.t_plasma_electron_alpha_thermal_equilibration_profile[
+                    i
+                ],
             )