diff --git a/desc/integrals/singularities.py b/desc/integrals/singularities.py
index 110e60f79..1679e8b13 100644
--- a/desc/integrals/singularities.py
+++ b/desc/integrals/singularities.py
@@ -708,7 +708,7 @@ def _kernel_biot_savart_A(eval_data, source_data, diag=False):
 
 
 def _kernel_Bn_over_r(eval_data, source_data, diag=False):
-    # Bn / |r|
+    # B dot n' / |r|
     source_x = jnp.atleast_2d(
         rpz2xyz(jnp.array([source_data["R"], source_data["phi"], source_data["Z"]]).T)
     )
@@ -727,7 +727,7 @@ def _kernel_Bn_over_r(eval_data, source_data, diag=False):
 
 
 def _kernel_Phi_dG_dn(eval_data, source_data, diag=False):
-    # Phi * dG(x,x')/dn' = -Phi * n dot r / r^3 where Phi has units tesla-meters.
+    # Phi * dG(x,x')/dn' = -Phi * n' dot r / r^3 where Phi has units tesla-meters.
     source_x = jnp.atleast_2d(
         rpz2xyz(jnp.array([source_data["R"], source_data["phi"], source_data["Z"]]).T)
     )
@@ -740,7 +740,7 @@ def _kernel_Phi_dG_dn(eval_data, source_data, diag=False):
         dx = eval_x[:, None] - source_x[None]
     n = rpz2xyz_vec(source_data["n_rho"], phi=source_data["phi"])
     return safediv(
-        source_data["Phi"] * dot(n, dx, axis=-1),
+        -source_data["Phi"] * dot(n, dx, axis=-1),
         safenorm(dx, axis=-1) ** 3,
     )
 
@@ -960,7 +960,9 @@ def compute_B_laplace(eq, eval_grid, source_grid=None, B0=None, G=1, R0=1, sym="
     source_data = eq.compute(data_keys, grid=source_grid)
     if B0 is None:
         B0 = ToroidalMagneticField(B0=G, R0=R0)
-    source_data["Bn"], _ = B0.compute_Bnormal(eq.surface, eval_grid=eval_grid)
+    source_data["Bn"], _ = B0.compute_Bnormal(
+        eq.surface, eval_grid=source_grid, source_grid=source_grid
+    )
 
     basis = DoubleFourierSeries(M=eq.M, N=eq.N, NFP=eq.NFP, sym=sym)
     trans_evl = Transform(eval_grid, basis)
@@ -1002,15 +1004,16 @@ def LHS(Phi_mn):
     RHS = -singular_integral(
         eval_data,
         source_data,
-        _kernel_Bn_over_r,
+        _kernel_Bn_over_r,  # flux surface assumption 𝐁₀ = −∇ϕ
         interpolator,
         loop=True,
     ).squeeze() / (2 * jnp.pi)
     # Fourier coefficients of Φ on boundary
     Phi_mn, _, _, _ = jnp.linalg.lstsq(A, RHS)
 
-    # 𝐁 - 𝐁₀ = ∇Φ = 𝐁_vacuum in the interior. Note that Merkel eq. 3.5 has a different
-    # sign on the first integral compared to eq. 1.4... Rederive see which is correct.
+    # 𝐁 - 𝐁₀ = ∇Φ = 𝐁_vacuum in the interior.
+    # Merkel eq. 1.4 is the Green's function solution to ∇²Φ = 0 in the interior.
+    # Note that 𝐁₀′ in eq. 3.5 has the wrong sign.
     grad_Phi = FourierCurrentPotentialField.from_surface(
         eq.surface, Phi_mn / mu_0, basis.modes[:, 1:]
     )