diff --git a/include/picongpu/particles/atomicPhysics/kernel/UpdateElectricField.kernel b/include/picongpu/particles/atomicPhysics/kernel/UpdateElectricField.kernel
new file mode 100644
index 0000000000..99ddffd57c
--- /dev/null
+++ b/include/picongpu/particles/atomicPhysics/kernel/UpdateElectricField.kernel
@@ -0,0 +1,130 @@
+/* Copyright 2024 Brian Marre
+ *
+ * This file is part of PIConGPU.
+ *
+ * PIConGPU is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * PIConGPU is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with PIConGPU.
+ * If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include "picongpu/defines.hpp"
+#include "picongpu/particles/atomicPhysics/KernelIndexation.hpp"
+
+#include <pmacc/lockstep/ForEach.hpp>
+
+#include <cstdint>
+
+namespace picongpu::particles::atomicPhysics::kernel
+{
+    /** update electric field due to field ionization energy use
+     *
+     * @tparam T_numberAtomicPhysicsIonSpecies specialization template parameter used to prevent compilation of all
+     *  atomicPhysics kernels if no atomic physics species is present.
+     */
+    template<uint32_t T_numberAtomicPhysicsIonSpecies>
+    struct UpdateElectricFieldKernel
+    {
+        /** call operator
+         *
+         * called by UpdateElectricField atomic physics sub-stage
+         *
+         * @param worker object containing the device and block information, passed by PMACC_KERNEL call
+         * @param areMapping mapping of blockIndex to block superCell index
+         * @param timeRemainingBox deviceDataBox giving access to the time remaining in the atomicPhysics step for each
+         *  local superCell
+         * @param eFieldBox deviceDataBox giving access to eField values for all local superCells
+         * @param fieldEnergyUseCacheBox deviceDataBox giving access to the field energy use cache for each local
+         *  superCell
+         */
+        template<
+            typename T_Worker,
+            typename T_AreaMapping,
+            typename T_TimeRemainingBox,
+            typename T_EFieldBox,
+            typename T_FieldEnergyUseCacheBox>
+        HDINLINE void operator()(
+            T_Worker const& worker,
+            T_AreaMapping const areaMapping,
+            T_TimeRemainingBox const timeRemainingBox,
+            T_EFieldBox eFieldBox,
+            T_FieldEnergyUseCacheBox const fieldEnergyUseCacheBox) const
+        {
+            auto const superCellIdx = KernelIndexation::getSuperCellIndex(worker, areaMapping);
+            auto const superCellFieldIdx = KernelIndexation::getSuperCellFieldIndex(worker, areaMapping, superCellIdx);
+
+            auto const timeRemaining = timeRemainingBox(superCellFieldIdx);
+            if(timeRemaining <= 0._X)
+                return;
+
+            T_FieldEnergyUseCache const& eFieldEnergyUseCache = fieldEnergyUseCacheBox(superCellFieldIdx);
+            DataSpace<picongpu::simDim> const superCellCellOffset = superCellIdx * picongpu::SuperCellSize::toRT();
+
+            auto forEachCell = pmacc::lockstep::makeForEach<T_FieldEnergyUseCache::numberCells, T_Worker>(worker);
+            forEachCell(
+                [&worker,
+                 &superCellCellOffset,
+                 &eFieldBox,
+                 &eFieldEnergyUseCache,
+                 &sharedResourcesOverSubscribed,
+                 &rejectionProbabilityCache](uint32_t const linearCellIndex)
+                {
+                    DataSpace<picongpu::simDim> const localCellIndex
+                        = pmacc::math::mapToND(picongpu::SuperCellSize::toRT(), static_cast<int>(linearCellIndex));
+                    DataSpace<picongpu::simDim> const cellIndex = localCellIndex + superCellCellOffset;
+
+                    // ((sim.unit.mass() * sim.unit.length())/(sim.unit.time()^2 * sim.unit.charge()))^2
+                    auto const oldEFieldVector = eFieldBox(cellIndex);
+                    auto const eFieldNormSquared = pmacc::math::l2norm2(oldEFieldVector);
+
+                    /// @note zero field also means at most zero energy use
+                    if(eFieldNormSquared <= 0)
+                        return;
+
+                    // eV * 1 = eV * sim.unit.energy()/sim.unit.energy() = (ev / sim.unit.energy()) * sim.unit.energy()
+                    // sim.unit.energy()
+                    float_X const eFieldEnergyUse
+                        = picongpu::sim.pic.get_eV() * eFieldEnergyUseCache.energyUsed(linearCellIndex);
+
+                    // sim.unit.charge()^2 * sim.unit.time()^2 / (sim.unit.mass() * sim.unit.length()^3)
+                    //  * sim.unit.length()^3
+                    // = sim.unit.charge()^2 * sim.unit.time()^2 / sim.unit.mass()
+                    constexpr float_X eps0HalfTimesCellVolume
+                        = (picongpu::sim.pic.getEps0() / 2._X) * picongpu::sim.pic.getCellSize().productOfComponents();
+
+                    /* (((sim.unit.mass() * sim.unit.length())/(sim.unit.time()^2 * sim.unit.charge()))^2
+                     *  - sim.unit.energy() / (sim.unit.charge()^2 * sim.unit.time()^2 / sim.unit.mass()))^1/2
+                     * = (((sim.unit.mass() * sim.unit.length())/(sim.unit.time()^2 * sim.unit.charge()))^2
+                     *  - (sim.unit.mass() * sim.unit.length()^2/sim.unit.time()^2
+                     *      / (sim.unit.charge()^2 * sim.unit.time()^2 / sim.unit.mass()))
+                     * )^1/2
+                     * = (((sim.unit.mass() * sim.unit.length())/(sim.unit.time()^2 * sim.unit.charge()))^2
+                     *  - (sim.unit.mass()^2 * sim.unit.length()^2/(sim.unit.time()^4 * sim.unit.charge()^2))
+                     * )^1/2
+                     * = ((sim.unit.mass()^2 * sim.unit.length()^2)/(sim.unit.time()^4 * sim.unit.charge()^2)
+                     *  - (sim.unit.mass()^2 * sim.unit.length()^2/(sim.unit.time()^4 * sim.unit.charge()^2))
+                     * )^1/2
+                     * = ((sim.unit.mass()^2 * sim.unit.length()^2)/(sim.unit.time()^4 * sim.unit.charge()^2))^1/2
+                     * = (sim.unit.mass() * sim.unit.length())/(sim.unit.time()^2 * sim.unit.charge())
+                     * = sim.unit.eField()
+                     */
+                    // sim.unit.eField()
+                    float_X const newEFieldNorm
+                        = pmacc::math::sqrt(eFieldNormSquared - eFieldEnergyUse / eps0HalfTimesCellVolume);
+
+                    eFieldBox(cellIndex) = oldEFieldVector / pmacc::math::sqrt(eFieldNormSquared) * newEFieldNorm;
+                });
+        }
+    };
+} // namespace picongpu::particles::atomicPhysics::kernel
diff --git a/include/picongpu/particles/atomicPhysics/stage/CheckForOverSubscription.hpp b/include/picongpu/particles/atomicPhysics/stage/CheckForOverSubscription.hpp
index cf708f5202..54323bb524 100644
--- a/include/picongpu/particles/atomicPhysics/stage/CheckForOverSubscription.hpp
+++ b/include/picongpu/particles/atomicPhysics/stage/CheckForOverSubscription.hpp
@@ -26,6 +26,7 @@
 #include "picongpu/fields/FieldE.hpp"
 #include "picongpu/particles/atomicPhysics/electronDistribution/LocalHistogramField.hpp"
 #include "picongpu/particles/atomicPhysics/kernel/CheckForOverSubscription.kernel"
+#include "picongpu/particles/atomicPhysics/localHelperFields/FieldEnergyUseCacheField.hpp"
 #include "picongpu/particles/atomicPhysics/localHelperFields/RejectionProbabilityCacheField.hpp"
 #include "picongpu/particles/atomicPhysics/localHelperFields/SharedResourcesOverSubscribedField.hpp"
 #include "picongpu/particles/atomicPhysics/localHelperFields/TimeRemainingField.hpp"
diff --git a/include/picongpu/particles/atomicPhysics/stage/UpdateElectricField.hpp b/include/picongpu/particles/atomicPhysics/stage/UpdateElectricField.hpp
new file mode 100644
index 0000000000..86c33084f3
--- /dev/null
+++ b/include/picongpu/particles/atomicPhysics/stage/UpdateElectricField.hpp
@@ -0,0 +1,78 @@
+/* Copyright 2024 Brian Marre
+ *
+ * This file is part of PIConGPU.
+ *
+ * PIConGPU is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * PIConGPU is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with PIConGPU.
+ * If not, see <http://www.gnu.org/licenses/>.
+ */
+
+//! @file update electric field with field ionization energy use
+
+#pragma once
+
+// need picongpu::atomicPhysics::ElectronHistogram from atomicPhysics.param
+#include "picongpu/defines.hpp"
+#include "picongpu/fields/FieldE.hpp"
+#include "picongpu/particles/atomicPhysics/kernel/UpdateElectricField.kernel"
+#include "picongpu/particles/atomicPhysics/localHelperFields/FieldEnergyUseCacheField.hpp"
+#include "picongpu/particles/atomicPhysics/localHelperFields/TimeRemainingField.hpp"
+
+#include <pmacc/Environment.hpp>
+#include <pmacc/mappings/kernel/AreaMapping.hpp>
+#include <pmacc/type/Area.hpp>
+
+#include <cstdint>
+
+namespace picongpu::particles::atomicPhysics::stage
+{
+    /** UpdateElectricField atomic physics sub-stage
+     *
+     * remove the used field energy from the electric field by reducing the norm of the local E-Field vector to match
+     *  the energy use
+     *
+     * @tparam T_numberAtomicPhysicsIonSpecies specialization template parameter used to prevent compilation of all
+     *  atomicPhysics kernels if no atomic physics species is present.
+     */
+    template<uint32_t T_numberAtomicPhysicsIonSpecies>
+    struct UpdateElectricField
+    {
+        //! call of kernel for every superCell
+        HINLINE void operator()(picongpu::MappingDesc const mappingDesc) const
+        {
+            // full local domain, no guards
+            pmacc::AreaMapping<CORE + BORDER, MappingDesc> mapper(mappingDesc);
+            pmacc::DataConnector& dc = pmacc::Environment<>::get().DataConnector();
+
+            auto& timeRemainingField = *dc.get<
+                picongpu::particles::atomicPhysics::localHelperFields::TimeRemainingField<picongpu::MappingDesc>>(
+                "TimeRemainingField");
+
+            auto& eField = *dc.get<FieldE>(FieldE::getName());
+
+            using FieldEnergyUseCacheField
+                = picongpu::particles::atomicPhysics::localHelperFields::FieldEnergyUseCacheField<
+                    picongpu::MappingDesc>;
+            auto& fieldEnergyUseCacheField = *dc.get<FieldEnergyUseCacheField>("FieldEnergyUseCacheField");
+
+            // macro for call of kernel for every superCell, see pull request #4321
+            PMACC_LOCKSTEP_KERNEL(
+                picongpu::particles::atomicPhysics::kernel::UpdateElectricField<T_numberAtomicPhysicsIonSpecies>())
+                .template config<FieldEnergyUseCacheField::ElementType::numberCells>(mapper.getGridDim())(
+                    mapper,
+                    timeRemainingField.getDeviceDataBox(),
+                    eField.getDeviceDataBox(),
+                    fieldEnergyUseCacheField.getDeviceDataBox());
+        }
+    };
+} // namespace picongpu::particles::atomicPhysics::stage
diff --git a/include/picongpu/simulation/stage/AtomicPhysics.x.cpp b/include/picongpu/simulation/stage/AtomicPhysics.x.cpp
index e0ba1c2cc2..f48266a811 100644
--- a/include/picongpu/simulation/stage/AtomicPhysics.x.cpp
+++ b/include/picongpu/simulation/stage/AtomicPhysics.x.cpp
@@ -46,6 +46,7 @@
 #include "picongpu/particles/atomicPhysics/stage/ResetTimeStepField.hpp"
 #include "picongpu/particles/atomicPhysics/stage/RollForOverSubscription.hpp"
 #include "picongpu/particles/atomicPhysics/stage/SpawnIonizationElectrons.hpp"
+#include "picongpu/particles/atomicPhysics/stage/UpdateElectricField.hpp"
 #include "picongpu/particles/atomicPhysics/stage/UpdateTimeRemaining.hpp"
 #include "picongpu/particles/filter/filter.hpp"
 
@@ -440,6 +441,12 @@ namespace picongpu::simulation::stage
                 ForEachIonSpeciesSpawnIonizationElectrons{}(mappingDesc, currentStep);
             }
 
+            HINLINE static void updateElectricField(picongpu::MappingDesc const& mappingDesc)
+            {
+                picongpu::particles::atomicPhysics::stage::UpdateElectricField<T_numberAtomicPhysicsIonSpecies>()(
+                    mappingDesc);
+            }
+
             template<typename T_DeviceReduce>
             HINLINE static void doIPDIonization(
                 picongpu::MappingDesc const& mappingDesc,
@@ -576,6 +583,7 @@ namespace picongpu::simulation::stage
 
                     recordChanges(mappingDesc);
                     updateElectrons(mappingDesc, currentStep);
+                    updateElectricField(mappingDesc);
                     doIPDIonization(mappingDesc, currentStep, deviceLocalReduce);
                     updateTimeRemaining(mappingDesc);
                     isSubSteppingComplete = isSubSteppingFinished(mappingDesc, deviceLocalReduce);