diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp
index 3e9a2324dd..6150301033 100644
--- a/src/particles/plasma/PlasmaParticleContainer.cpp
+++ b/src/particles/plasma/PlasmaParticleContainer.cpp
@@ -629,7 +629,7 @@ LaserIonization (const int islice,
 
         // This kernel supports multiple deposition orders (0, 1, 2, 3) at compile time.
         // It calculates the momentum of ionized electrons based on equations (B8) and (B9)
-        // from the Massimo (2020) article and equation (14) from the C. Schroeder (2014) article.
+        // from the F. Massimo (2020) article and equation (12) from the P. Tomassini (2021) article.
         // Additionally, it computes the energy of emitted electrons and assigns their properties
         // (momentum, position, etc.) to newly created electrons in the plasma container.
         amrex::AnyCTO(
@@ -675,8 +675,11 @@ LaserIonization (const int islice,
                     Ep *= E0;
                     amrex::Real delta = std::sqrt(Ep) * laser_dp_prefactor[ion_lev_loc];
                     amrex::Real delta2 = delta * delta;
-                    amrex::Real width_p = amrex::abs(A) * delta * (1._rt - (3._rt/4._rt) * delta2
-                        - (3._rt/2._rt) * delta2 + laser_dp_second_prefactor[ion_lev_loc] * delta2);
+                    amrex::Real delta4 = delta2 * delta2;
+                    amrex::Real alpha = laser_dp_second_prefactor[ion_lev_loc];
+                    amrex::Real s1 = - (7._rt/4._rt) + alpha / 2._rt;
+                    amrex::Real s2 = (1._rt/16._rt) * ( 8._rt * (alpha*alpha) - 68._rt*alpha + 131._rt );
+                    amrex::Real width_p = amrex::abs(A) * delta * (1._rt + s1*delta2 + s2*delta4);
                     amrex::Real p_pol = amrex::RandomNormal(0.0, width_p, engine);
                     ux = p_pol; // linear polarization is assumed along x.
                     uz = (amrex::abs(A * A) * 0.25_rt + p_pol * p_pol * 0.5_rt);
diff --git a/src/particles/plasma/PlasmaParticleContainerInit.cpp b/src/particles/plasma/PlasmaParticleContainerInit.cpp
index 60be0acb8d..54b7140b0a 100644
--- a/src/particles/plasma/PlasmaParticleContainerInit.cpp
+++ b/src/particles/plasma/PlasmaParticleContainerInit.cpp
@@ -453,7 +453,7 @@ InitIonizationModule (const amrex::Geometry& geom, const amrex::Real background_
         h_adk_exp_prefactor[i] = -2./3. * std::pow( Uion/UH,3./2.) * Ea;
         h_laser_adk_prefactor[i] = (3./MathConst::pi) * std::pow(Uion/UH, -3./2.) / Ea;
         h_laser_dp_prefactor[i] = std::sqrt(3./2./Ea) * std::pow(UH/Uion, 3./4.);
-        h_laser_dp_second_prefactor[i] = 2 * ion_atomic_number * std::sqrt(UH/Uion) - 1./2.;
+        h_laser_dp_second_prefactor[i] = 2.*ion_atomic_number * std::sqrt(UH/Uion) - 1.;
     }
 
     amrex::Gpu::copy(amrex::Gpu::hostToDevice,
diff --git a/tests/checksum/benchmarks_json/laser_ionization.1Rank.json b/tests/checksum/benchmarks_json/laser_ionization.1Rank.json
index b32e1976dd..264eb92ece 100644
--- a/tests/checksum/benchmarks_json/laser_ionization.1Rank.json
+++ b/tests/checksum/benchmarks_json/laser_ionization.1Rank.json
@@ -1,23 +1,23 @@
 {
   "lev=0": {
-    "Bx": 13.149192792269,
-    "By": 9.6319466943698,
-    "Bz": 2214.7586533665,
+    "Bx": 13.128118563716,
+    "By": 9.6402246361225,
+    "Bz": 2204.1621631093,
     "ExmBy": 0.0,
     "EypBx": 0.0,
-    "Ez": 779215198153.58,
+    "Ez": 775483620764.12,
     "Psi": 0.0,
-    "Sx": 1114892448780.4,
-    "Sy": 1516523075637.6,
+    "Sx": 1115240346817.3,
+    "Sy": 1513788006780.0,
     "aabs": 0.46877495720552,
-    "chi": 265993671876340.0,
-    "jx": 704455285849490.0,
+    "chi": 265993671876490.0,
+    "jx": 700841379585790.0,
     "jx_beam": 0.0,
     "jy": 0.0,
     "jy_beam": 0.0,
     "jz_beam": 0.0,
     "laserEnvelope": 74.898701643489,
-    "rho_elec": 1202832873.9999,
+    "rho_elec": 1202832836.9171,
     "rho_ion": 1202806069.8845,
     "rhomjz": 0.0
   }