diff --git a/Regression/Checksum/benchmarks_json/test_rz_magnetostatic_eb_picmi.json b/Regression/Checksum/benchmarks_json/test_rz_magnetostatic_eb_picmi.json
index 5d55368cd9e..d73ddc70731 100644
--- a/Regression/Checksum/benchmarks_json/test_rz_magnetostatic_eb_picmi.json
+++ b/Regression/Checksum/benchmarks_json/test_rz_magnetostatic_eb_picmi.json
@@ -1,10 +1,10 @@
 {
   "lev=0": {
-    "Az": 1.7305570206582026,
-    "Bt": 15.521136797767891,
-    "Er": 4333723925.211036,
+    "Az": 1.7305570206581988,
+    "Bt": 15.521136535283418,
+    "Er": 4333723925.211021,
     "jz": 207219733.49323815,
-    "phi": 479458088.8864637,
+    "phi": 479458088.8864618,
     "rho": 0.6946928324123767
   }
 }
\ No newline at end of file
diff --git a/Source/ablastr/fields/VectorPoissonSolver.H b/Source/ablastr/fields/VectorPoissonSolver.H
index 26e2a0cfb3c..3f32bca51ae 100644
--- a/Source/ablastr/fields/VectorPoissonSolver.H
+++ b/Source/ablastr/fields/VectorPoissonSolver.H
@@ -197,6 +197,12 @@ computeVectorPotential   (  amrex::Vector<amrex::Array<amrex::MultiFab*, 3> > co
 
 #ifdef WARPX_DIM_RZ
             linop[adim]->setRZ(true);
+            if (adim < 2) {
+                // In cylindrical coordinates (with azimuthal symmetry), the vector Laplacian
+                // along r and theta has an extra `-1/r^2` term (not present for Az).
+                // The alpha term below adds the 1/r^2 part.
+                linop[adim]->setAlpha(1._rt);
+            }
 #endif
 
             linop[adim]->setDomainBC( boundary_handler.lobc[adim], boundary_handler.hibc[adim] );