diff --git a/pyproject.toml b/pyproject.toml
index 6334815bf..19566a1b4 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -176,4 +176,4 @@ markers = [
     "kinetic",
     "hybrid",
     "single",
-]
+]
\ No newline at end of file
diff --git a/src/struphy/kinetic_background/base.py b/src/struphy/kinetic_background/base.py
index 765ee1508..5bccda827 100644
--- a/src/struphy/kinetic_background/base.py
+++ b/src/struphy/kinetic_background/base.py
@@ -580,3 +580,281 @@ def add_perturbation(self) -> bool:
     def add_perturbation(self, new):
         assert isinstance(new, bool)
         self._add_perturbation = new
+
+
+class CanonicalMaxwellian(KineticBackground):
+    r"""canonical Maxwellian distribution function.
+    It is defined by three constants of motion in the axissymmetric toroidal system:
+
+    - Shifted canonical toroidal momentum
+    .. math::
+
+        \psi_c = \psi + \frac{m_s F}{q_s B}v_\parallel - \text{sign}(v_\parallel)\sqrt{2(\epsilon - \mu B)}\frac{m_sF}{q_sB} \mathcal{H}(\epsilon - \mu B),
+
+    - Energy
+    .. math::
+
+        \epsilon = \frac{1}{2}m_sv_\parallel² + \mu B,
+
+    - Magnetic moment
+    .. math::
+
+        \mu = \frac{m_s v_\perp²}{2B},
+
+    where :math:`\psi` is the poloidal magnetic flux function, :math:`F=F(\psi)` is the poloidal current function and :math:`\mathcal{H}` is the Heaviside function.
+    With the three constants of motion, a canonical Maxwellian distribution function is defined as
+
+    .. math::
+
+        F(\psi_c, \epsilon, \mu) = \frac{n(\psi_c)}{(2\pi)^{3/2}v_\text{th}³(\psi_c)} \text{exp}\left[ - \frac{\epsilon}{v_\text{th}²(\psi_c)}\right].
+
+    """
+
+    @abstractmethod
+    def vth(self, *etas):
+        """Thermal velocities (0-forms).
+
+        Parameters
+        ----------
+        etas : numpy.arrays
+            Evaluation points. All arrays must be of same shape (can be 1d for flat evaluation).
+
+        Returns
+        -------
+        A list[float] (background values) or a list[numpy.array] of the evaluated thermal velocities.
+        """
+        pass
+
+    @abstractmethod
+    def eval_energy(self, *etas):
+        """Energy evaluated at given particle positions and velocities."""
+        pass
+
+    @abstractmethod
+    def eval_psic(self, *etas):
+        """Shifted canonical toroidal momentum evaluated at given particle positions and velocities."""
+        pass
+
+    @abstractmethod
+    def eval_rc(self, *etas):
+        """Callable function return evaluation points ( :math:`r_c` ) from particle positions ( :math:`\eta_1, \eta_2, \eta_3` ).
+
+        Parameters
+        ----------
+        etas : numpy.arrays
+            Evaluation points. All arrays must be of same shape (can be 1d for flat evaluation).
+
+        Returns
+        -------
+        A list[float] (background values) or a list[numpy.array] of the evaluation points (etas).
+        """
+        pass
+
+    @property
+    @abstractmethod
+    def maxw_params(self) -> dict:
+        """Parameters dictionary defining moments of the Maxwellian."""
+
+    def check_maxw_params(self):
+        for k, v in self.maxw_params.items():
+            assert isinstance(k, str)
+            assert isinstance(v, tuple), f"Maxwallian parameter {k} must be tuple, but is {v}"
+            assert len(v) == 2
+
+            assert isinstance(v[0], (float, int, Callable))
+            assert isinstance(v[1], Perturbation) or v[1] is None
+
+    @classmethod
+    def gaussian(self, e, vth=1.0):
+        """3-dim. normal distribution, to which array-valued thermal velocities can be passed.
+
+        Parameters
+        ----------
+        v : float | array-like
+            Velocity coordinate(s).
+
+        vth : float | array-like
+            Thermal velocity evaluated at position array.
+
+        Returns
+        -------
+        An array of size(e).
+        """
+
+        if isinstance(vth, xp.ndarray):
+            assert e.shape == vth.shape, f"{e.shape = } but {vth.shape = }"
+
+        out = 2.0 * xp.sqrt(e / xp.pi) / vth**3 * xp.exp(-e / vth**2)
+
+        return out
+
+    def __call__(self, *args):
+        """Evaluates the canonical Maxwellian distribution function.
+
+        There are two use-cases for this function in the code:
+
+        1. Evaluating for particles ("flat evaluation", inputs are all 1D of length N_p)
+        2. Evaluating the function on a meshgrid (constants of motion).
+
+        Hence all arguments must always have
+
+        1. the same shape
+        2. either ndim = 1 or ndim = 3 (energy, mu, canonical toroidal momentum).
+
+        Parameters
+        ----------
+        *args : array_like
+            Constants of motion arguments in the order eta1, eta2, eta3, v1, ..., vn.
+
+        Returns
+        -------
+        f : xp.ndarray
+            The evaluated Maxwellian.
+        """
+
+        # Check that all args have the same shape
+        shape0 = xp.shape(args[0])
+        for i, arg in enumerate(args):
+            assert xp.shape(arg) == shape0, f"Argument {i} has {xp.shape(arg) = }, but must be {shape0 = }."
+            assert xp.ndim(arg) == 1 or xp.ndim(arg) == 3, (
+                f"{xp.ndim(arg) = } not allowed for canonical Maxwellian evaluation."
+            )  # flat or meshgrid evaluation
+
+        # Get result evaluated at eta's
+        res = self.n(*args)
+        vths = self.vth(*args)
+
+        # take care of correct broadcasting, assuming args come from constants of motion meshgrid
+        if xp.ndim(args[0]) == 3:
+            # move eta axes to the back
+            arg_t = xp.moveaxis(args[0], 0, -1)
+            arg_t = xp.moveaxis(arg_t, 0, -1)
+            arg_t = xp.moveaxis(arg_t, 0, -1)
+
+            # broadcast
+            res_broad = res + 0.0 * arg_t
+
+            # move eta axes to the front
+            res = xp.moveaxis(res_broad, -1, 0)
+            res = xp.moveaxis(res, -1, 0)
+            res = xp.moveaxis(res, -1, 0)
+
+        # Multiply result with gaussian in energy
+        # correct broadcasting
+        if xp.ndim(args[0]) == 3:
+            vth_broad = vths + 0.0 * arg_t
+            vth = xp.moveaxis(vth_broad, -1, 0)
+            vth = xp.moveaxis(vth, -1, 0)
+            vth = xp.moveaxis(vth, -1, 0)
+        else:
+            vth = vths
+
+        e = self.eval_energy(*args)
+        res *= self.gaussian(e, vth=vth)
+
+        return res
+
+    def _evaluate_moment(self, eta1, eta2, eta3, vparallel, mu, *, name: str = "n", add_perturbation: bool = None):
+        """Scalar moment evaluation as background + perturbation.
+
+        Parameters
+        ----------
+        eta1, eta2, eta3 : numpy.arrays
+            Evaluation points. All arrays must be of same shape (can be 1d for flat evaluation).
+
+        vparallel : numpy.array
+            Parallel velocity.
+
+        mu : numpy.array
+            Magnetic moment.
+
+        name : str
+            Which moment to evaluate (see varaible "dct" below).
+
+        add_perturbation : bool | None
+            Whether to add the perturbation defined in maxw_params. If None, is taken from self.add_perturbation.
+
+        Returns
+        -------
+        A float (background value) or a numpy.array of the evaluated scalar moment.
+        """
+
+        # collect arguments
+        assert isinstance(eta1, xp.ndarray)
+        assert isinstance(eta2, xp.ndarray)
+        assert isinstance(eta3, xp.ndarray)
+        assert isinstance(vparallel, xp.ndarray)
+        assert isinstance(mu, xp.ndarray)
+
+        params = self.maxw_params[name]
+        assert isinstance(params, tuple)
+        assert len(params) == 2
+
+        # flat evaluation for markers
+        if eta1.ndim == 1:
+            etas = [
+                xp.concatenate(
+                    (eta1[:, None], eta2[:, None], eta3[:, None]),
+                    axis=1,
+                ),
+            ]
+        # assuming that input comes from meshgrid.
+        elif eta1.ndim == 4:
+            etas = (
+                eta1[:, :, :, 0],
+                eta2[:, :, :, 0],
+                eta3[:, :, :, 0],
+            )
+        elif eta1.ndim == 5:
+            etas = (
+                eta1[:, :, :, 0, 0],
+                eta2[:, :, :, 0, 0],
+                eta3[:, :, :, 0, 0],
+            )
+        elif eta1.ndim == 6:
+            etas = (
+                eta1[:, :, :, 0, 0, 0],
+                eta2[:, :, :, 0, 0, 0],
+                eta3[:, :, :, 0, 0, 0],
+            )
+        else:
+            etas = (eta1, eta2, eta3)
+
+        # initialize output
+        if eta1.ndim == 1:
+            out = 0.0 * eta1
+        else:
+            out = 0.0 * etas[0]
+
+        # evaluate background
+        background = params[0]
+        if isinstance(background, (float, int)):
+            out += background
+        else:
+            assert callable(background)
+            # if eta1.ndim == 1:
+            #     out += background(eta1, eta2, eta3)
+            # else:
+            out += background(self.eval_rc(eta1, eta2, eta3, vparallel, mu))
+
+        # add perturbation
+        if add_perturbation is None:
+            add_perturbation = self.add_perturbation
+
+        perturbation = params[1]
+        if perturbation is not None and add_perturbation:
+            assert isinstance(perturbation, Perturbation)
+            out += perturbation(self.eval_rc(eta1, eta2, eta3, vparallel, mu))
+
+        return out
+
+    @property
+    def add_perturbation(self) -> bool:
+        if not hasattr(self, "_add_perturbation"):
+            self._add_perturbation = True
+        return self._add_perturbation
+
+    @add_perturbation.setter
+    def add_perturbation(self, new):
+        assert isinstance(new, bool)
+        self._add_perturbation = new
diff --git a/src/struphy/kinetic_background/maxwellians.py b/src/struphy/kinetic_background/maxwellians.py
index d9e34dcab..ab53b2ac6 100644
--- a/src/struphy/kinetic_background/maxwellians.py
+++ b/src/struphy/kinetic_background/maxwellians.py
@@ -7,7 +7,7 @@
 from struphy.fields_background.base import FluidEquilibriumWithB
 from struphy.fields_background.equils import set_defaults
 from struphy.initial.base import Perturbation
-from struphy.kinetic_background.base import Maxwellian
+from struphy.kinetic_background.base import CanonicalMaxwellian, Maxwellian
 
 
 class Maxwellian3D(Maxwellian):
@@ -245,14 +245,23 @@ def velocity_jacobian_det(self, eta1, eta2, eta3, *v):
         -------
         """
 
-        assert eta1.ndim == 1
-        assert eta2.ndim == 1
-        assert eta3.ndim == 1
-        assert len(v) == 2
-
-        # call equilibrium
-        etas = (xp.vstack((eta1, eta2, eta3)).T).copy()
-        absB0 = self.equil.absB0(etas)
+        # collect arguments
+        assert isinstance(eta1, xp.ndarray)
+        assert isinstance(eta2, xp.ndarray)
+        assert isinstance(eta3, xp.ndarray)
+        assert isinstance(v[0], xp.ndarray)
+        assert isinstance(v[1], xp.ndarray)
+        assert eta1.shape == eta2.shape == eta3.shape == v[0].shape == v[1].shape
+        assert eta1.ndim == 1, "Input arguments must be a marker array."
+
+        etas = [
+            xp.concatenate(
+                (eta1[:, None], eta2[:, None], eta3[:, None]),
+                axis=1,
+            ),
+        ]
+
+        absB0 = self.equil.absB0(*etas)
 
         # J = v_perp/B
         jacobian_det = v[1] / absB0
@@ -301,38 +310,16 @@ def vth(self, eta1, eta2, eta3):
         return [ou * mom_fac for ou, mom_fac in zip(out, self.moment_factors["vth"])]
 
 
-class CanonicalMaxwellian:
-    r"""canonical Maxwellian distribution function.
-    It is defined by three constants of motion in the axissymmetric toroidal system:
-
-    - Shifted canonical toroidal momentum
-
-    .. math::
-
-        \psi_c = \psi + \frac{m_s F}{q_s B}v_\parallel - \text{sign}(v_\parallel)\sqrt{2(\epsilon - \mu B)}\frac{m_sF}{q_sB} \mathcal{H}(\epsilon - \mu B),
-
-    - Energy
-
-    .. math::
-
-        \epsilon = \frac{1}{2}m_sv_\parallel² + \mu B,
-
-    - Magnetic moment
-
-    .. math::
-
-        \mu = \frac{m_s v_\perp²}{2B},
-
-    where :math:`\psi` is the poloidal magnetic flux function, :math:`F=F(\psi)` is the poloidal current function and :math:`\mathcal{H}` is the Heaviside function.
-
-    With the three constants of motion, a canonical Maxwellian distribution function is defined as
-
-    .. math::
-
-        F(\psi_c, \epsilon, \mu) = \frac{n(\psi_c)}{(2\pi)^{3/2}v_\text{th}³(\psi_c)} \text{exp}\left[ - \frac{\epsilon}{v_\text{th}²(\psi_c)}\right].
+class CanonicalMaxwellian2D(CanonicalMaxwellian):
+    r"""Guiding-center :class:`~struphy.kinetic_background.base.CanonicalMaxwellian` depending on
+    three constants of motion :math:`\epsilon, \mu, \psi_c`, in a axissymmetric toroidal system:
 
     Parameters
     ----------
+    n, vth : tuple
+        Moments of the canonical Maxwellian as tuples. The first entry defines the background
+        (float for constant background or callable), the second entry defines a Perturbation (can be None).
+
     maxw_params : dict
         Parameters for the kinetic background.
 
@@ -354,10 +341,12 @@ def __init__(
         vth: tuple[float | Callable, Perturbation] = (1.0, None),
         equil: FluidEquilibriumWithB = None,
         volume_form: bool = True,
+        epsilon: float = 1.,
     ):
         self._maxw_params = {}
         self._maxw_params["n"] = n
         self._maxw_params["vth"] = vth
+        self._epsilon = epsilon
 
         self.check_maxw_params()
 
@@ -376,6 +365,16 @@ def coords(self):
         r"""Coordinates of the CanonicalMaxwellian, :math:`(\epsilon, \mu, \psi_c)`."""
         return "constants_of_motion"
 
+    @property
+    def vdim(self):
+        """Dimension of the velocity space."""
+        return 2
+
+    @property
+    def is_polar(self):
+        """List of booleans of length vdim. True for a velocity coordinate that is a radial polar coordinate (v_perp)."""
+        return [False, True]
+
     @property
     def maxw_params(self):
         """Parameters dictionary defining constant moments of the Maxwellian."""
@@ -397,116 +396,38 @@ def check_maxw_params(self):
             assert isinstance(v[0], (float, int, Callable))
             assert isinstance(v[1], Perturbation) or v[1] is None
 
-    def velocity_jacobian_det(self, eta1, eta2, eta3, energy):
+    def velocity_jacobian_det(self, eta1, eta2, eta3, vparallel, mu):
         r"""TODO"""
+        # collect arguments
+        assert isinstance(eta1, xp.ndarray)
+        assert isinstance(eta2, xp.ndarray)
+        assert isinstance(eta3, xp.ndarray)
+        assert isinstance(vparallel, xp.ndarray)
+        assert isinstance(mu, xp.ndarray)
+        assert eta1.shape == eta2.shape == eta3.shape == vparallel.shape == mu.shape
+        assert eta1.ndim == 1, "Input arguments must be a marker array."
 
-        assert eta1.ndim == 1
-        assert eta2.ndim == 1
-        assert eta3.ndim == 1
-
-        if self.maxw_params["type"] == "Particles6D":
-            return xp.sqrt(2.0 * energy) * 4.0 * xp.pi
-
-        else:
-            # call equilibrium
-            etas = (xp.vstack((eta1, eta2, eta3)).T).copy()
-            absB0 = self.equil.absB0(etas)
-
-            return xp.sqrt(energy) * 2.0 * xp.sqrt(2.0) / absB0
-
-    def gaussian(self, e, vth=1.0):
-        """3-dim. normal distribution, to which array-valued thermal velocities can be passed.
-
-        Parameters
-        ----------
-        e : float | array-like
-            Energy.
-
-        vth : float | array-like
-            Thermal velocity evaluated at psic.
-
-        Returns
-        -------
-        An array of size(e).
-        """
-
-        if isinstance(vth, xp.ndarray):
-            assert e.shape == vth.shape, f"{e.shape =} but {vth.shape =}"
-
-        return 2.0 * xp.sqrt(e / xp.pi) / vth**3 * xp.exp(-e / vth**2)
-
-    def __call__(self, *args):
-        """Evaluates the canonical Maxwellian distribution function.
-
-        There are two use-cases for this function in the code:
-
-        1. Evaluating for particles ("flat evaluation", inputs are all 1D of length N_p)
-        2. Evaluating the function on a meshgrid (in phase space).
-
-        Hence all arguments must always have
-
-        1. the same shape
-        2. either ndim = 1 or ndim = 3.
-
-        Parameters
-        ----------
-        *args : array_like
-            Position-velocity arguments in the order energy, magnetic moment, canonical toroidal momentum.
-
-        Returns
-        -------
-        f : xp.ndarray
-            The evaluated Maxwellian.
-        """
+        etas = [
+            xp.concatenate(
+                (eta1[:, None], eta2[:, None], eta3[:, None]),
+                axis=1,
+            ),
+        ]
 
-        # Check that all args have the same shape
-        shape0 = xp.shape(args[0])
-        for i, arg in enumerate(args):
-            assert xp.shape(arg) == shape0, f"Argument {i} has {xp.shape(arg) =}, but must be {shape0 =}."
-            assert xp.ndim(arg) == 1 or xp.ndim(arg) == 3, (
-                f"{xp.ndim(arg) =} not allowed for canonical Maxwellian evaluation."
-            )  # flat or meshgrid evaluation
-
-        # Get result evaluated with each particles' psic
-        res = self.n(args[2])
-        vths = self.vth(args[2])
-
-        # take care of correct broadcasting, assuming args come from phase space meshgrid
-        if xp.ndim(args[0]) == 3:
-            # move eta axes to the back
-            arg_t = xp.moveaxis(args[0], 0, -1)
-            arg_t = xp.moveaxis(arg_t, 0, -1)
-            arg_t = xp.moveaxis(arg_t, 0, -1)
-
-            # broadcast
-            res_broad = res + 0.0 * arg_t
-
-            # move eta axes to the front
-            res = xp.moveaxis(res_broad, -1, 0)
-            res = xp.moveaxis(res, -1, 0)
-            res = xp.moveaxis(res, -1, 0)
-
-        # Multiply result with gaussian in energy
-        if xp.ndim(args[0]) == 3:
-            vth_broad = vths + 0.0 * arg_t
-            vth = xp.moveaxis(vth_broad, -1, 0)
-            vth = xp.moveaxis(vth, -1, 0)
-            vth = xp.moveaxis(vth, -1, 0)
-        else:
-            vth = vths
+        absB0 = self.equil.absB0(*etas)
 
-        res *= self.gaussian(args[0], vth=vth)
+        energy = self.eval_energy(eta1, eta2, eta3, vparallel, mu)
 
-        return res
+        return xp.sqrt(energy) * 2.0 * xp.sqrt(2.0) / absB0
 
     @property
-    def volume_form(self):
+    def volume_form(self) -> bool:
         """Boolean. True if the background is represented as a volume form (thus including the velocity Jacobian |v_perp|)."""
         return self._volume_form
 
     @property
     def moment_factors(self):
-        """Collection of factors multiplied onto the defined moments n, u, and vth."""
+        """Collection of factors multiplied onto the defined moments n and vth."""
         return self._moment_factors
 
     @moment_factors.setter
@@ -514,7 +435,62 @@ def moment_factors(self, **kwargs):
         for kw, arg in kwargs:
             self._moment_factors[kw] = arg
 
-    def rc(self, psic):
+    def eval_energy(self, eta1, eta2, eta3, vparallel, mu):
+        r"""Energy evaluated at given particle positions and velocities."""
+        # call domain and equilibrium information
+        if eta1.ndim == 1:
+            etas = [
+                xp.concatenate(
+                    (eta1[:, None], eta2[:, None], eta3[:, None]),
+                    axis=1,
+                ),
+            ]
+            absB0 = self.equil.absB0(*etas)
+        else:
+            absB0 = self.equil.absB0(eta1, eta2, eta3)
+
+        # calculate energy
+        energy = 1/2 * vparallel**2 + mu * absB0
+
+        return energy
+
+    def eval_psic(self, eta1, eta2, eta3, vparallel, mu):
+        r"""Shifted canonical toroidal momentum evaluated at given particle positions and velocities."""
+        # call domain and equilibrium information
+        a1 = self.equil.domain.params["a1"]
+        B0 = self.equil.params["B0"]
+        R0 = self.equil.params["R0"]
+        if eta1.ndim == 1:
+            etas = [
+                xp.concatenate(
+                    (eta1[:, None], eta2[:, None], eta3[:, None]),
+                    axis=1,
+                ),
+            ]
+            absB0 = self.equil.absB0(*etas)
+        else:
+            absB0 = self.equil.absB0(eta1, eta2, eta3)
+        psi = self.equil.psi_r(eta1 * (1-a1) + a1)
+
+        # calculate energy
+        energy = self.eval_energy(eta1, eta2, eta3, vparallel, mu)
+
+        # calculate psic
+        psic = psi - self._epsilon * B0 * R0 / absB0 * vparallel
+
+        positive_mask = (energy - mu * B0) > 0
+        correction = xp.zeros_like(psic)
+        correction[positive_mask] = (
+            self._epsilon
+            * xp.sign(vparallel[positive_mask])
+            * xp.sqrt(2 * (energy[positive_mask] - mu[positive_mask] * B0))
+            * R0
+        )
+        psic += correction
+
+        return psic
+
+    def eval_rc(self, eta1, eta2, eta3, vparallel, mu):
         r""" Square root of radially normalized canonical toroidal momentum.
 
         .. math::
@@ -528,17 +504,9 @@ def rc(self, psic):
             \end{aligned}
 
         where :math:`\psi_\text{axis}` and :math:`\psi_\text{edge}` are poloidal magnetic flux function at the center and edge of poloidal plane respectively.
-
-        Parameters
-        ----------
-        psic : numpy.arrays
-            Evaluation points. All arrays must be of same shape (can be 1d for flat evaluation).
-
-        Returns
-        -------
-        A numpy.array of the evaluated :math:`r_c`.
-
         """
+        # calculate psic
+        psic = self.eval_psic(eta1, eta2, eta3, vparallel, mu)
 
         # calculate rc²
         rc_squared = (psic - self.equil.psi_range[0]) / (self.equil.psi_range[1] - self.equil.psi_range[0])
@@ -555,74 +523,30 @@ def rc(self, psic):
 
         return rc
 
-    def n(self, psic, add_perturbation: bool = None):
-        """Density as background + perturbation.
-
-        Parameters
-        ----------
-        psic : numpy.array
-            Evaluation points. All arrays must be of same shape (can be 1d for flat evaluation).
-
-        Returns
-        -------
-        A float (background value) or a numpy.array of the evaluated density.
-        """
-        # collect arguments
-        assert isinstance(psic, xp.ndarray)
-
-        # assuming that input comes from meshgrid.
-        if psic.ndim == 3:
-            psic = psic[0, 0, :]
-
-        # set background density
-        if isinstance(self.maxw_params["n"][0], (float, int)):
-            res = self.maxw_params["n"][0] + 0.0 * psic
-        else:
-            nfun = self.maxw_params["n"][1]
-            # for typ, params in mom_funcs.items():
-            #     nfun = getattr(moment_functions, typ)(**params)
-            res = nfun(eta1=self.rc(psic))
-
-        # add perturbation
-        if add_perturbation is None:
-            add_perturbation = self.add_perturbation
-
-        perturbation = self.maxw_params["n"][1]
-        if perturbation is not None and add_perturbation:
-            assert isinstance(perturbation, Perturbation)
-            res = perturbation(eta1=self.rc(psic))
-            # if eta1.ndim == 1:
-            #     out += perturbation(eta1, eta2, eta3)
-            # else:
-            #     out += perturbation(*etas)
-
-        return res * self.moment_factors["n"]
-
-    def vth(self, psic):
-        """Thermal velocities as background + perturbation.
-
-        Parameters
-        ----------
-        psic : numpy.arrays
-            Evaluation points. All arrays must be of same shape (can be 1d for flat evaluation).
-
-        Returns
-        -------
-        A list[float] (background value) or a list[numpy.array] of the evaluated thermal velocities.
-        """
-
-        # collect arguments
-        assert isinstance(psic, xp.ndarray)
+    def n(self, eta1, eta2, eta3, vparallel, mu):
+        """Zero-th moment (density)."""
+        out = self._evaluate_moment(eta1, eta2, eta3, vparallel, mu, name="n")
+        return out * self.moment_factors["n"]
 
-        # assuming that input comes from meshgrid.
-        if psic.ndim == 3:
-            psic = psic[0, 0, :]
+    def u(self, eta1, eta2, eta3, vparallel, mu):
+        """Mean velocities."""
+        pass
 
-        res = self.maxw_params["vth"][0] + 0.0 * psic
+    def vth(self, eta1, eta2, eta3, vparallel, mu):
+        """Thermal velocities."""
+        out = self._evaluate_moment(eta1, eta2, eta3, vparallel, mu, name="vth")
+        return out * self.moment_factors["vth"]
 
-        # TODO: add perturbation
+    @property
+    def add_perturbation(self) -> bool:
+        if not hasattr(self, "_add_perturbation"):
+            self._add_perturbation = True
+        return self._add_perturbation
 
-        return res * self.moment_factors["vth"]
+    @add_perturbation.setter
+    def add_perturbation(self, new):
+        assert isinstance(new, bool)
+        self._add_perturbation = new
 
     @property
     def add_perturbation(self) -> bool:
diff --git a/src/struphy/kinetic_background/tests/test_maxwellians.py b/src/struphy/kinetic_background/tests/test_maxwellians.py
index d63b7d30a..85506b117 100644
--- a/src/struphy/kinetic_background/tests/test_maxwellians.py
+++ b/src/struphy/kinetic_background/tests/test_maxwellians.py
@@ -1503,7 +1503,7 @@ def test_canonical_maxwellian_uniform(Nel, show_plot=False):
     from struphy.fields_background import equils
     from struphy.geometry import domains
     from struphy.initial import perturbations
-    from struphy.kinetic_background.maxwellians import CanonicalMaxwellian
+    from struphy.kinetic_background.maxwellians import CanonicalMaxwellian2D
 
     e1 = xp.linspace(0.0, 1.0, Nel[0])
     e2 = xp.linspace(0.0, 1.0, Nel[1])
@@ -1516,7 +1516,7 @@ def test_canonical_maxwellian_uniform(Nel, show_plot=False):
 
     epsilon = 1.0
 
-    # evaluate three constants of motions at AdhocTorus equilibrium
+    # evaluate three constants of motion at AdhocTorus equilibrium
     AdhocTorus_params = {
         "a": 1.0,
         "R0": 10.0,
@@ -1563,7 +1563,7 @@ def test_canonical_maxwellian_uniform(Nel, show_plot=False):
     # ===========================================================
     maxw_params = {"n": 2.0, "vth": 1.0}
 
-    maxwellian = CanonicalMaxwellian(n=(2.0, None), vth=(1.0, None))
+    maxwellian = CanonicalMaxwellian2D(n=(2.0, None), vth=(1.0, None))
 
     # Test constant value at v_para = v_perp = 0.01
     res = maxwellian(energy, mu, psic).squeeze()
@@ -1673,7 +1673,7 @@ def test_canonical_maxwellian_uniform(Nel, show_plot=False):
     }
     pert = perturbations.ITPA_density(n0=n0, c=c)
 
-    maxwellian = CanonicalMaxwellian(n=(0.0, pert), equil=mhd_equil, volume_form=False)
+    maxwellian = CanonicalMaxwellian2D(n=(0.0, pert), equil=mhd_equil, volume_form=False)
 
     e1 = xp.linspace(0.0, 1.0, Nel[0])
     e2 = xp.linspace(0.0, 1.0, Nel[1])
@@ -1704,7 +1704,7 @@ def test_canonical_maxwellian_uniform(Nel, show_plot=False):
     res = maxwellian(energy, mu, psic).squeeze()
 
     # calculate rc
-    rc = maxwellian.rc(psic)
+    rc = maxwellian.psic_to_rc(psic)
 
     ana_res = n0 * c[3] * xp.exp(-c[2] / c[1] * xp.tanh((rc - c[0]) / c[2]))
     ana_res *= 2 * xp.sqrt(energy / xp.pi) / maxw_params["vth"] ** 3 * xp.exp(-energy / maxw_params["vth"] ** 2)
diff --git a/src/struphy/models/base.py b/src/struphy/models/base.py
index 0fcc46401..757c9d337 100644
--- a/src/struphy/models/base.py
+++ b/src/struphy/models/base.py
@@ -758,14 +758,14 @@ def update_distr_functions(self):
                 obj = var.particles
                 assert isinstance(obj, Particles)
 
-                if obj.n_cols_diagnostics > 0:
-                    for i in range(obj.n_cols_diagnostics):
+                if obj.n_cols_diag > 0:
+                    for i in range(obj.n_cols_diag):
                         str_dn = f"d{i + 1}"
                         dim_to_int[str_dn] = 3 + obj.vdim + 3 + i
 
                 for bin_plot in species.binning_plots:
                     comps = bin_plot.slice.split("_")
-                    components = [False] * (3 + obj.vdim + 3 + obj.n_cols_diagnostics)
+                    components = [False] * (3 + obj.vdim + 3 + obj.n_cols_diag)
 
                     for comp in comps:
                         components[dim_to_int[comp]] = True
diff --git a/src/struphy/models/species.py b/src/struphy/models/species.py
index 5d6936107..6715f2505 100644
--- a/src/struphy/models/species.py
+++ b/src/struphy/models/species.py
@@ -151,6 +151,8 @@ def set_markers(
         weights_params: WeightsParameters = None,
         boundary_params: BoundaryParameters = None,
         bufsize: float = 1.0,
+        n_cols_diag: int = 0,
+        n_cols_aux: int = 5,
     ):
         """Set marker parameters for loading, weight calculation, kernel density reconstruction
         and boundary conditions.
@@ -164,7 +166,15 @@ def set_markers(
         boundary_params : BoundaryParameters
 
         bufsize : float
-            Size of buffer (as multiple of total size, default=.25) in markers array."""
+            Size of buffer (as multiple of total size, default=.25) in markers array.
+
+        n_cols_diag : int
+            Number of diagnostics columns.
+
+        n_cols_aux : int
+            Number of auxiliary columns.
+        ----------
+        """
 
         # defaults
         if loading_params is None:
@@ -180,6 +190,8 @@ def set_markers(
         self.weights_params = weights_params
         self.boundary_params = boundary_params
         self.bufsize = bufsize
+        self.n_cols_diag = n_cols_diag
+        self.n_cols_aux = n_cols_aux
 
     def set_sorting_boxes(
         self,
diff --git a/src/struphy/models/variables.py b/src/struphy/models/variables.py
index f1a4d9f47..128b41362 100644
--- a/src/struphy/models/variables.py
+++ b/src/struphy/models/variables.py
@@ -232,6 +232,8 @@ def allocate(
             weights_params=self.species.weights_params,
             boundary_params=self.species.boundary_params,
             bufsize=self.species.bufsize,
+            n_cols_diag=self.species.n_cols_diag,
+            n_cols_aux=self.species.n_cols_aux,
             domain=domain,
             equil=equil,
             projected_equil=projected_equil,
@@ -383,7 +385,6 @@ def allocate(
             background=self.backgrounds,
             n_as_volume_form=self.n_as_volume_form,
             perturbations=self.perturbations,
-            equation_params=self.species.equation_params,
             verbose=verbose,
         )
 
diff --git a/src/struphy/pic/accumulation/accum_kernels_gc.py b/src/struphy/pic/accumulation/accum_kernels_gc.py
index fecf6a255..e7e13b77c 100644
--- a/src/struphy/pic/accumulation/accum_kernels_gc.py
+++ b/src/struphy/pic/accumulation/accum_kernels_gc.py
@@ -640,6 +640,7 @@ def cc_lin_mhd_5d_M(
     "curl_norm_b",
     "norm_b1",
     "grad_PB",
+    "grad_PBeq",
 )
 def cc_lin_mhd_5d_gradB(
     args_markers: "MarkerArguments",
@@ -668,6 +669,9 @@ def cc_lin_mhd_5d_gradB(
     grad_PB1: "float[:,:,:]",
     grad_PB2: "float[:,:,:]",
     grad_PB3: "float[:,:,:]",
+    grad_PBeq1: "float[:,:,:]",
+    grad_PBeq2: "float[:,:,:]",
+    grad_PBeq3: "float[:,:,:]",
     basis_u: "int",
 ):
     r"""Accumulation kernel for the propagator :class:`~struphy.propagators.propagators_coupling.CurrentCoupling5DGradB`.
@@ -700,7 +704,9 @@ def cc_lin_mhd_5d_gradB(
     """
 
     markers = args_markers.markers
+    n_markers = args_markers.n_markers
     Np = args_markers.Np
+    first_init_idx = args_markers.first_init_idx
 
     # allocate for magnetic field evaluation
     b = empty(3, dtype=float)
@@ -710,6 +716,7 @@ def cc_lin_mhd_5d_gradB(
     curl_norm_b = empty(3, dtype=float)
     norm_b1 = empty(3, dtype=float)
     grad_PB = empty(3, dtype=float)
+    grad_PBeq = empty(3, dtype=float)
 
     # allocate for metric coeffs
     dfm = empty((3, 3), dtype=float)
@@ -723,14 +730,15 @@ def cc_lin_mhd_5d_gradB(
 
     tmp_v = empty(3, dtype=float)
 
-    # get number of markers
-    n_markers_loc = shape(markers)[0]
-
-    for ip in range(n_markers_loc):
+    for ip in range(n_markers):
         # only do something if particle is a "true" particle (i.e. not a hole)
         if markers[ip, 0] == -1.0:
             continue
 
+        # if particle is refilled
+        if markers[ip, first_init_idx] == -1.0:
+            continue
+
         # marker positions
         eta1 = markers[ip, 0]
         eta2 = markers[ip, 1]
@@ -766,6 +774,9 @@ def cc_lin_mhd_5d_gradB(
         # grad_PB; 1form
         eval_1form_spline_mpi(span1, span2, span3, args_derham, grad_PB1, grad_PB2, grad_PB3, grad_PB)
 
+        # grad_PBeq; 1form
+        eval_1form_spline_mpi(span1, span2, span3, args_derham, grad_PBeq1, grad_PBeq2, grad_PBeq3, grad_PBeq)
+
         # b_star; 2form transformed into H1vec
         b_star[:] = b + curl_norm_b * v * epsilon
 
@@ -795,19 +806,12 @@ def cc_lin_mhd_5d_gradB(
 
             # call the appropriate matvec filler
             particle_to_mat_kernels.vec_fill_v0vec(
-                args_derham,
-                span1,
-                span2,
-                span3,
-                vec1,
-                vec2,
-                vec3,
-                filling_v[0],
-                filling_v[1],
-                filling_v[2],
+                args_derham, span1, span2, span3, vec1, vec2, vec3, filling_v[0], filling_v[1], filling_v[2]
             )
 
         elif basis_u == 2:
+            grad_PB += grad_PBeq
+
             linalg_kernels.matrix_matrix(b_prod, norm_b_prod, tmp)
             linalg_kernels.matrix_vector(tmp, grad_PB, tmp_v)
 
@@ -815,16 +819,7 @@ def cc_lin_mhd_5d_gradB(
 
             # call the appropriate matvec filler
             particle_to_mat_kernels.vec_fill_v2(
-                args_derham,
-                span1,
-                span2,
-                span3,
-                vec1,
-                vec2,
-                vec3,
-                filling_v[0],
-                filling_v[1],
-                filling_v[2],
+                args_derham, span1, span2, span3, vec1, vec2, vec3, filling_v[0], filling_v[1], filling_v[2]
             )
     vec1 /= Np
     vec2 /= Np
@@ -1252,16 +1247,7 @@ def cc_lin_mhd_5d_gradB_dg(
 
             # call the appropriate matvec filler
             particle_to_mat_kernels.vec_fill_v0vec(
-                args_derham,
-                span1,
-                span2,
-                span3,
-                vec1,
-                vec2,
-                vec3,
-                filling_v[0],
-                filling_v[1],
-                filling_v[2],
+                args_derham, span1, span2, span3, vec1, vec2, vec3, filling_v[0], filling_v[1], filling_v[2]
             )
 
         elif basis_u == 2:
@@ -1299,16 +1285,7 @@ def cc_lin_mhd_5d_gradB_dg(
 
             # call the appropriate matvec filler
             particle_to_mat_kernels.vec_fill_v2(
-                args_derham,
-                span1,
-                span2,
-                span3,
-                vec1,
-                vec2,
-                vec3,
-                filling_v[0],
-                filling_v[1],
-                filling_v[2],
+                args_derham, span1, span2, span3, vec1, vec2, vec3, filling_v[0], filling_v[1], filling_v[2]
             )
 
     vec1 /= Np
diff --git a/src/struphy/pic/accumulation/filter.py b/src/struphy/pic/accumulation/filter.py
index 42bf0885c..70aafb834 100644
--- a/src/struphy/pic/accumulation/filter.py
+++ b/src/struphy/pic/accumulation/filter.py
@@ -1,7 +1,6 @@
 from dataclasses import dataclass
 
 import cunumpy as xp
-import numpy as np
 from scipy.fft import irfft, rfft
 
 from struphy.feec.psydac_derham import Derham
@@ -165,7 +164,7 @@ def _apply_toroidal_fourier_filter(self, vec, modes: tuple[int, ...]):
         else:
             vec_temp = xp.zeros(int((tor_Nel - 1) / 2) + 1, dtype=complex)
 
-        for axis, comp, starts, ends in self._yield_dir_components(vec):
+        for _, comp, starts, ends in self._yield_dir_components(vec):
             for i in range(3):
                 ir[i] = int(ends[i] + 1 - starts[i])
 
@@ -183,4 +182,4 @@ def _apply_toroidal_fourier_filter(self, vec, modes: tuple[int, ...]):
                     # inverse FFT back to real space, write in-place
                     comp._data[ii, jj, pn[2] : pn[2] + ir[2]] = irfft(vec_temp, n=tor_Nel)
 
-        vec.update_ghost_regions()
+            comp.update_ghost_regions()
diff --git a/src/struphy/pic/accumulation/particles_to_grid.py b/src/struphy/pic/accumulation/particles_to_grid.py
index 099faced3..23268ef72 100644
--- a/src/struphy/pic/accumulation/particles_to_grid.py
+++ b/src/struphy/pic/accumulation/particles_to_grid.py
@@ -219,7 +219,7 @@ def __call__(self, *optional_args, **args_control):
                 vec.update_ghost_regions()
 
                 self.accfilter(vec)
-                vec_finished = True
+            vec_finished = True
 
         if self.particles.clone_config is None:
             num_clones = 1
diff --git a/src/struphy/pic/base.py b/src/struphy/pic/base.py
index 999188782..b0f961c78 100644
--- a/src/struphy/pic/base.py
+++ b/src/struphy/pic/base.py
@@ -167,6 +167,8 @@ def __init__(
         perturbations: dict[str, Perturbation] = None,
         n_as_volume_form: bool = False,
         equation_params: dict = None,
+        n_cols_diag: int = None,
+        n_cols_aux: int = None,
         verbose: bool = False,
     ):
         self._clone_config = clone_config
@@ -198,6 +200,8 @@ def __init__(
         self._equil = equil
         self._projected_equil = projected_equil
         self._equation_params = equation_params
+        self._n_cols_diag = n_cols_diag
+        self._n_cols_aux = n_cols_aux
 
         # check for mpi communicator (i.e. sub_comm of clone)
         if self.mpi_comm is None:
@@ -372,16 +376,19 @@ def vdim(self):
         pass
 
     @property
-    @abstractmethod
-    def n_cols_diagnostics(self):
+    def n_cols_diag(self):
         """Number of columns for storing diagnostics for each marker."""
-        pass
+        return self._n_cols_diag
 
     @property
-    @abstractmethod
     def n_cols_aux(self):
         """Number of auxiliary columns for each marker (e.g. for storing evaluation data)."""
-        pass
+        return self._n_cols_aux
+
+    @property
+    def equation_params(self):
+        """Parameters appearing in model equation due to Struphy normalization."""
+        return self._equation_params
 
     @property
     def first_diagnostics_idx(self):
@@ -392,7 +399,7 @@ def first_diagnostics_idx(self):
     @property
     def first_pusher_idx(self):
         """Starting index for storing initial conditions for a Pusher call."""
-        return self.first_diagnostics_idx + self.n_cols_diagnostics
+        return self.first_diagnostics_idx + self.n_cols_diag
 
     @property
     def n_cols_pusher(self):
@@ -567,11 +574,6 @@ def verbose(self):
         """Show some more particles info."""
         return self._verbose
 
-    @property
-    def equation_params(self):
-        """Parameters appearing in model equation due to Struphy normalization."""
-        return self._equation_params
-
     @property
     def initial_condition(self) -> KineticBackground:
         """Kinetic initial condition"""
@@ -1021,36 +1023,13 @@ def _set_background_function(self):
         #                 self._f0 = self._f0 + bckgr
 
     def _set_background_coordinates(self):
-        if self.type != "sph" and self.f0.coords == "constants_of_motion":
-            # Particles6D
-            if self.vdim == 3:
-                assert self.n_cols_diagnostics >= 7, (
-                    f"In case of the distribution '{self.f0}' with Particles6D, minimum number of n_cols_diagnostics is 7!"
-                )
-
-                self._f_coords_index = self.index["com"]["6D"]
-                self._f_jacobian_coords_index = self.index["pos+energy"]["6D"]
-
-            # Particles5D
-            elif self.vdim == 2:
-                assert self.n_cols_diagnostics >= 3, (
-                    f"In case of the distribution '{self.f0}' with Particles5D, minimum number of n_cols_diagnostics is 3!"
-                )
-
-                self._f_coords_index = self.index["com"]["5D"]
-                self._f_jacobian_coords_index = self.index["pos+energy"]["5D"]
-
         if self.type == "sph":
             self._f_coords_index = self.index["coords"]
             self._f_jacobian_coords_index = self.index["coords"]
-        else:
-            if self.f0.coords == "constants_of_motion":
-                self._f_coords_index = self.index["com"]
-                self._f_jacobian_coords_index = self.index["pos+energy"]
 
-            else:
-                self._f_coords_index = self.index["coords"]
-                self._f_jacobian_coords_index = self.index["coords"]
+        else:
+            self._f_coords_index = self.index["coords"]
+            self._f_jacobian_coords_index = self.index["coords"]
 
     def _n_mks_load_and_Np_per_clone(self):
         """Return two arrays: 1) an array of sub_comm.size where the i-th entry corresponds to the number of markers drawn on process i,
@@ -1810,6 +1789,9 @@ def initialize_weights(
         else:
             assert self.domain is not None, "A domain is needed to initialize weights."
 
+            if xp.size(self.markers_wo_holes_and_ghost) == 0:
+                return
+
             # set initial condition
             if bckgr_params is not None:
                 self._bckgr_params = bckgr_params
@@ -1820,6 +1802,9 @@ def initialize_weights(
             if self.type != "sph":
                 self._set_initial_condition()
 
+                if self.f_init.coords == "constants_of_motion":
+                    self.save_constants_of_motion()
+
             # evaluate initial distribution function
             if self.type == "sph":
                 f_init = self.f_init(self.positions)
@@ -1864,6 +1849,9 @@ def update_weights(self):
         The background :attr:`~struphy.pic.base.Particles.f0` is used for this.
         """
 
+        if xp.size(self.markers_wo_holes_and_ghost) == 0:
+            return
+
         if self.type == "sph":
             f0 = self.f0.n0(self.positions)
         else:
@@ -2223,7 +2211,7 @@ def gyro_transfer(self, outside_inds):
         assert xp.all(xp.isclose(v_perp, v - norm_b_cart * v_parallel))
 
         # calculate Larmor radius
-        Larmor_r = xp.cross(norm_b_cart, v_perp, axis=0) / absB0 * self._epsilon
+        Larmor_r = xp.cross(norm_b_cart, v_perp, axis=0) / absB0 * self.equation_params.epsilon
 
         # transform cartesian coordinates to logical coordinates
         # TODO: currently only possible with the geomoetry where its inverse map is defined.
diff --git a/src/struphy/pic/particles.py b/src/struphy/pic/particles.py
index 6c818b3ee..ba1ce2055 100644
--- a/src/struphy/pic/particles.py
+++ b/src/struphy/pic/particles.py
@@ -43,8 +43,8 @@ def __init__(
             kwargs["background"] = self.default_background()
 
         # default number of diagnostics and auxiliary columns
-        self._n_cols_diagnostics = kwargs.pop("n_cols_diagn", 0)
-        self._n_cols_aux = kwargs.pop("n_cols_aux", 5)
+        kwargs.setdefault("n_cols_diag", 0)
+        kwargs.setdefault("n_cols_aux", 5)
 
         super().__init__(**kwargs)
 
@@ -56,23 +56,12 @@ def __init__(
             self._absB0_h = self.projected_equil.absB0
             self._b2_h = self.projected_equil.b2
             self._derham = self.projected_equil.derham
-            self._epsilon = self.equation_params["epsilon"]
 
     @property
     def vdim(self):
         """Dimension of the velocity space."""
         return 3
 
-    @property
-    def n_cols_diagnostics(self):
-        """Number of the diagnostics columns."""
-        return self._n_cols_diagnostics
-
-    @property
-    def n_cols_aux(self):
-        """Number of the auxiliary columns."""
-        return self._n_cols_aux
-
     @property
     def coords(self):
         """Coordinates of the Particles6D, :math:`(v_1, v_2, v_3)`."""
@@ -156,7 +145,7 @@ def s0(self, eta1, eta2, eta3, *v, flat_eval=False, remove_holes=True):
 
     def save_constants_of_motion(self):
         """
-        Calculate each markers' guiding center constants of motions
+        Calculate each markers' guiding center constants of motion
         and assign them into diagnostics columns of marker array:
 
         ================= ============== ======= ============ ============= ==============
@@ -188,7 +177,7 @@ def save_constants_of_motion(self):
             self._derham.args_derham,
             self.domain.args_domain,
             self.first_diagnostics_idx,
-            self._epsilon,
+            self.equation_params.epsilon,
             self._b2_h[0]._data,
             self._b2_h[1]._data,
             self._b2_h[2]._data,
@@ -229,7 +218,7 @@ def save_constants_of_motion(self):
             self.markers,
             self._derham.args_derham,
             self.first_diagnostics_idx,
-            self._epsilon,
+            self.equation_params.epsilon,
             B0,
             R0,
             self._absB0_h._data,
@@ -299,24 +288,6 @@ class Particles5D(Particles):
     ===== ============== ========== ====== ======= ====== ====== ==========
     value position (eta) v_parallel v_perp  weight   s0     w0   buffer
     ===== ============== ========== ====== ======= ====== ====== ==========
-
-    Parameters
-    ----------
-    name : str
-        Name of particle species.
-
-    Np : int
-        Number of particles.
-
-    bc : list
-        Either 'remove', 'reflect', 'periodic' or 'refill' in each direction.
-
-    loading : str
-        Drawing of markers; either 'pseudo_random', 'sobol_standard',
-        'sobol_antithetic', 'external' or 'restart'.
-
-    **kwargs : dict
-        Parameters for markers, see :class:`~struphy.pic.base.Particles`.
     """
 
     @classmethod
@@ -336,8 +307,8 @@ def __init__(
         #     kwargs["bckgr_params"] = self.default_bckgr_params()
 
         # default number of diagnostics and auxiliary columns
-        self._n_cols_diagnostics = kwargs.pop("n_cols_diagn", 3)
-        self._n_cols_aux = kwargs.pop("n_cols_aux", 12)
+        kwargs.setdefault("n_cols_diag", 3)
+        kwargs.setdefault("n_cols_aux", 12)
 
         super().__init__(
             projected_equil=projected_equil,
@@ -347,7 +318,6 @@ def __init__(
         # magnetic background
         if self.equil is not None:
             assert isinstance(self.equil, FluidEquilibriumWithB), "Particles5D needs background with magnetic field."
-        self._magn_bckgr = self.equil
 
         self._absB0_h = self.projected_equil.absB0
         self._unit_b1_h = self.projected_equil.unit_b1
@@ -361,21 +331,6 @@ def vdim(self):
         """Dimension of the velocity space."""
         return 2
 
-    @property
-    def n_cols_diagnostics(self):
-        """Number of the diagnostics columns."""
-        return self._n_cols_diagnostics
-
-    @property
-    def n_cols_aux(self):
-        """Number of the auxiliary columns."""
-        return self._n_cols_aux
-
-    @property
-    def magn_bckgr(self):
-        """Fluid equilibrium with B."""
-        return self._magn_bckgr
-
     @property
     def absB0_h(self):
         """Discrete 0-form coefficients of |B_0|."""
@@ -386,11 +341,6 @@ def unit_b1_h(self):
         """Discrete 1-form coefficients of B/|B|."""
         return self._unit_b1_h
 
-    @property
-    def epsilon(self):
-        """One of equation params, epsilon"""
-        return self._epsilon
-
     @property
     def coords(self):
         r"""Coordinates of the Particles5D, :math:`(v_\parallel, \mu)`."""
@@ -421,22 +371,14 @@ def svol(self, eta1, eta2, eta3, *v):
         """
         # load sampling density svol (normalized to 1 in logical space)
         maxw_params = {
-            "n": 1.0,
-            "u_para": self.loading_params.moments[0],
-            "u_perp": self.loading_params.moments[1],
-            "vth_para": self.loading_params.moments[2],
-            "vth_perp": self.loading_params.moments[3],
+            "n": (1.0, None),
+            "u_para": (self.loading_params.moments[0], None),
+            "u_perp": (self.loading_params.moments[1], None),
+            "vth_para": (self.loading_params.moments[2], None),
+            "vth_perp": (self.loading_params.moments[3], None),
         }
 
-        self._svol = maxwellians.GyroMaxwellian2D(
-            n=(1.0, None),
-            u_para=(self.loading_params.moments[0], None),
-            u_perp=(self.loading_params.moments[1], None),
-            vth_para=(self.loading_params.moments[2], None),
-            vth_perp=(self.loading_params.moments[3], None),
-            volume_form=True,
-            equil=self._magn_bckgr,
-        )
+        self._svol = maxwellians.GyroMaxwellian2D(**maxw_params, equil=self.equil)
 
         if self.spatial == "uniform":
             out = self._svol(eta1, eta2, eta3, *v)
@@ -551,13 +493,11 @@ def save_constants_of_motion(self):
         r = self.markers[~self.holes, 0] * (1 - a1) + a1
         self.markers[~self.holes, idx_can_momentum] = self.equil.psi_r(r)
 
-        self._epsilon = self.equation_params["epsilon"]
-
         utilities_kernels.eval_canonical_toroidal_moment_5d(
             self.markers,
             self.derham.args_derham,
             self.first_diagnostics_idx,
-            self.epsilon,
+            self.equation_params.epsilon,
             B0,
             R0,
             self.absB0_h._data,
@@ -661,8 +601,8 @@ def __init__(
             kwargs["background"] = self.default_background()
 
         # default number of diagnostics and auxiliary columns
-        self._n_cols_diagnostics = kwargs.pop("n_cols_diagn", 0)
-        self._n_cols_aux = kwargs.pop("n_cols_aux", 5)
+        kwargs.setdefault("n_cols_diag", 0)
+        kwargs.setdefault("n_cols_aux", 5)
 
         super().__init__(**kwargs)
 
@@ -671,16 +611,6 @@ def vdim(self):
         """Dimension of the velocity space."""
         return 0
 
-    @property
-    def n_cols_diagnostics(self):
-        """Number of the diagnostics columns."""
-        return self._n_cols_diagnostics
-
-    @property
-    def n_cols_aux(self):
-        """Number of the auxiliary columns."""
-        return self._n_cols_aux
-
     @property
     def coords(self):
         """Coordinates of the Particles3D."""
@@ -802,8 +732,8 @@ def __init__(
         #         kwargs["sorting_params"]["communicate"] = True
 
         # default number of diagnostics and auxiliary columns
-        self._n_cols_diagnostics = kwargs.pop("n_cols_diagn", 0)
-        self._n_cols_aux = kwargs.pop("n_cols_aux", 24)
+        kwargs.setdefault("n_cols_diag", 0)
+        kwargs.setdefault("n_cols_aux", 24)
 
         clone_config = kwargs.get("clone_config", None)
         assert clone_config is None, "SPH can only be launched with --nclones 1"
@@ -815,16 +745,6 @@ def vdim(self):
         """Dimension of the velocity space."""
         return 3
 
-    @property
-    def n_cols_diagnostics(self):
-        """Number of the diagnostics columns."""
-        return self._n_cols_diagnostics
-
-    @property
-    def n_cols_aux(self):
-        """Number of the auxiliary columns."""
-        return self._n_cols_aux
-
     @property
     def coords(self):
         """Coordinates of the Particles6D, :math:`(v_1, v_2, v_3)`."""
diff --git a/src/struphy/pic/pushing/pusher_kernels_gc.py b/src/struphy/pic/pushing/pusher_kernels_gc.py
index 6e5df8034..d10157f2e 100644
--- a/src/struphy/pic/pushing/pusher_kernels_gc.py
+++ b/src/struphy/pic/pushing/pusher_kernels_gc.py
@@ -2470,7 +2470,11 @@ def push_gc_cc_J2_stage_Hdiv(
         last = 0.0
 
     for ip in range(n_markers):
-        # check if marker is a hole
+        # only do something if particle is a "true" particle (i.e. not a hole)
+        if markers[ip, 0] == -1.0:
+            continue
+
+        # if particle is refilled
         if markers[ip, first_init_idx] == -1.0:
             continue
 
diff --git a/src/struphy/propagators/propagators_coupling.py b/src/struphy/propagators/propagators_coupling.py
index a4bd37925..a098d15bb 100644
--- a/src/struphy/propagators/propagators_coupling.py
+++ b/src/struphy/propagators/propagators_coupling.py
@@ -1291,7 +1291,7 @@ def allocate(self):
 
         # magnetic equilibrium field
         unit_b1 = self.projected_equil.unit_b1
-        curl_unit_b1 = self.projected_equil.curl_unit_b1
+        curl_unit_b2 = self.projected_equil.curl_unit_b2
         self._b2 = self.projected_equil.b2
 
         # magnetic field
@@ -1326,9 +1326,9 @@ def allocate(self):
             unit_b1[0]._data,
             unit_b1[1]._data,
             unit_b1[2]._data,
-            curl_unit_b1[0]._data,
-            curl_unit_b1[1]._data,
-            curl_unit_b1[2]._data,
+            curl_unit_b2[0]._data,
+            curl_unit_b2[1]._data,
+            curl_unit_b2[2]._data,
             self._u_form_int,
         )
 
@@ -1341,7 +1341,7 @@ def allocate(self):
             pusher_kernel = Pyccelkernel(pusher_kernels_gc.push_gc_cc_J1_H1vec)
         else:
             raise ValueError(
-                f'{self.options.u_space  =} not valid, choose from "Hcurl", "Hdiv" or "H1vec.',
+                f'{self.options.u_space  = } not valid, choose from "Hcurl", "Hdiv" or "H1vec.',
             )
 
         args_pusher_kernel = (
@@ -1353,9 +1353,9 @@ def allocate(self):
             unit_b1[0]._data,
             unit_b1[1]._data,
             unit_b1[2]._data,
-            curl_unit_b1[0]._data,
-            curl_unit_b1[1]._data,
-            curl_unit_b1[2]._data,
+            curl_unit_b2[0]._data,
+            curl_unit_b2[1]._data,
+            curl_unit_b2[2]._data,
             self._u_avg[0]._data,
             self._u_avg[1]._data,
             self._u_avg[2]._data,
@@ -1571,10 +1571,11 @@ def allocate(self):
             tol=self.options.solver_params.tol,
             maxiter=self.options.solver_params.maxiter,
             verbose=self.options.solver_params.verbose,
+            recycle=self.options.solver_params.recycle,
         )
         # magnetic equilibrium field
         unit_b1 = self.projected_equil.unit_b1
-        curl_unit_b1 = self.projected_equil.curl_unit_b1
+        curl_unit_b2 = self.projected_equil.curl_unit_b2
         self._b2 = self.projected_equil.b2
         gradB1 = self.projected_equil.gradB1
         absB0 = self.projected_equil.absB0
@@ -1615,12 +1616,15 @@ def allocate(self):
                 unit_b1[0]._data,
                 unit_b1[1]._data,
                 unit_b1[2]._data,
-                curl_unit_b1[0]._data,
-                curl_unit_b1[1]._data,
-                curl_unit_b1[2]._data,
+                curl_unit_b2[0]._data,
+                curl_unit_b2[1]._data,
+                curl_unit_b2[2]._data,
                 self._grad_PB_b[0]._data,
                 self._grad_PB_b[1]._data,
                 self._grad_PB_b[2]._data,
+                gradB1[0]._data,
+                gradB1[1]._data,
+                gradB1[2]._data,
                 self._u_form_int,
             )
 
@@ -1631,7 +1635,7 @@ def allocate(self):
                 self._pusher_kernel = pusher_kernels_gc.push_gc_cc_J2_stage_H1vec
             else:
                 raise ValueError(
-                    f'{self.options.u_space  =} not valid, choose from "Hdiv" or "H1vec.',
+                    f'{self.options.u_space  = } not valid, choose from "Hdiv" or "H1vec.',
                 )
 
             # temp fix due to refactoring of ButcherTableau:
@@ -1651,9 +1655,9 @@ def allocate(self):
                 unit_b1[0]._data,
                 unit_b1[1]._data,
                 unit_b1[2]._data,
-                curl_unit_b1[0]._data,
-                curl_unit_b1[1]._data,
-                curl_unit_b1[2]._data,
+                curl_unit_b2[0]._data,
+                curl_unit_b2[1]._data,
+                curl_unit_b2[2]._data,
                 self._u_temp[0]._data,
                 self._u_temp[1]._data,
                 self._u_temp[2]._data,
@@ -1692,9 +1696,9 @@ def allocate(self):
                 unit_b1[0]._data,
                 unit_b1[1]._data,
                 unit_b1[2]._data,
-                curl_unit_b1[0]._data,
-                curl_unit_b1[1]._data,
-                curl_unit_b1[2]._data,
+                curl_unit_b2[0]._data,
+                curl_unit_b2[1]._data,
+                curl_unit_b2[2]._data,
                 self._grad_PB_b[0]._data,
                 self._grad_PB_b[1]._data,
                 self._grad_PB_b[2]._data,
@@ -1743,9 +1747,9 @@ def allocate(self):
                 unit_b1[0]._data,
                 unit_b1[1]._data,
                 unit_b1[2]._data,
-                curl_unit_b1[0]._data,
-                curl_unit_b1[1]._data,
-                curl_unit_b1[2]._data,
+                curl_unit_b2[0]._data,
+                curl_unit_b2[1]._data,
+                curl_unit_b2[2]._data,
                 self.variables.u.spline.vector[0]._data,
                 self.variables.u.spline.vector[1]._data,
                 self.variables.u.spline.vector[2]._data,
@@ -1761,9 +1765,9 @@ def allocate(self):
                 unit_b1[0]._data,
                 unit_b1[1]._data,
                 unit_b1[2]._data,
-                curl_unit_b1[0]._data,
-                curl_unit_b1[1]._data,
-                curl_unit_b1[2]._data,
+                curl_unit_b2[0]._data,
+                curl_unit_b2[1]._data,
+                curl_unit_b2[2]._data,
                 self._u_mid[0]._data,
                 self._u_mid[1]._data,
                 self._u_mid[2]._data,
@@ -1800,12 +1804,13 @@ def __call__(self, dt):
         b_full.update_ghost_regions()
 
         if self.options.algo == "explicit":
-            PB_b = self._PB.dot(b_full, out=self._PB_b)
+            PB_b = self._PB.dot(self._b_tilde, out=self._PB_b)
             grad_PB_b = self.derham.grad.dot(PB_b, out=self._grad_PB_b)
             grad_PB_b.update_ghost_regions()
 
             # save old u
             u_new = un.copy(out=self._u_new)
+            u_temp = un.copy(out=self._u_temp)
 
             for stage in range(self.options.butcher.n_stages):
                 # accumulate
@@ -1839,6 +1844,8 @@ def __call__(self, dt):
                 # calculate u^{n+1}
                 u_new += ku * dt * self.options.butcher.b[stage]
 
+                u_new.update_ghost_regions()
+
                 if self.options.solver_params.info and MPI.COMM_WORLD.Get_rank() == 0:
                     print("Stage: ", stage)
                     print("Status     for CurrentCoupling5DGradB:", info["success"])
@@ -1978,7 +1985,7 @@ def __call__(self, dt):
                 sum_u_diff_loc = xp.sum((u_diff.toarray() ** 2))
 
                 sum_H_diff_loc = xp.sum(
-                    (markers[~holes, :3] - markers[~holes, first_init_idx : first_init_idx + 3]) ** 2,
+                    (markers[~holes, :3] - markers[~holes, first_init_idx : first_init_idx + 3]) ** 2
                 )
 
                 buffer_array = xp.array([sum_u_diff_loc])
@@ -2072,7 +2079,7 @@ def __call__(self, dt):
                 )
 
                 sum_H_diff_loc = xp.sum(
-                    xp.abs(markers[~holes, 0:3] - markers[~holes, first_free_idx : first_free_idx + 3]),
+                    xp.abs(markers[~holes, 0:3] - markers[~holes, first_free_idx : first_free_idx + 3])
                 )
 
                 if particles.mpi_comm is not None:
diff --git a/src/struphy/propagators/propagators_fields.py b/src/struphy/propagators/propagators_fields.py
index a54aebd47..cd20f2ac5 100644
--- a/src/struphy/propagators/propagators_fields.py
+++ b/src/struphy/propagators/propagators_fields.py
@@ -2166,6 +2166,10 @@ def __call__(self, dt):
         else:
             self._ode_solver(0.0, dt)
 
+        # update_weights
+        if self.options.energetic_ions.species.weights_params.control_variate:
+            self.options.energetic_ions.particles.update_weights()
+
         if self._info and MPI.COMM_WORLD.Get_rank() == 0:
             if self.options.algo == "implicit":
                 print("Status     for ShearAlfvenCurrentCoupling5D:", info["success"])
@@ -2356,7 +2360,7 @@ class Options:
         b_tilde: FEECVariable = None
         ep_scale: float = 1.0
         u_space: OptsVecSpace = "Hdiv"
-        solver: OptsSymmSolver = "pcg"
+        solver: OptsGenSolver = "pbicgstab"
         precond: OptsMassPrecond = "MassMatrixPreconditioner"
         solver_params: SolverParameters = None
         filter_params: FilterParameters = None
@@ -2364,7 +2368,7 @@ class Options:
         def __post_init__(self):
             # checks
             check_option(self.u_space, OptsVecSpace)
-            check_option(self.solver, OptsSymmSolver)
+            check_option(self.solver, OptsGenSolver)
             check_option(self.precond, OptsMassPrecond)
             assert isinstance(self.energetic_ions, PICVariable)
             assert self.energetic_ions.space == "Particles5D"
@@ -2406,7 +2410,7 @@ def allocate(self):
 
         # magnetic equilibrium field
         unit_b1 = self.projected_equil.unit_b1
-        curl_unit_b1 = self.projected_equil.curl_unit_b1
+        curl_unit_b2 = self.projected_equil.curl_unit_b2
         self._b2 = self.projected_equil.b2
 
         # scaling factor
@@ -2438,9 +2442,9 @@ def allocate(self):
             unit_b1[0]._data,
             unit_b1[1]._data,
             unit_b1[2]._data,
-            curl_unit_b1[0]._data,
-            curl_unit_b1[1]._data,
-            curl_unit_b1[2]._data,
+            curl_unit_b2[0]._data,
+            curl_unit_b2[1]._data,
+            curl_unit_b2[2]._data,
             self._u_form_int,
         )