BLAST-WarpX · lucafedeli88 · Nov 27, 2025 · Nov 27, 2025 · Nov 27, 2025 · Nov 27, 2025
diff --git a/Docs/source/refs.bib b/Docs/source/refs.bib
@@ -2791,3 +2791,14 @@ @article{Zhang2017
 numpages = {28},
 pages = {011304},
 }
+
+@article{NarozhnyPRA2004,
+  title={On e+ e- pair production by colliding electromagnetic pulses},
+  author={Narozhny, Nikolay Borisovich and Bulanov, SS and Mur, Vadim Davydovich and Popov, Vladimir Stepanovich},
+  journal={Journal of Experimental and Theoretical Physics Letters},
+  volume={80},
+  number={6},
+  pages={382--385},
+  year={2004},
+  publisher={Springer}
+}
diff --git a/Docs/source/theory/multiphysics/qed.rst b/Docs/source/theory/multiphysics/qed.rst
@@ -20,7 +20,31 @@ Breit-Wheeler
 Schwinger process
 -----------------
 
-If the code is compiled with QED and the user activates the Schwinger process in the input file,
+In the Schwinger process, electron-positron pairs are created in vacuum by a sufficiently strong electromagnetic field.
+An expression for the Schwinger pair production rate can be found in :cite:p:`mqed-NarozhnyPRA2004`.
+Rewritten in SI units, the pair production rate per unit volume is:
+
+.. math::
+
+    \dfrac{d^2N}{dt dV} = \dfrac{e^2 E_s^2}{4 \pi^2 \hbar^2 c} \epsilon \eta \coth{\left({\dfrac{\pi \eta}{\epsilon}}\right)}\exp{\left({-\dfrac{\pi}{\epsilon}}\right)}
+
+
+where $e$ is the elementary charge, :math:`E_s`` is the Schwinger field, :math:`\epsilon = \mathcal{E}/E_s` and :math:`\eta = \mathcal{H}/E_s`.
+:math:`\mathcal{E}` and :math:`\mathcal{H}$` are given by:
+
+.. math::
+
+    \mathcal{E} = \sqrt{\sqrt{\mathcal{F}^2 + \mathcal{G}^2} + \mathcal{F}} \\
+    \mathcal{H} = \sqrt{\sqrt{\mathcal{F}^2 + \mathcal{G}^2} - \mathcal{F}} \\
+
+:math:`\mathcal{F}` and :math:`\mathcal{G}$` are the invariants of the electromagnetic field and are equal to:
+
+.. math::
+
+    \mathcal{F} = (\mathbf{E}^2 - c^2 \mathbf{B}^2)/2 \\
+    \mathcal{G} = c \mathbf{E} \cdot \mathbf{B} \\
+
+If the user activates the Schwinger process in the input file,
 electron-positron pairs can be created in vacuum in the function
 ``MultiParticleContainer::doQEDSchwinger``:
 
@@ -44,3 +68,6 @@ the total number of physical pairs created.
 So far the Schwinger module requires using ``warpx.grid_type = collocated`` or
 ``algo.field_gathering = momentum-conserving`` (so that the auxiliary fields are calculated on the nodes)
 and is not compatible with either mesh refinement, RZ, RCYLINDER, and RSPHERE coordinates or single precision.
+
+.. bibliography::
+    :keyprefix: mqed-
diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
@@ -3902,7 +3902,7 @@ Lookup tables and other settings for QED modules
 ------------------------------------------------
 
 Lookup tables store pre-computed values for functions used by the QED modules.
-**This feature requires to compile with QED=TRUE (and also with QED_TABLE_GEN=TRUE for table generation)**
+Unless otherwise specified, **these feature requires to compile with QED=TRUE (and also with QED_TABLE_GEN=TRUE for table generation)**
 
 * ``qed_bw.lookup_table_mode`` (`string`)
     There are three options to prepare the lookup table required by the Breit-Wheeler module:
@@ -3991,7 +3991,7 @@ Lookup tables store pre-computed values for functions used by the QED modules.
 
 * ``warpx.do_qed_schwinger`` (`bool`) optional (default `0`)
     If this is 1, Schwinger electron-positron pairs can be generated in vacuum in the cells where the EM field is high enough.
-    Activating the Schwinger process requires the code to be compiled with ``QED=TRUE`` and ``PICSAR``.
+    Activating the Schwinger process does not require the code to be compiled with ``QED=TRUE``.
     If ``warpx.do_qed_schwinger = 1``, Schwinger product species must be specified with
     ``qed_schwinger.ele_product_species`` and ``qed_schwinger.pos_product_species``.
     Schwinger process requires either ``warpx.grid_type = collocated`` or

diff --git a/Regression/Checksum/benchmarks_json/test_3d_qed_schwinger_4.json b/Regression/Checksum/benchmarks_json/test_3d_qed_schwinger_4.json
@@ -1,31 +1,31 @@
 {
-  "lev=0": {
-    "Bx": 10.282332458788613,
-    "By": 1707847966720000.0,
-    "Bz": 0.5093466319155224,
-    "Ex": 3460290602.1197195,
-    "Ey": 0.0,
-    "Ez": 5.119999999999012e+23,
-    "jx": 354749197616543.3,
-    "jy": 0.0,
-    "jz": 8899048833947625.0
-  },
   "ele_schwinger": {
     "particle_momentum_x": 3.695301335661867e-12,
-    "particle_momentum_y": 1.0868696459770461e-26,
+    "particle_momentum_y": 1.0868696454896924e-26,
     "particle_momentum_z": 2.8810249733308265e-12,
     "particle_position_x": 0.00028800000000000006,
     "particle_position_y": 0.000216,
     "particle_position_z": 0.0002540000004700524,
-    "particle_weight": 63968.698937700814
+    "particle_weight": 63968.69690641211
+  },
+  "lev=0": {
+    "Bx": 10.282332458788613,
+    "By": 1707847966720000.0,
+    "Bz": 0.5093466195466327,
+    "Ex": 3460290503.614944,
+    "Ey": 0.0,
+    "Ez": 5.119999999999012e+23,
+    "jx": 354749188059060.5,
+    "jy": 0.0,
+    "jz": 8899048551363524.0
   },
   "pos_schwinger": {
     "particle_momentum_x": 3.695301335661867e-12,
-    "particle_momentum_y": 1.0881422563719172e-26,
+    "particle_momentum_y": 1.0881422558277477e-26,
     "particle_momentum_z": 2.881024973330826e-12,
     "particle_position_x": 0.000288,
     "particle_position_y": 0.000216,
     "particle_position_z": 0.00023799999952994755,
-    "particle_weight": 63968.698937700814
+    "particle_weight": 63968.696906412115
   }
 }
diff --git a/Source/Diagnostics/ReducedDiags/ColliderRelevant.cpp b/Source/Diagnostics/ReducedDiags/ColliderRelevant.cpp
@@ -9,9 +9,6 @@
 
 #include "Diagnostics/ReducedDiags/ReducedDiags.H"
 #include "Fields.H"
-#if (defined WARPX_QED)
-#   include "Particles/ElementaryProcess/QEDInternals/QedChiFunctions.H"
-#endif
 #include "Particles/Gather/FieldGather.H"
 #include "Particles/Gather/GetExternalFields.H"
 #include "Particles/MultiParticleContainer.H"
@@ -47,6 +44,7 @@
 #include <AMReX_Vector.H>
 
 #include <ablastr/coarsen/sample.H>
+#include <ablastr/qed/ChiFunctions.H>
 #include <ablastr/warn_manager/WarnManager.H>
 
 #include <algorithm>
@@ -521,10 +519,10 @@ void ColliderRelevant::ComputeDiags (int step)
                     // compute chi
                     amrex::Real chi = 0.0_rt;
                     if (is_photon) {
-                        chi = QedUtils::chi_photon(ux[i]*m, uy[i]*m, uz[i]*m,
+                        chi = ablastr::qed::chi_photon(ux[i]*m, uy[i]*m, uz[i]*m,
                                             ex, ey, ez, bx, by, bz);
                     } else {
-                        chi = QedUtils::chi_ele_pos(ux[i]*m, uy[i]*m, uz[i]*m,
+                        chi = ablastr::qed::chi_ele_pos(ux[i]*m, uy[i]*m, uz[i]*m,
                                             ex, ey, ez, bx, by, bz);
                     }
                     return {chi, chi, chi*w[i]};

diff --git a/Source/Diagnostics/ReducedDiags/ParticleExtrema.cpp b/Source/Diagnostics/ReducedDiags/ParticleExtrema.cpp
@@ -8,9 +8,6 @@
 #include "ParticleExtrema.H"
 
 #include "Diagnostics/ReducedDiags/ReducedDiags.H"
-#if (defined WARPX_QED)
-#   include "Particles/ElementaryProcess/QEDInternals/QedChiFunctions.H"
-#endif
 #include "Fields.H"
 #include "Particles/Gather/FieldGather.H"
 #include "Particles/Gather/GetExternalFields.H"
@@ -22,6 +19,7 @@
 #include "WarpX.H"
 
 #include <ablastr/fields/MultiFabRegister.H>
+#include <ablastr/qed/ChiFunctions.H>
 
 #include <AMReX_Algorithm.H>
 #include <AMReX_Array.H>
@@ -345,10 +343,10 @@ void ParticleExtrema::ComputeDiags (int step)
                         // compute chi
                         amrex::Real chi = 0.0_rt;
                         if ( is_photon ) {
-                            chi = QedUtils::chi_photon(ux[i]*m, uy[i]*m, uz[i]*m,
+                            chi = ablastr::qed::chi_photon(ux[i]*m, uy[i]*m, uz[i]*m,
                                              ex, ey, ez, bx, by, bz);
                         } else {
-                            chi = QedUtils::chi_ele_pos(ux[i]*m, uy[i]*m, uz[i]*m,
+                            chi = ablastr::qed::chi_ele_pos(ux[i]*m, uy[i]*m, uz[i]*m,
                                              ex, ey, ez, bx, by, bz);
                         }
                         return {chi,chi};

diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp
@@ -218,11 +218,12 @@ WarpX::Evolve (int numsteps)
         // multi-physics: field ionization
         doFieldIonization();
 
-#ifdef WARPX_QED
         // multi-physics: QED effects
+#ifdef WARPX_QED
         doQEDEvents();
-        mypc->doQEDSchwinger();
 #endif
+        mypc->doQEDSchwinger();
+
 
         // perform particle injection
         ExecutePythonCallback("particleinjection");

diff --git a/Source/Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper.H b/Source/Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper.H
@@ -9,7 +9,6 @@
 
 #include "BreitWheelerEngineWrapper_fwd.H"
 
-#include "QedChiFunctions.H"
 #include "QedWrapperCommons.H"
 #include "Utils/WarpXConst.H"
 
@@ -18,6 +17,8 @@
 #include <AMReX_REAL.H>
 #include <AMReX_Random.H>
 
+#include <ablastr/qed/ChiFunctions.H>
+
 #include <picsar_qed/containers/picsar_array.hpp>
 #include <picsar_qed/math/cmath_overloads.hpp>
 #include <picsar_qed/math/math_constants.h>
@@ -160,7 +161,7 @@ public:
         const auto py = m_e*uy;
         const auto pz = m_e*uz;
 
-        const auto chi_phot = QedUtils::chi_photon(
+        const auto chi_phot = ablastr::qed::chi_photon(
             px, py, pz, ex, ey, ez, bx, by, bz);
 
         //Optical depth is not evolved for photons having less energy than what is
@@ -255,7 +256,7 @@ public:
         auto py = uy*me;
         auto pz = uz*me;
 
-        const auto chi_photon = QedUtils::chi_photon(
+        const auto chi_photon = ablastr::qed::chi_photon(
             px, py, pz, ex, ey, ez, bx, by, bz);
 
         const auto momentum_photon = pxr_m::vec3<amrex::ParticleReal>{px, py, pz};

diff --git a/Source/Particles/ElementaryProcess/QEDInternals/QedChiFunctions.H b/Source/Particles/ElementaryProcess/QEDInternals/QedChiFunctions.H
diff --git a/Source/Particles/ElementaryProcess/QEDInternals/QuantumSyncEngineWrapper.H b/Source/Particles/ElementaryProcess/QEDInternals/QuantumSyncEngineWrapper.H
@@ -9,7 +9,6 @@
 
 #include "BreitWheelerEngineWrapper_fwd.H"
 
-#include "QedChiFunctions.H"
 #include "QedWrapperCommons.H"
 #include "Utils/WarpXConst.H"
 
@@ -18,6 +17,8 @@
 #include <AMReX_REAL.H>
 #include <AMReX_Random.H>
 
+#include <ablastr/qed/ChiFunctions.H>
+
 #include <picsar_qed/containers/picsar_array.hpp>
 #include <picsar_qed/math/cmath_overloads.hpp>
 #include <picsar_qed/math/vec_functions.hpp>
@@ -150,7 +151,7 @@ public:
         const amrex::ParticleReal gamma = std::sqrt(1._rt + (ux*ux + uy*uy + uz*uz)*inv_c2);
         const auto energy = gamma*m_e*PhysConst::c*PhysConst::c;
 
-        const auto chi_part = QedUtils::chi_ele_pos(
+        const auto chi_part = ablastr::qed::chi_ele_pos(
             m_e*ux, m_e*uy, m_e*uz, ex, ey, ez, bx, by, bz);
 
         if (chi_part < m_qs_minimum_chi_part) {
@@ -234,7 +235,7 @@ public:
         auto py = uy*me;
         auto pz = uz*me;
 
-        const auto chi_particle = QedUtils::chi_ele_pos(
+        const auto chi_particle = ablastr::qed::chi_ele_pos(
             px, py, pz, ex, ey, ez, bx, by, bz);
 
         auto momentum_particle = pxr_m::vec3<amrex::ParticleReal>{px, py, pz};