diff --git a/src/struphy/propagators/gyrokinetic_poisson_solve.py b/src/struphy/propagators/gyrokinetic_poisson_solve.py
index f279d4ef7..b5bc5a89e 100644
--- a/src/struphy/propagators/gyrokinetic_poisson_solve.py
+++ b/src/struphy/propagators/gyrokinetic_poisson_solve.py
@@ -180,6 +180,8 @@ class Options:
         multigrid: MultiGridOptions = None
         solver_params: SolverParameters = None
         filter_params: dict[PICVariable | SPHVariable, FilterParameters] = None
+        # the adiabatic term is physical, not a stabilization: the rhs is not projected
+        enforce_compatibility: bool = False
 
         def __post_init__(self):
             # checks
diff --git a/src/struphy/propagators/implicit_diffusion.py b/src/struphy/propagators/implicit_diffusion.py
index 86f23feb4..fff7b4e9d 100644
--- a/src/struphy/propagators/implicit_diffusion.py
+++ b/src/struphy/propagators/implicit_diffusion.py
@@ -135,7 +135,7 @@ def __init__(
 
     @dataclass(repr=False)
     class Options(OptionsBase):
-        """Configuration options for :class:`ImplicitDiffusion`.
+        r"""Configuration options for :class:`ImplicitDiffusion`.
 
         Parameters
         ----------
@@ -222,6 +222,22 @@ class Options(OptionsBase):
             ``verbose``, ``info``, ``recycle``).
             If ``None``, defaults to ``SolverParameters()``.
 
+        enforce_compatibility : bool, default=False
+            If ``True`` and the diffusion operator is singular (no Dirichlet boundary conditions
+            and no polar splines, i.e. the constant functions are in its kernel), the right-hand side
+            is made compatible before the solve by removing its constant ("net charge") component,
+
+            .. math::
+
+                \mathbf b \leftarrow \mathbf b - \frac{\mathbf 1^\top \mathbf b}{\mathbf 1^\top \mathbb M^0 \mathbf 1}\, \mathbb M^0 \mathbf 1\,,
+
+            where :math:`\mathbf 1` (all ones) is the coefficient vector of the constant function
+            (partition of unity). This is equivalent to adding a uniform neutralizing background.
+            It is meant for (nearly) singular problems where ``sigma_1`` is only a small
+            stabilization; then ``phi`` stays of order one instead of growing like
+            ``(net charge) / sigma_1``, which would destroy the accuracy of ``grad phi``.
+            For ``stab_mat="M0"`` and ``sigma_2 = 0`` it changes ``phi`` only by a constant, for any ``sigma_1``.
+
         filter_params : dict[PICVariable | SPHVariable, FilterParameters], default=None
             If not None, specifies a filter to the accumulation of a specific variable.
             Keyed by the ``pic_variable`` of the corresponding
@@ -251,6 +267,7 @@ class Options(OptionsBase):
         precond_params: dict = None
         multigrid: MultiGridOptions = None
         solver_params: SolverParameters = None
+        enforce_compatibility: bool = False
         filter_params: dict[PICVariable | SPHVariable, FilterParameters] = None
 
         def __post_init__(self):
@@ -435,6 +452,28 @@ def verify_rhs(rho) -> StencilVector | FEECVariable | AccumulatorVector:
         self._rhs2 = phi.space.zeros()
         self._tmp_src = phi.space.zeros()
 
+        # constant function and its weak form M0 @ 1, for removing the kernel component of the rhs
+        self._ones = None
+        if self.options.enforce_compatibility and self.operator_is_singular:
+            self._ones = phi.space.zeros()
+            owned = tuple(slice(s, e + 1) for s, e in zip(self._ones.starts, self._ones.ends))
+            self._ones[owned] = 1.0
+            self._ones.update_ghost_regions()
+            self._m0_ones = self.mass_ops.M0.dot(self._ones)
+            self._ones_m0_ones = self._ones.inner(self._m0_ones)
+
+    @property
+    def operator_is_singular(self) -> bool:
+        """True if the diffusion operator ``grad.T @ diffusion_mat @ grad`` is known to be singular
+        with the constant functions (coefficient vector of all ones) in its kernel: no Dirichlet boundary
+        conditions and tensor-product splines. Polar splines are conservatively treated as regular."""
+        no_dirichlet = not any(any(bc) for bc in self.derham.dirichlet_bc)
+        return (
+            no_dirichlet
+            and not self.derham.polar_splines
+            and isinstance(self.variables.phi.spline.vector, StencilVector)
+        )
+
     def _build_stiffness_pc(self, sig_1: float) -> preconditioner.StiffnessPreconditioner:
         """StiffnessPreconditioner for ``sig_1 * stab_mat (I - A) + grad.T @ diffusion_mat @ grad``."""
         return preconditioner.StiffnessPreconditioner(
@@ -531,6 +570,10 @@ def __call__(self, dt):
 
         rhs += self._rhs2
 
+        # remove the kernel component (net charge) of the rhs: 1^T rhs = 0 (compatibility condition)
+        if self._ones is not None:
+            rhs.mul_iadd(-self._ones.inner(rhs) / self._ones_m0_ones, self._m0_ones)
+
         if self.diagnostic is not None:
             proj = L2Projector("H1", self.mass_ops)
             self.diagnostic.spline.vector = proj.solve(rhs)
diff --git a/src/struphy/propagators/poisson_solve.py b/src/struphy/propagators/poisson_solve.py
index d23253045..79137400c 100644
--- a/src/struphy/propagators/poisson_solve.py
+++ b/src/struphy/propagators/poisson_solve.py
@@ -36,7 +36,7 @@ class PoissonSolve(ImplicitDiffusion):
 
     @dataclass(repr=False)
     class Options(OptionsBase):
-        """Configuration options for :class:`Poisson`.
+        r"""Configuration options for :class:`Poisson`.
 
         Parameters
         ----------
@@ -93,6 +93,23 @@ class Options(OptionsBase):
             ``verbose``, ``info``, ``recycle``).
             If ``None``, defaults to ``SolverParameters()``.
 
+        enforce_compatibility : bool, default=True
+            If ``True`` and the Poisson operator is singular (no Dirichlet boundary conditions
+            and no polar splines, e.g. fully periodic), the constant ("net charge") component of the
+            right-hand side is removed before the solve,
+
+            .. math::
+
+                \mathbf b \leftarrow \mathbf b - \frac{\mathbf 1^\top \mathbf b}{\mathbf 1^\top \mathbb M^0 \mathbf 1}\, \mathbb M^0 \mathbf 1\,,
+
+            such that the discrete compatibility condition :math:`\mathbf 1^\top \mathbf b = 0` holds
+            (:math:`\mathbf 1` is the coefficient vector of the constant function). This is equivalent to
+            a uniform neutralizing background; for particle sources, the net charge is Monte-Carlo noise.
+            Without it, the constant mode of ``phi`` becomes ``(net charge) / stab_eps``, and
+            ``grad phi`` loses digits to cancellation. ``stab_eps`` then only fixes the constant of ``phi``.
+            Set to ``False`` if ``stab_eps * stab_mat`` is a physical term (e.g. adiabatic electrons with
+            ``stab_mat="M0ad"``) and not just a stabilization.
+
         filter_params : dict[PICVariable | SPHVariable, FilterParameters], default=None
             If not None, specifies a filter to the accumulation of a specific variable.
 
@@ -118,6 +135,7 @@ class Options(OptionsBase):
         precond_params: dict = None
         multigrid: MultiGridOptions = None
         solver_params: SolverParameters = None
+        enforce_compatibility: bool = True
         filter_params: dict[PICVariable | SPHVariable, FilterParameters] = None
 
         def __post_init__(self):
diff --git a/src/struphy/propagators/tests/test_poisson.py b/src/struphy/propagators/tests/test_poisson.py
index 9c15ee437..b686ed354 100644
--- a/src/struphy/propagators/tests/test_poisson.py
+++ b/src/struphy/propagators/tests/test_poisson.py
@@ -880,6 +880,99 @@ def test_implicit_diffusion_multigrid_dt():
         assert prop._solver._info["niter"] <= 15
 
 
+@pytest.mark.parametrize("stab_mat", ["Id", "M0"])
+def test_poisson_net_charge_periodic_1d(stab_mat):
+    """Periodic Poisson problem with a non-neutral source (net charge, e.g. Monte-Carlo noise).
+
+    The constant functions span the kernel of the periodic Laplacian. With the default options the
+    net charge is removed from the rhs (compatibility condition), such that phi stays O(1) and
+    E = -grad phi is accurate. Without it, the constant mode of phi is (net charge) / stab_eps.
+    """
+    Lx = 12.56
+    k = 2 * xp.pi / Lx
+    amp = 0.1
+    net = 0.3
+
+    domain = domains.Cuboid(r1=Lx)
+    derham = Derham(TensorProductGrid(num_elements=(16, 1, 1)), DerhamOptions(degree=(3, 1, 1)), comm=comm)
+    mass_ops = WeightedMassOperators(derham, domain)
+    Propagator.derham = derham
+    Propagator.domain = domain
+    Propagator.mass_ops = mass_ops
+
+    def rho(e1, e2, e3):
+        return net + amp * xp.cos(2 * xp.pi * e1)
+
+    # uniform points without the duplicate endpoint, such that the mean of cos(k x) vanishes
+    e1 = xp.linspace(0.0, 1.0, 100, endpoint=False)
+    x = domain(e1, 0.0, 0.0)[0][:, 0, 0]
+    phi_exact = amp / k**2 * xp.cos(k * x)
+    e_exact = amp / k * xp.sin(k * x)
+
+    def solve(**kwargs):
+        phi = FEECVariable(space="H1")
+        phi.allocate(derham=derham, domain=domain)
+        poisson = PoissonSolve(rho=rho)
+        poisson.variables.phi = phi
+        poisson.options = poisson.Options(
+            stab_mat=stab_mat,
+            solver_params=SolverParameters(tol=1e-13, maxiter=3000, recycle=False),
+            **kwargs,
+        )
+        poisson.allocate()
+        poisson(1.0)
+        e = FEECVariable(space="Hcurl")
+        e.allocate(derham=derham, domain=domain)
+        derham.grad.dot(-phi.spline.vector, out=e.spline.vector)
+        phi_vals = domain.push(phi.spline, e1, 0.0, 0.0, kind="0")[:, 0, 0]
+        e_vals = domain.push(e.spline, e1, 0.0, 0.0, kind="1")[0][:, 0, 0]
+        return poisson, phi_vals, e_vals
+
+    # default options (stab_eps -> 1e-14): phi is O(1) and E is accurate
+    poisson, phi_vals, e_vals = solve()
+    assert poisson.operator_is_singular
+    e_err = xp.max(xp.abs(e_vals - e_exact)) / xp.max(xp.abs(e_exact))
+    phi_err = xp.max(xp.abs(phi_vals - xp.mean(phi_vals) - phi_exact)) / xp.max(xp.abs(phi_exact))
+    logger.info(f"{stab_mat=}: {xp.max(xp.abs(phi_vals))=}, {e_err=}, {phi_err=}")
+    assert xp.max(xp.abs(phi_vals)) < 2.0 * xp.max(xp.abs(phi_exact))
+    assert e_err < 1e-3
+    assert phi_err < 1e-3
+
+    # without compatibility, the constant mode of phi is huge: (net charge) / stab_eps
+    _, phi_vals_nc, _ = solve(enforce_compatibility=False)
+    logger.info(f"{stab_mat=}, enforce_compatibility=False: {xp.max(xp.abs(phi_vals_nc))=}")
+    assert xp.max(xp.abs(phi_vals_nc)) > 1e10
+
+    # for stab_mat="M0" and a finite stab_eps, the projection changes phi only by a constant
+    if stab_mat == "M0":
+        _, phi1, e1_vals = solve(stab_eps=1e-2)
+        _, phi2, e2_vals = solve(stab_eps=1e-2, enforce_compatibility=False)
+        assert xp.max(xp.abs(e1_vals - e2_vals)) < 1e-10 * xp.max(xp.abs(e_exact))
+        shift = phi2 - phi1
+        assert xp.max(xp.abs(shift - xp.mean(shift))) < 1e-10 * xp.max(xp.abs(phi_exact))
+
+
+def test_poisson_dirichlet_not_singular():
+    """With Dirichlet boundary conditions the Poisson operator is regular and the rhs is not modified."""
+    domain = domains.Cuboid(r1=4.0)
+    derham = Derham(
+        TensorProductGrid(num_elements=(8, 1, 1)),
+        DerhamOptions(degree=(2, 1, 1), bcs=(("dirichlet", "dirichlet"), None, None)),
+        comm=comm,
+    )
+    Propagator.derham = derham
+    Propagator.domain = domain
+    Propagator.mass_ops = WeightedMassOperators(derham, domain)
+
+    phi = FEECVariable(space="H1")
+    phi.allocate(derham=derham, domain=domain)
+    poisson = PoissonSolve(rho=lambda e1, e2, e3: 1.0 + 0.0 * e1)
+    poisson.variables.phi = phi
+    poisson.allocate()
+    assert not poisson.operator_is_singular
+    assert poisson._ones is None
+
+
 @pytest.mark.parametrize("stab_mat", ["M0", "Id"])
 @pytest.mark.parametrize("bc_type", ["periodic", "dirichlet", "neumann"])
 def test_poisson_2d_kronecker_preconds(bc_type, stab_mat):