diff --git a/pyphare/pyphare/pharein/diagnostics.py b/pyphare/pyphare/pharein/diagnostics.py
index 2721c9633..3bd164232 100644
--- a/pyphare/pyphare/pharein/diagnostics.py
+++ b/pyphare/pyphare/pharein/diagnostics.py
@@ -222,6 +222,7 @@ def population_in_model(population):
 class FluidDiagnostics_(Diagnostics):
     fluid_quantities = [
         "density",
+        "charge_density",
         "mass_density",
         "flux",
         "bulkVelocity",
diff --git a/pyphare/pyphare/pharesee/hierarchy/__init__.py b/pyphare/pyphare/pharesee/hierarchy/__init__.py
index 6bc7e43c8..280428513 100644
--- a/pyphare/pyphare/pharesee/hierarchy/__init__.py
+++ b/pyphare/pyphare/pharesee/hierarchy/__init__.py
@@ -11,10 +11,11 @@
 
 
 def hierarchy_from(
-    simulator=None, qty=None, pop="", h5_filename=None, times=None, hier=None, **kwargs
+    simulator=None, qty=None, pop="", h5_filename=None, times=None, hier=None, func=None, **kwargs
 ):
     from .fromh5 import hierarchy_fromh5
     from .fromsim import hierarchy_from_sim
+    from .fromfunc import hierarchy_from_func
 
     """
     this function reads an HDF5 PHARE file and returns a PatchHierarchy from
@@ -39,4 +40,9 @@ def hierarchy_from(
     if simulator is not None and qty is not None:
         return hierarchy_from_sim(simulator, qty, pop=pop)
 
+    if func is not None and hier is not None:
+        if not callable(func):
+            raise TypeError("func must be callable")
+        return hierarchy_from_func(func, hier, **kwargs)
+
     raise ValueError("can't make hierarchy")
diff --git a/pyphare/pyphare/pharesee/hierarchy/fromfunc.py b/pyphare/pyphare/pharesee/hierarchy/fromfunc.py
new file mode 100644
index 000000000..ac313c298
--- /dev/null
+++ b/pyphare/pyphare/pharesee/hierarchy/fromfunc.py
@@ -0,0 +1,63 @@
+from pyphare.pharesee.hierarchy.hierarchy_utils import compute_hier_from
+
+
+def hierarchy_from_func1d(func, hier, **kwargs):
+    assert hier.ndim == 1
+
+    def compute_(patch_datas, **kwargs):
+        ref_name = next(iter(patch_datas.keys()))
+        x_ = patch_datas[ref_name].x
+
+        return (
+            {"name": "value", "data": func(x_, **kwargs), "centering": patch_datas[ref_name].centerings},
+        )
+
+    return compute_hier_from(compute_, hier, **kwargs)
+
+
+def hierarchy_from_func2d(func, hier, **kwargs):
+    from pyphare.pharesee.hierarchy.hierarchy_utils import compute_hier_from
+
+    def compute_(patch_datas, **kwargs):
+        ref_name = next(iter(patch_datas.keys()))
+        x_ = patch_datas[ref_name].x
+        y_ = patch_datas[ref_name].y
+
+        return (
+            {"name": "value", "data": func(x_, y_, **kwargs), "centering": patch_datas[ref_name].centerings},
+        )
+
+    return compute_hier_from(compute_, hier, **kwargs)
+
+
+def hierarchy_from_func(func, hier, **kwargs):
+
+    """
+    Route hierarchical computation to appropriate dimension handler.
+
+    Parameters
+    ----------
+    func : callable
+        Function to apply to coordinates of the hierarchy
+    hier : Hierarchy
+        Hierarchy object (1D or 2D)
+    **kwargs : dict
+        Additional arguments passed to func
+
+    Returns
+    -------
+    dict
+        Computed hierarchical data
+
+    Raises
+    ------
+    ValueError
+        If hierarchy dimension is not supported
+    """
+
+    if hier.ndim == 1:
+        return hierarchy_from_func1d(func, hier, **kwargs)
+    if hier.ndim == 2:
+        return hierarchy_from_func2d(func, hier, **kwargs)
+    else:
+        raise ValueError(f"Unsupported hierarchy dimension: {hier.ndim}")
diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py
index cab07023c..ea575368e 100644
--- a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py
+++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py
@@ -430,8 +430,7 @@ def plot1d(self, **kwargs):
                 if qty is None:
                     qty = pdata_names[0]
 
-                layout = patch.patch_datas[qty].layout
-                nbrGhosts = layout.nbrGhostFor(qty)
+                nbrGhosts = patch.patch_datas[qty].ghosts_nbr
                 val = patch.patch_datas[qty][patch.box]
                 x = patch.patch_datas[qty].x[nbrGhosts[0] : -nbrGhosts[0]]
                 label = "L{level}P{patch}".format(level=lvl_nbr, patch=ip)
@@ -542,9 +541,10 @@ def plot2d(self, **kwargs):
         if "ylim" in kwargs:
             ax.set_ylim(kwargs["ylim"])
 
-        divider = make_axes_locatable(ax)
-        cax = divider.append_axes("right", size="5%", pad=0.08)
-        plt.colorbar(im, ax=ax, cax=cax)
+        if kwargs.get("cbar", True):
+            divider = make_axes_locatable(ax)
+            cax = divider.append_axes("right", size="5%", pad=0.08)
+            fig.colorbar(im, ax=ax, cax=cax)
 
         if kwargs.get("legend", None) is not None:
             ax.legend()
diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py
index 487dc8605..b7fbcc220 100644
--- a/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py
+++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py
@@ -36,6 +36,7 @@
     "momentum_tensor_zz": "Mzz",
     "density": "rho",
     "mass_density": "rho",
+    "charge_density": "rho",
     "tags": "tags",
 }
 
diff --git a/pyphare/pyphare/pharesee/plotting.py b/pyphare/pyphare/pharesee/plotting.py
index 180db4375..8f7614919 100644
--- a/pyphare/pyphare/pharesee/plotting.py
+++ b/pyphare/pyphare/pharesee/plotting.py
@@ -274,7 +274,7 @@ def finest_field_plot(run_path, qty, **kwargs):
             time = times[0]
         interpolator, finest_coords = r.GetVi(time, merged=True, interp=interp)[qty]
     elif qty == "rho":
-        file = os.path.join(run_path, "ions_density.h5")
+        file = os.path.join(run_path, "ions_charge_density.h5")
         if time is None:
             times = get_times_from_h5(file)
             time = times[0]
diff --git a/pyphare/pyphare/pharesee/run/run.py b/pyphare/pyphare/pharesee/run/run.py
index 003943497..30f6a25f8 100644
--- a/pyphare/pyphare/pharesee/run/run.py
+++ b/pyphare/pyphare/pharesee/run/run.py
@@ -27,7 +27,7 @@
     "EM_B": "B",
     "EM_E": "E",
     "ions_bulkVelocity": "Vi",
-    "ions_density": "Ni",
+    "ions_charge_density": "Ni",
     "particle_count": "nppc",
 }
 
@@ -114,7 +114,7 @@ def GetMassDensity(self, time, merged=False, interp="nearest", **kwargs):
         return ScalarField(self._get(hier, time, merged, interp))
 
     def GetNi(self, time, merged=False, interp="nearest", **kwargs):
-        hier = self._get_hierarchy(time, "ions_density.h5", **kwargs)
+        hier = self._get_hierarchy(time, "ions_charge_density.h5", **kwargs)
         return ScalarField(self._get(hier, time, merged, interp))
 
     def GetN(self, time, pop_name, merged=False, interp="nearest", **kwargs):
@@ -151,7 +151,7 @@ def GetPi(self, time, merged=False, interp="nearest", **kwargs):
         return self._get(Pi, time, merged, interp)  # should later be a TensorField
 
     def GetPe(self, time, merged=False, interp="nearest", all_primal=True):
-        hier = self._get_hierarchy(time, "ions_density.h5")
+        hier = self._get_hierarchy(time, "ions_charge_density.h5")
 
         Te = hier.sim.electrons.closure.Te
 
diff --git a/pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py b/pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py
index 29a8249e3..c31bb3095 100644
--- a/pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py
+++ b/pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py
@@ -177,11 +177,11 @@ def test_large_patchoverlaps(self, expected):
             collections=[
                 {
                     "boxes": overlap_boxes,
-                    "facecolor": "yellow",
+                    "facecolor": "grey",
                 },
                 {
                     "boxes": [p.box for p in hierarchy.level(ilvl).patches],
-                    "facecolor": "grey",
+                    "facecolor": "yellow",
                 },
             ],
         )
diff --git a/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py b/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py
index f65bb39f0..908484607 100644
--- a/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py
+++ b/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py
@@ -145,7 +145,7 @@ def vthz(x, y):
             for quantity in ["E", "B"]:
                 ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
-            for quantity in ["density", "bulkVelocity"]:
+            for quantity in ["charge_density", "bulkVelocity"]:
                 ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
             return sim
diff --git a/src/amr/level_initializer/hybrid_level_initializer.hpp b/src/amr/level_initializer/hybrid_level_initializer.hpp
index e0ef8386a..3e112a34e 100644
--- a/src/amr/level_initializer/hybrid_level_initializer.hpp
+++ b/src/amr/level_initializer/hybrid_level_initializer.hpp
@@ -98,7 +98,7 @@ namespace solver
                     core::depositParticles(ions, layout, interpolate_, core::LevelGhostDeposit{});
                 }
 
-                ions.computeDensity();
+                ions.computeChargeDensity();
                 ions.computeBulkVelocity();
             }
 
diff --git a/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp b/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp
index 3b6a90393..e47bc5779 100644
--- a/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp
+++ b/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp
@@ -385,20 +385,21 @@ namespace amr
                 {
                     // first thing to do is to project patchGhostParitcles moments
                     auto& patchGhosts = pop.patchGhostParticles();
-                    auto& density     = pop.density();
+                    auto& particleDensity     = pop.particleDensity();
+                    auto& chargeDensity     = pop.chargeDensity();
                     auto& flux        = pop.flux();
 
-                    interpolate_(makeRange(patchGhosts), density, flux, layout);
+                    interpolate_(makeRange(patchGhosts), particleDensity, chargeDensity, flux, layout);
 
                     if (level.getLevelNumber() > 0) // no levelGhost on root level
                     {
                         // then grab levelGhostParticlesOld and levelGhostParticlesNew
                         // and project them with alpha and (1-alpha) coefs, respectively
                         auto& levelGhostOld = pop.levelGhostParticlesOld();
-                        interpolate_(makeRange(levelGhostOld), density, flux, layout, 1. - alpha);
+                        interpolate_(makeRange(levelGhostOld), particleDensity, chargeDensity, flux, layout, 1. - alpha);
 
                         auto& levelGhostNew = pop.levelGhostParticlesNew();
-                        interpolate_(makeRange(levelGhostNew), density, flux, layout, alpha);
+                        interpolate_(makeRange(levelGhostNew), particleDensity, chargeDensity, flux, layout, alpha);
                     }
                 }
             }
@@ -538,7 +539,7 @@ namespace amr
 
                 auto& J  = hybridModel.state.J;
                 auto& Vi = hybridModel.state.ions.velocity();
-                auto& Ni = hybridModel.state.ions.density();
+                auto& Ni = hybridModel.state.ions.chargeDensity();
 
                 Jold_.copyData(J);
                 ViOld_.copyData(Vi);
diff --git a/src/amr/physical_models/hybrid_model.hpp b/src/amr/physical_models/hybrid_model.hpp
index 449bafc22..6076836f5 100644
--- a/src/amr/physical_models/hybrid_model.hpp
+++ b/src/amr/physical_models/hybrid_model.hpp
@@ -154,7 +154,7 @@ void HybridModel<GridLayoutT, Electromag, Ions, Electrons, AMR_Types, Grid_t>::f
 
     hybridInfo.modelMagnetic        = core::VecFieldNames{state.electromag.B};
     hybridInfo.modelElectric        = core::VecFieldNames{state.electromag.E};
-    hybridInfo.modelIonDensity      = state.ions.densityName();
+    hybridInfo.modelIonDensity      = state.ions.chargeDensityName();
     hybridInfo.modelIonBulkVelocity = core::VecFieldNames{state.ions.velocity()};
     hybridInfo.modelCurrent         = core::VecFieldNames{state.J};
 
diff --git a/src/amr/tagging/default_hybrid_tagger_strategy.hpp b/src/amr/tagging/default_hybrid_tagger_strategy.hpp
index b9a016f7d..1bf22f723 100644
--- a/src/amr/tagging/default_hybrid_tagger_strategy.hpp
+++ b/src/amr/tagging/default_hybrid_tagger_strategy.hpp
@@ -36,7 +36,7 @@ void DefaultHybridTaggerStrategy<HybridModel>::tag(HybridModel& model,
     auto& By = model.state.electromag.B.getComponent(PHARE::core::Component::Y);
     auto& Bz = model.state.electromag.B.getComponent(PHARE::core::Component::Z);
 
-    auto& N = model.state.ions.density();
+    // auto& N = model.state.ions.chargeDensity();
 
     // we loop on cell indexes for all qties regardless of their centering
     auto const& [start_x, _]
diff --git a/src/core/data/electrons/electrons.hpp b/src/core/data/electrons/electrons.hpp
index a079022da..b18d0c290 100644
--- a/src/core/data/electrons/electrons.hpp
+++ b/src/core/data/electrons/electrons.hpp
@@ -65,7 +65,7 @@ class StandardHybridElectronFluxComputer
     {
         if (isUsable())
         {
-            return ions_.density();
+            return ions_.chargeDensity();
         }
         else
         {
@@ -78,7 +78,7 @@ class StandardHybridElectronFluxComputer
     {
         if (isUsable())
         {
-            return ions_.density();
+            return ions_.chargeDensity();
         }
         else
         {
@@ -112,23 +112,23 @@ class StandardHybridElectronFluxComputer
         auto const& Vix = ions_.velocity()(Component::X);
         auto const& Viy = ions_.velocity()(Component::Y);
         auto const& Viz = ions_.velocity()(Component::Z);
-        auto const& Ni  = ions_.density();
+        auto const& Ne  = ions_.chargeDensity();
 
         auto& Vex = Ve_(Component::X);
         auto& Vey = Ve_(Component::Y);
         auto& Vez = Ve_(Component::Z);
 
         // from Ni because all components defined on primal
-        layout.evalOnBox(Ni, [&](auto const&... args) {
+        layout.evalOnBox(Ne, [&](auto const&... args) {
             auto arr = std::array{args...};
 
             auto const JxOnVx = GridLayout::project(Jx, arr, GridLayout::JxToMoments());
             auto const JyOnVy = GridLayout::project(Jy, arr, GridLayout::JyToMoments());
             auto const JzOnVz = GridLayout::project(Jz, arr, GridLayout::JzToMoments());
 
-            Vex(arr) = Vix(arr) - JxOnVx / Ni(arr);
-            Vey(arr) = Viy(arr) - JyOnVy / Ni(arr);
-            Vez(arr) = Viz(arr) - JzOnVz / Ni(arr);
+            Vex(arr) = Vix(arr) - JxOnVx / Ne(arr);
+            Vey(arr) = Viy(arr) - JyOnVy / Ne(arr);
+            Vez(arr) = Viz(arr) - JzOnVz / Ne(arr);
         });
     }
 
@@ -211,7 +211,7 @@ class IsothermalElectronPressureClosure
         if (!Pe_.isUsable())
             throw std::runtime_error("Error - isothermal closure pressure not usable");
 
-        auto const& Ne_ = ions_.density();
+        auto const& Ne_ = ions_.chargeDensity();
         std::transform(std::begin(Ne_), std::end(Ne_), std::begin(Pe_),
                        [this](auto n) { return n * Te_; });
     }
diff --git a/src/core/data/ions/ion_population/ion_population.hpp b/src/core/data/ions/ion_population/ion_population.hpp
index aa6554e47..33dfea9d8 100644
--- a/src/core/data/ions/ion_population/ion_population.hpp
+++ b/src/core/data/ions/ion_population/ion_population.hpp
@@ -31,7 +31,8 @@ namespace core
             , mass_{initializer["mass"].template to<double>()}
             , flux_{name_ + "_flux", HybridQuantity::Vector::V}
             , momentumTensor_{name_ + "_momentumTensor", HybridQuantity::Tensor::M}
-            , rho_{name_ + "_rho", HybridQuantity::Scalar::rho}
+            , particleDensity_{name_ + "_particleDensity", HybridQuantity::Scalar::rho}
+            , chargeDensity_{name_ + "_chargeDensity", HybridQuantity::Scalar::rho}
             , particles_{name_}
             , particleInitializerInfo_{initializer["particle_initializer"]}
         {
@@ -42,20 +43,16 @@ namespace core
 
         NO_DISCARD std::string const& name() const { return name_; }
 
-
         NO_DISCARD auto const& particleInitializerInfo() const { return particleInitializerInfo_; }
 
-
-
         NO_DISCARD bool isUsable() const
         {
-            return core::isUsable(particles_, rho_, flux_, momentumTensor_);
+            return core::isUsable(particles_, particleDensity_, chargeDensity_, flux_, momentumTensor_);
         }
 
-
         NO_DISCARD bool isSettable() const
         {
-            return core::isSettable(particles_, rho_, flux_, momentumTensor_);
+            return core::isSettable(particles_, particleDensity_, chargeDensity_, flux_, momentumTensor_);
         }
 
         NO_DISCARD auto& domainParticles() const { return particles_.domainParticles(); }
@@ -79,9 +76,11 @@ namespace core
             return particles_.levelGhostParticlesNew();
         }
 
+        NO_DISCARD field_type const& particleDensity() const { return particleDensity_; }
+        NO_DISCARD field_type& particleDensity() { return particleDensity_; }
 
-        NO_DISCARD field_type const& density() const { return rho_; }
-        NO_DISCARD field_type& density() { return rho_; }
+        NO_DISCARD field_type const& chargeDensity() const { return chargeDensity_; }
+        NO_DISCARD field_type& chargeDensity() { return chargeDensity_; }
 
         NO_DISCARD VecField const& flux() const { return flux_; }
         NO_DISCARD VecField& flux() { return flux_; }
@@ -100,7 +99,8 @@ namespace core
 
         NO_DISCARD auto getCompileTimeResourcesViewList()
         {
-            return std::forward_as_tuple(flux_, momentumTensor_, rho_, particles_);
+            return std::forward_as_tuple(flux_, momentumTensor_, particleDensity_, chargeDensity_,
+                                         particles_);
         }
 
 
@@ -124,7 +124,8 @@ namespace core
         double mass_;
         VecField flux_;
         TensorField momentumTensor_;
-        field_type rho_;
+        field_type particleDensity_;
+        field_type chargeDensity_;
         ParticlesPack<ParticleArray> particles_;
         initializer::PHAREDict const& particleInitializerInfo_;
     };
diff --git a/src/core/data/ions/ions.hpp b/src/core/data/ions/ions.hpp
index 2701ebdac..a9a536802 100644
--- a/src/core/data/ions/ions.hpp
+++ b/src/core/data/ions/ions.hpp
@@ -41,15 +41,14 @@ namespace core
 
 
         explicit Ions(PHARE::initializer::PHAREDict const& dict)
-            : rho_{densityName(), HybridQuantity::Scalar::rho}
-            , massDensity_{massDensityName(), HybridQuantity::Scalar::rho}
+            : massDensity_{massDensityName(), HybridQuantity::Scalar::rho}
+            , chargeDensity_{chargeDensityName(), HybridQuantity::Scalar::rho}
             , bulkVelocity_{"bulkVel", HybridQuantity::Vector::V}
             , populations_{generate(
                   [&dict](auto ipop) { //
                       return IonPopulation{dict["pop" + std::to_string(ipop)]};
                   },
                   dict["nbrPopulations"].template to<std::size_t>())}
-            , sameMasses_{allSameMass_()}
             , momentumTensor_{"momentumTensor", HybridQuantity::Tensor::M}
         {
         }
@@ -58,29 +57,24 @@ namespace core
         NO_DISCARD auto nbrPopulations() const { return populations_.size(); }
         NO_DISCARD auto size() const { return nbrPopulations(); }
 
+        NO_DISCARD field_type const& massDensity() const { return massDensity_; }
+        NO_DISCARD field_type const& massDensity() { return massDensity_; }
 
-        NO_DISCARD field_type const& density() const { return rho_; }
-        NO_DISCARD field_type& density() { return rho_; }
-
-        NO_DISCARD field_type const& massDensity() const
-        {
-            return sameMasses_ ? rho_ : massDensity_;
-        }
-        NO_DISCARD field_type& massDensity() { return sameMasses_ ? rho_ : massDensity_; }
-
+        NO_DISCARD field_type const& chargeDensity() const { return chargeDensity_; }
+        NO_DISCARD field_type& chargeDensity() { return chargeDensity_; }
 
         NO_DISCARD vecfield_type const& velocity() const { return bulkVelocity_; }
         NO_DISCARD vecfield_type& velocity() { return bulkVelocity_; }
 
-        NO_DISCARD std::string static densityName() { return "rho"; }
+        NO_DISCARD std::string static chargeDensityName() { return "chargeDensity"; }
         NO_DISCARD std::string static massDensityName() { return "massDensity"; }
 
         tensorfield_type const& momentumTensor() const { return momentumTensor_; }
         tensorfield_type& momentumTensor() { return momentumTensor_; }
 
-        void computeDensity()
+        void computeChargeDensity()
         {
-            rho_.zero();
+            chargeDensity_.zero();
 
             for (auto const& pop : populations_)
             {
@@ -88,11 +82,13 @@ namespace core
                 // nodes. This is more efficient and easier to code as we don't
                 // have to account for the field dimensionality.
 
-                auto& popDensity = pop.density();
-                std::transform(std::begin(rho_), std::end(rho_), std::begin(popDensity),
-                               std::begin(rho_), std::plus<Float>{});
+                auto& popDensity = pop.chargeDensity();
+                std::transform(std::begin(chargeDensity_), std::end(chargeDensity_),
+                               std::begin(popDensity), std::begin(chargeDensity_),
+                               std::plus<Float>{});
             }
         }
+
         void computeMassDensity()
         {
             massDensity_.zero();
@@ -103,7 +99,7 @@ namespace core
                 // nodes. This is more efficient and easier to code as we don't
                 // have to account for the field dimensionality.
 
-                auto& popDensity = pop.density();
+                auto& popDensity = pop.particleDensity();
                 std::transform(
                     std::begin(massDensity_), std::end(massDensity_), std::begin(popDensity),
                     std::begin(massDensity_),
@@ -114,14 +110,7 @@ namespace core
 
         void computeBulkVelocity()
         {
-            // the bulk velocity is sum(pop_mass * pop_flux) / sum(pop_mass * pop_density)
-            // if all populations have the same mass, this is equivalent to sum(pop_flux) /
-            // sum(pop_density) sum(pop_density) is rho_ and already known by the time we get here.
-            // sum(pop_mass * pop_flux) is massDensity_ and is computed by computeMassDensity() if
-            // needed
-            if (!sameMasses_)
-                computeMassDensity();
-            auto const& density = (sameMasses_) ? rho_ : massDensity_;
+            computeMassDensity();
 
             bulkVelocity_.zero();
             auto& vx = bulkVelocity_.getComponent(Component::X);
@@ -131,28 +120,26 @@ namespace core
             for (auto& pop : populations_)
             {
                 // account for mass only if populations have different masses
-                std::function<Float(Float, Float)> plus = std::plus<Float>{};
                 std::function<Float(Float, Float)> plusMass
                     = [&pop](Float const& v, Float const& f) { return v + f * pop.mass(); };
 
                 auto const& flux    = pop.flux();
                 auto&& [fx, fy, fz] = flux();
 
-
                 std::transform(std::begin(vx), std::end(vx), std::begin(fx), std::begin(vx),
-                               sameMasses_ ? plus : plusMass);
+                               plusMass);
                 std::transform(std::begin(vy), std::end(vy), std::begin(fy), std::begin(vy),
-                               sameMasses_ ? plus : plusMass);
+                               plusMass);
                 std::transform(std::begin(vz), std::end(vz), std::begin(fz), std::begin(vz),
-                               sameMasses_ ? plus : plusMass);
+                               plusMass);
             }
 
 
-            std::transform(std::begin(vx), std::end(vx), std::begin(density), std::begin(vx),
+            std::transform(std::begin(vx), std::end(vx), std::begin(massDensity_), std::begin(vx),
                            std::divides<Float>{});
-            std::transform(std::begin(vy), std::end(vy), std::begin(density), std::begin(vy),
+            std::transform(std::begin(vy), std::end(vy), std::begin(massDensity_), std::begin(vy),
                            std::divides<Float>{});
-            std::transform(std::begin(vz), std::end(vz), std::begin(density), std::begin(vz),
+            std::transform(std::begin(vz), std::end(vz), std::begin(massDensity_), std::begin(vz),
                            std::divides<Float>{});
         }
 
@@ -185,11 +172,8 @@ namespace core
         // because it is for internal use only so no object will ever need to access it.
         NO_DISCARD bool isUsable() const
         {
-            bool usable
-                = rho_.isUsable() and bulkVelocity_.isUsable() and momentumTensor_.isUsable();
-
-            // if all populations have the same mass, we don't need the massDensity_
-            usable &= (sameMasses_) ? true : massDensity_.isUsable();
+            bool usable = chargeDensity_.isUsable() and bulkVelocity_.isUsable()
+                          and momentumTensor_.isUsable() and massDensity_.isUsable();
 
             for (auto const& pop : populations_)
             {
@@ -202,11 +186,8 @@ namespace core
 
         NO_DISCARD bool isSettable() const
         {
-            bool settable
-                = rho_.isSettable() and bulkVelocity_.isSettable() and momentumTensor_.isSettable();
-
-            // if all populations have the same mass, we don't need the massDensity_
-            settable &= (sameMasses_) ? true : massDensity_.isSettable();
+            bool settable = massDensity_.isSettable() and chargeDensity_.isSettable()
+                            and bulkVelocity_.isSettable() and momentumTensor_.isSettable();
 
             for (auto const& pop : populations_)
             {
@@ -230,7 +211,8 @@ namespace core
 
         NO_DISCARD auto getCompileTimeResourcesViewList()
         {
-            return std::forward_as_tuple(bulkVelocity_, momentumTensor_, rho_, massDensity_);
+            return std::forward_as_tuple(bulkVelocity_, momentumTensor_, chargeDensity_,
+                                         massDensity_);
         }
 
 
@@ -251,29 +233,12 @@ namespace core
             return ss.str();
         }
 
-        NO_DISCARD bool sameMasses() const { return sameMasses_; }
-
 
     private:
-        bool allSameMass_() const
-        {
-            return all(populations_, [this](auto const& pop) { // arbitrary small diff
-                return float_equals(pop.mass(), populations_.front().mass(), /*abs_tol=*/1e-10);
-            });
-        }
-
-
-
-
-        field_type rho_;
         field_type massDensity_;
+        field_type chargeDensity_;
         vecfield_type bulkVelocity_;
         std::vector<IonPopulation> populations_;
-
-        // note this is set at construction when all populations are added
-        // in the future if some populations are dynamically created during the simulation
-        // this should be updated accordingly
-        bool sameMasses_{false};
         tensorfield_type momentumTensor_;
     };
 } // namespace core
diff --git a/src/core/numerics/interpolator/interpolator.hpp b/src/core/numerics/interpolator/interpolator.hpp
index a0a6a5a96..e8858d6a4 100644
--- a/src/core/numerics/interpolator/interpolator.hpp
+++ b/src/core/numerics/interpolator/interpolator.hpp
@@ -469,8 +469,9 @@ class Interpolator : private Weighter<interpOrder>
      * onto the particle.
      */
     template<typename ParticleRange, typename VecField, typename GridLayout, typename Field>
-    inline void operator()(ParticleRange& particleRange, Field& density, VecField& flux,
-                           GridLayout const& layout, double coef = 1.)
+    inline void operator()(ParticleRange& particleRange, Field& particleDensity,
+                           Field& chargeDensity, VecField& flux, GridLayout const& layout,
+                           double coef = 1.)
     {
         auto begin                        = particleRange.begin();
         auto end                          = particleRange.end();
@@ -487,8 +488,11 @@ class Interpolator : private Weighter<interpOrder>
                                                                  currPart->delta);
 
             particleToMesh_(
-                density, *currPart, [](auto const& part) { return 1.; }, startIndex_, weights_,
-                coef);
+                particleDensity, *currPart, [](auto const& part) { return 1.; }, startIndex_,
+                weights_, coef);
+            particleToMesh_(
+                chargeDensity, *currPart, [](auto const& part) { return part.charge; }, startIndex_,
+                weights_, coef);
             particleToMesh_(
                 xFlux, *currPart, [](auto const& part) { return part.v[0]; }, startIndex_, weights_,
                 coef);
@@ -502,10 +506,10 @@ class Interpolator : private Weighter<interpOrder>
         PHARE_LOG_STOP(3, "ParticleToMesh::operator()");
     }
     template<typename ParticleRange, typename VecField, typename GridLayout, typename Field>
-    inline void operator()(ParticleRange&& range, Field& density, VecField& flux,
-                           GridLayout const& layout, double coef = 1.)
+    inline void operator()(ParticleRange&& range, Field& particleDensity, Field& chargeDensity,
+                           VecField& flux, GridLayout const& layout, double coef = 1.)
     {
-        (*this)(range, density, flux, layout, coef);
+        (*this)(range, particleDensity, chargeDensity, flux, layout, coef);
     }
 
 
diff --git a/src/core/numerics/ion_updater/ion_updater.hpp b/src/core/numerics/ion_updater/ion_updater.hpp
index dbae99b75..6f51a6147 100644
--- a/src/core/numerics/ion_updater/ion_updater.hpp
+++ b/src/core/numerics/ion_updater/ion_updater.hpp
@@ -107,7 +107,7 @@ void IonUpdater<Ions, Electromag, GridLayout>::updatePopulations(Ions& ions, Ele
 template<typename Ions, typename Electromag, typename GridLayout>
 void IonUpdater<Ions, Electromag, GridLayout>::updateIons(Ions& ions)
 {
-    ions.computeDensity();
+    ions.computeChargeDensity();
     ions.computeBulkVelocity();
 }
 
@@ -158,7 +158,7 @@ void IonUpdater<Ions, Electromag, GridLayout>::updateAndDepositDomain_(Ions& ion
             inRange, outRange, em, pop.mass(), interpolator_, layout,
             [](auto& particleRange) { return particleRange; }, inDomainBox);
 
-        interpolator_(inDomain, pop.density(), pop.flux(), layout);
+        interpolator_(inDomain, pop.particleDensity(), pop.chargeDensity(), pop.flux(), layout);
 
         // TODO : we can erase here because we know we are working on a state
         // that has been saved in the solverPPC
@@ -181,7 +181,8 @@ void IonUpdater<Ions, Electromag, GridLayout>::updateAndDepositDomain_(Ions& ion
             auto enteredInDomain = pusher_->move(inRange, outRange, em, pop.mass(), interpolator_,
                                                  layout, inGhostBox, inDomainBox);
 
-            interpolator_(enteredInDomain, pop.density(), pop.flux(), layout);
+            interpolator_(enteredInDomain, pop.particleDensity(), pop.chargeDensity(), pop.flux(),
+                          layout);
 
             if (copyInDomain)
             {
@@ -269,7 +270,8 @@ void IonUpdater<Ions, Electromag, GridLayout>::updateAndDepositAll_(Ions& ions,
         pushAndCopyInDomain(makeIndexRange(pop.patchGhostParticles()));
         pushAndCopyInDomain(makeIndexRange(pop.levelGhostParticles()));
 
-        interpolator_(makeIndexRange(domainParticles), pop.density(), pop.flux(), layout);
+        interpolator_(makeIndexRange(domainParticles), pop.particleDensity(), pop.chargeDensity(),
+                      pop.flux(), layout);
     }
 }
 
diff --git a/src/core/numerics/moments/moments.hpp b/src/core/numerics/moments/moments.hpp
index d71135614..87b805fa1 100644
--- a/src/core/numerics/moments/moments.hpp
+++ b/src/core/numerics/moments/moments.hpp
@@ -15,7 +15,8 @@ namespace core
     {
         for (auto& pop : ions)
         {
-            pop.density().zero();
+            pop.particleDensity().zero();
+            pop.chargeDensity().zero();
             pop.flux().zero();
         }
     }
@@ -40,23 +41,24 @@ namespace core
     {
         for (auto& pop : ions)
         {
-            auto& density = pop.density();
-            auto& flux    = pop.flux();
+            auto& particleDensity = pop.particleDensity();
+            auto& chargeDensity   = pop.chargeDensity();
+            auto& flux            = pop.flux();
 
             if constexpr (std::is_same_v<DepositTag, DomainDeposit>)
             {
                 auto& partArray = pop.domainParticles();
-                interpolate(partArray, density, flux, layout);
+                interpolate(partArray, particleDensity, chargeDensity, flux, layout);
             }
             else if constexpr (std::is_same_v<DepositTag, PatchGhostDeposit>)
             {
                 auto& partArray = pop.patchGhostParticles();
-                interpolate(partArray, density, flux, layout);
+                interpolate(partArray, particleDensity, chargeDensity, flux, layout);
             }
             else if constexpr (std::is_same_v<DepositTag, LevelGhostDeposit>)
             {
                 auto& partArray = pop.levelGhostParticlesOld();
-                interpolate(partArray, density, flux, layout);
+                interpolate(partArray, particleDensity, chargeDensity, flux, layout);
             }
         }
     }
diff --git a/src/diagnostic/detail/types/fluid.hpp b/src/diagnostic/detail/types/fluid.hpp
index ce8d98c78..89839c5d7 100644
--- a/src/diagnostic/detail/types/fluid.hpp
+++ b/src/diagnostic/detail/types/fluid.hpp
@@ -146,12 +146,13 @@ void FluidDiagnosticWriter<H5Writer>::createFiles(DiagnosticProperties& diagnost
     for (auto const& pop : this->h5Writer_.modelView().getIons())
     {
         std::string tree{"/ions/pop/" + pop.name() + "/"};
-        checkCreateFileFor_(diagnostic, fileData_, tree, "density", "flux", "momentum_tensor");
+        checkCreateFileFor_(diagnostic, fileData_, tree, "density", "charge_density", "flux",
+                            "momentum_tensor");
     }
 
     std::string tree{"/ions/"};
-    checkCreateFileFor_(diagnostic, fileData_, tree, "density", "mass_density", "bulkVelocity",
-                        "momentum_tensor");
+    checkCreateFileFor_(diagnostic, fileData_, tree, "charge_density", "mass_density",
+                        "bulkVelocity", "momentum_tensor");
 }
 
 
@@ -198,7 +199,9 @@ void FluidDiagnosticWriter<H5Writer>::getDataSetInfo(DiagnosticProperties& diagn
         std::string tree{"/ions/pop/" + pop.name() + "/"};
         auto& popAttr = patchAttributes[lvlPatchID]["fluid_" + pop.name()];
         if (isActiveDiag(diagnostic, tree, "density"))
-            infoDS(pop.density(), "density", popAttr);
+            infoDS(pop.particleDensity(), "density", popAttr);
+        if (isActiveDiag(diagnostic, tree, "charge_density"))
+            infoDS(pop.chargeDensity(), "charge_density", popAttr);
         if (isActiveDiag(diagnostic, tree, "flux"))
             infoVF(pop.flux(), "flux", popAttr);
         if (isActiveDiag(diagnostic, tree, "momentum_tensor"))
@@ -206,8 +209,8 @@ void FluidDiagnosticWriter<H5Writer>::getDataSetInfo(DiagnosticProperties& diagn
     }
 
     std::string tree{"/ions/"};
-    if (isActiveDiag(diagnostic, tree, "density"))
-        infoDS(ions.density(), "density", patchAttributes[lvlPatchID]["ion"]);
+    if (isActiveDiag(diagnostic, tree, "charge_density"))
+        infoDS(ions.chargeDensity(), "charge_density", patchAttributes[lvlPatchID]["ion"]);
     if (isActiveDiag(diagnostic, tree, "mass_density"))
         infoDS(ions.massDensity(), "mass_density", patchAttributes[lvlPatchID]["ion"]);
     if (isActiveDiag(diagnostic, tree, "bulkVelocity"))
@@ -268,6 +271,8 @@ void FluidDiagnosticWriter<H5Writer>::initDataSets(
             std::string popPath(path + "pop/" + pop.name() + "/");
             if (isActiveDiag(diagnostic, tree, "density"))
                 initDS(path, attr[popId], "density", null);
+            if (isActiveDiag(diagnostic, tree, "charge_density"))
+                initDS(path, attr[popId], "charge_density", null);
             if (isActiveDiag(diagnostic, tree, "flux"))
                 initVF(path, attr[popId], "flux", null);
             if (isActiveDiag(diagnostic, tree, "momentum_tensor"))
@@ -275,8 +280,8 @@ void FluidDiagnosticWriter<H5Writer>::initDataSets(
         }
 
         std::string tree{"/ions/"};
-        if (isActiveDiag(diagnostic, tree, "density"))
-            initDS(path, attr["ion"], "density", null);
+        if (isActiveDiag(diagnostic, tree, "charge_density"))
+            initDS(path, attr["ion"], "charge_density", null);
         if (isActiveDiag(diagnostic, tree, "mass_density"))
             initDS(path, attr["ion"], "mass_density", null);
         if (isActiveDiag(diagnostic, tree, "bulkVelocity"))
@@ -307,7 +312,9 @@ void FluidDiagnosticWriter<H5Writer>::write(DiagnosticProperties& diagnostic)
     {
         std::string tree{"/ions/pop/" + pop.name() + "/"};
         if (isActiveDiag(diagnostic, tree, "density"))
-            writeDS(path + "density", pop.density());
+            writeDS(path + "density", pop.particleDensity());
+        if (isActiveDiag(diagnostic, tree, "charge_density"))
+            writeDS(path + "charge_density", pop.chargeDensity());
         if (isActiveDiag(diagnostic, tree, "flux"))
             writeTF(path + "flux", pop.flux());
         if (isActiveDiag(diagnostic, tree, "momentum_tensor"))
@@ -315,8 +322,8 @@ void FluidDiagnosticWriter<H5Writer>::write(DiagnosticProperties& diagnostic)
     }
 
     std::string tree{"/ions/"};
-    if (isActiveDiag(diagnostic, tree, "density"))
-        writeDS(path + "density", ions.density());
+    if (isActiveDiag(diagnostic, tree, "charge_density"))
+        writeDS(path + "charge_density", ions.chargeDensity());
     if (isActiveDiag(diagnostic, tree, "mass_density"))
         writeDS(path + "mass_density", ions.massDensity());
     if (isActiveDiag(diagnostic, tree, "bulkVelocity"))
@@ -345,6 +352,7 @@ void FluidDiagnosticWriter<H5Writer>::writeAttributes(
     {
         std::string tree = "/ions/pop/" + pop.name() + "/";
         checkWrite(tree, "density", pop);
+        checkWrite(tree, "charge_density", pop);
         checkWrite(tree, "flux", pop);
         checkWrite(tree, "momentum_tensor", pop);
     }
diff --git a/src/phare/phare_init.py b/src/phare/phare_init.py
index 2c90bce3a..d79de42b7 100644
--- a/src/phare/phare_init.py
+++ b/src/phare/phare_init.py
@@ -101,7 +101,7 @@ def vthz(x):
     ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
 
-for quantity in ["density", "bulkVelocity"]:
+for quantity in ["charge_density", "bulkVelocity"]:
     ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
 pops = [
diff --git a/src/phare/phare_init_small.py b/src/phare/phare_init_small.py
index b4bccc1cb..f3f8e76e1 100644
--- a/src/phare/phare_init_small.py
+++ b/src/phare/phare_init_small.py
@@ -104,7 +104,7 @@ def vthz(x):
     ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
 
-for quantity in ["density", "bulkVelocity"]:
+for quantity in ["charge_density", "bulkVelocity"]:
     ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
 pops = [
diff --git a/src/python3/patch_level.hpp b/src/python3/patch_level.hpp
index 517c03baa..b7600a44a 100644
--- a/src/python3/patch_level.hpp
+++ b/src/python3/patch_level.hpp
@@ -37,7 +37,7 @@ class PatchLevel
         auto& ions = model_.state.ions;
 
         auto visit = [&](GridLayout& grid, std::string patchID, std::size_t /*iLevel*/) {
-            setPatchDataFromField(patchDatas.emplace_back(), ions.density(), grid, patchID);
+            setPatchDataFromField(patchDatas.emplace_back(), ions.chargeDensity(), grid, patchID);
         };
 
         PHARE::amr::visitLevel<GridLayout>(*hierarchy_.getPatchLevel(lvl_),
@@ -59,7 +59,7 @@ class PatchLevel
                 if (!pop_data.count(pop.name()))
                     pop_data.emplace(pop.name(), Inner());
 
-                setPatchDataFromField(pop_data.at(pop.name()).emplace_back(), pop.density(), grid,
+                setPatchDataFromField(pop_data.at(pop.name()).emplace_back(), pop.chargeDensity(), grid,
                                       patchID);
             }
         };
diff --git a/tests/core/data/electrons/test_electrons.cpp b/tests/core/data/electrons/test_electrons.cpp
index d69df65ff..5f5ffc700 100644
--- a/tests/core/data/electrons/test_electrons.cpp
+++ b/tests/core/data/electrons/test_electrons.cpp
@@ -126,9 +126,8 @@ class NDlayout
 template<typename TypeInfo /*= std::pair<DimConst<1>, InterpConst<1>>*/>
 struct ElectronsTest : public ::testing::Test
 {
-    static constexpr auto dim         = typename TypeInfo::first_type{}();
-    static constexpr auto interp      = typename TypeInfo::second_type{}();
-    static constexpr auto densityName = std::string_view{"rho"};
+    static constexpr auto dim    = typename TypeInfo::first_type{}();
+    static constexpr auto interp = typename TypeInfo::second_type{}();
 
 
     using GridYee = GridLayout<GridLayoutImplYee<dim, interp>>;
@@ -149,9 +148,9 @@ struct ElectronsTest : public ::testing::Test
     Electromag<VecFieldND> electromag;
 
     UsableVecField<dim> J, F, Ve, Vi;
-    UsableTensorField<dim> M, protons_M;
+    UsableTensorField<dim> ionTensor, protonTensor;
 
-    GridND Nibuffer, NiProtons, Pe;
+    GridND ionChargeDensity, ionMassDensity, protonParticleDensity, protonChargeDensity, Pe;
 
     ParticleArray_t domainParticles{layout.AMRBox()};
     ParticleArray_t patchGhostParticles = domainParticles;
@@ -164,22 +163,25 @@ struct ElectronsTest : public ::testing::Test
     template<typename... Args>
     auto static _ions(Args&... args)
     {
-        auto const& [Fi, Nibuffer, NiProtons, Vi, M, protons_M, pack]
+        auto const& [ionFlux, ionChargeDensity, ionMassDensity, protonParticleDensity,
+                     protonChargeDensity, Vi, ionTensor, protonTensor, pack]
             = std::forward_as_tuple(args...);
         IonsT ions{createDict<dim>()["ions"]};
         {
-            auto const& [V, m, d, md] = ions.getCompileTimeResourcesViewList();
-            d.setBuffer(&Nibuffer);
+            auto const& [V, m, d_c, d_m] = ions.getCompileTimeResourcesViewList();
+            d_c.setBuffer(&ionChargeDensity);
+            d_m.setBuffer(&ionMassDensity);
             Vi.set_on(V);
-            M.set_on(m);
+            ionTensor.set_on(m);
         }
         auto& pops = ions.getRunTimeResourcesViewList();
         assert(pops.size() == 1);
 
-        auto const& [F, m, d, poppack] = pops[0].getCompileTimeResourcesViewList();
-        d.setBuffer(&NiProtons);
-        Fi.set_on(F);
-        protons_M.set_on(m);
+        auto const& [F, m, Np, Nc, poppack] = pops[0].getCompileTimeResourcesViewList();
+        Np.setBuffer(&protonParticleDensity);
+        Nc.setBuffer(&protonChargeDensity);
+        ionFlux.set_on(F);
+        protonTensor.set_on(m);
         poppack.setBuffer(&pack);
         return ions;
     }
@@ -191,14 +193,19 @@ struct ElectronsTest : public ::testing::Test
         , F{"protons_flux", layout, HybridQuantity::Vector::V}
         , Ve{"StandardHybridElectronFluxComputer_Ve", layout, HybridQuantity::Vector::V}
         , Vi{"bulkVel", layout, HybridQuantity::Vector::V}
-        , M{"momentumTensor", layout, HybridQuantity::Tensor::M}
-        , protons_M{"protons_momentumTensor", layout, HybridQuantity::Tensor::M}
-        , Nibuffer{std::string{densityName}, HybridQuantity::Scalar::rho,
-                   layout.allocSize(HybridQuantity::Scalar::rho)}
-        , NiProtons{"protons_rho", HybridQuantity::Scalar::rho,
-                    layout.allocSize(HybridQuantity::Scalar::rho)}
+        , ionTensor{"momentumTensor", layout, HybridQuantity::Tensor::M}
+        , protonTensor{"protons_momentumTensor", layout, HybridQuantity::Tensor::M}
+        , ionChargeDensity{"chargeDensity", HybridQuantity::Scalar::rho,
+                           layout.allocSize(HybridQuantity::Scalar::rho)}
+        , ionMassDensity{"massDensity", HybridQuantity::Scalar::rho,
+                         layout.allocSize(HybridQuantity::Scalar::rho)}
+        , protonParticleDensity{"protons_particleDensity", HybridQuantity::Scalar::rho,
+                                layout.allocSize(HybridQuantity::Scalar::rho)}
+        , protonChargeDensity{"protons_chargeDensity", HybridQuantity::Scalar::rho,
+                              layout.allocSize(HybridQuantity::Scalar::rho)}
         , Pe{"Pe", HybridQuantity::Scalar::P, layout.allocSize(HybridQuantity::Scalar::P)}
-        , ions{_ions(F, Nibuffer, NiProtons, Vi, M, protons_M, pack)}
+        , ions{_ions(F, ionChargeDensity, ionMassDensity, protonParticleDensity,
+                     protonChargeDensity, Vi, ionTensor, protonTensor, pack)}
         , electrons{createDict<dim>()["electrons"], ions, J}
     {
         auto&& emm = std::get<0>(electrons.getCompileTimeResourcesViewList());
@@ -237,7 +244,7 @@ struct ElectronsTest : public ::testing::Test
             fill(Jy, [](double x) { return std::sinh(0.3 * x); });
             fill(Jz, [](double x) { return std::sinh(0.4 * x); });
 
-            fill(Nibuffer, [](double x) { return std::cosh(0.1 * x); });
+            fill(ionChargeDensity, [](double x) { return std::cosh(0.1 * x); });
         }
         else if constexpr (dim == 2)
         {
@@ -266,7 +273,7 @@ struct ElectronsTest : public ::testing::Test
             fill(Jy, [](double x, double y) { return std::sinh(0.3 * x) * std::sinh(0.3 * y); });
             fill(Jz, [](double x, double y) { return std::sinh(0.4 * x) * std::sinh(0.4 * y); });
 
-            fill(Nibuffer,
+            fill(ionChargeDensity,
                  [](double x, double y) { return std::cosh(0.1 * x) * std::cosh(0.1 * y); });
         }
         else if constexpr (dim == 3)
@@ -313,7 +320,7 @@ struct ElectronsTest : public ::testing::Test
                 return std::sinh(0.4 * x) * std::sinh(0.4 * y) * std::sinh(0.4 * z);
             });
 
-            fill(Nibuffer, [](double x, double y, double z) {
+            fill(ionChargeDensity, [](double x, double y, double z) {
                 return std::cosh(0.1 * x) * std::cosh(0.1 * y) * std::cosh(0.1 * z);
             });
         }
@@ -363,7 +370,7 @@ TYPED_TEST(ElectronsTest, ThatElectronsDensityEqualIonDensity)
     electrons.update(layout);
 
     auto& Ne = electrons.density();
-    auto& Ni = ions.density();
+    auto& Ni = ions.chargeDensity();
 
     if constexpr (dim == 1)
     {
diff --git a/tests/core/data/ion_population/test_ion_population_fixtures.hpp b/tests/core/data/ion_population/test_ion_population_fixtures.hpp
index e26ee34e3..1fe25487d 100644
--- a/tests/core/data/ion_population/test_ion_population_fixtures.hpp
+++ b/tests/core/data/ion_population/test_ion_population_fixtures.hpp
@@ -50,7 +50,7 @@ class UsableIonsPopulation : public _defaults::IonPopulation_t
 public:
     UsableIonsPopulation(initializer::PHAREDict const& dict, GridLayout_t const& layout)
         : Super{dict}
-        , rho{this->name() + "_rho", layout, HybridQuantity::Scalar::rho}
+        , rho{this->name() + "_particleDensity", layout, HybridQuantity::Scalar::rho}
         , F{this->name() + "_flux", layout, HybridQuantity::Vector::V}
         , M{this->name() + "_momentumTensor", layout, HybridQuantity::Tensor::M}
         , particles{this->name(), layout.AMRBox()}
diff --git a/tests/core/data/ions/test_ions.cpp b/tests/core/data/ions/test_ions.cpp
index 7fb1b1426..11769fac2 100644
--- a/tests/core/data/ions/test_ions.cpp
+++ b/tests/core/data/ions/test_ions.cpp
@@ -154,7 +154,7 @@ TEST_F(theIons, areSettableUponConstruction)
 #ifndef NDEBUG // no throw in release mode! JUST SEGFAULTS! :D
 TEST_F(theIons, throwIfAccessingDensityWhileNotUsable)
 {
-    EXPECT_ANY_THROW(auto& n = ions.density()(0));
+    EXPECT_ANY_THROW(auto& n = ions.chargeDensity()(0));
 }
 #endif
 
diff --git a/tests/core/numerics/interpolator/test_main.cpp b/tests/core/numerics/interpolator/test_main.cpp
index 28e2fff14..9d1577d39 100644
--- a/tests/core/numerics/interpolator/test_main.cpp
+++ b/tests/core/numerics/interpolator/test_main.cpp
@@ -516,6 +516,7 @@ class ACollectionOfParticles_1d : public ::testing::Test
     ParticleArray_t particles;
 
     Grid_t rho;
+    Grid_t rho_c;
 
     UsableVecFieldND v;
     std::array<double, nbrPointsSupport(Interpolator::interp_order)> weights;
@@ -526,6 +527,7 @@ class ACollectionOfParticles_1d : public ::testing::Test
         : part{}
         , particles{grow(layout.AMRBox(), safeLayer)}
         , rho{"field", HybridQuantity::Scalar::rho, nx}
+        , rho_c{"field", HybridQuantity::Scalar::rho, nx}
         , v{"v", layout, HybridQuantity::Vector::V}
     {
         if constexpr (Interpolator::interp_order == 1)
@@ -632,7 +634,7 @@ class ACollectionOfParticles_1d : public ::testing::Test
             part.v[2]     = +1.;
             particles.push_back(part);
         }
-        interpolator(makeIndexRange(particles), rho, v, layout);
+        interpolator(makeIndexRange(particles), rho, rho_c, v, layout);
     }
 
 
@@ -678,12 +680,14 @@ struct ACollectionOfParticles_2d : public ::testing::Test
     GridLayout_t layout{ConstArray<double, dim>(.1), {nx, ny}, ConstArray<double, dim>(0)};
     ParticleArray_t particles;
     Grid_t rho;
+    Grid_t rho_c;
     UsableVecFieldND v;
     Interpolator interpolator;
 
     ACollectionOfParticles_2d()
         : particles{grow(layout.AMRBox(), safeLayer)}
         , rho{"field", HybridQuantity::Scalar::rho, nx, ny}
+        , rho_c{"field", HybridQuantity::Scalar::rho, nx, ny}
         , v{"v", layout, HybridQuantity::Vector::V}
     {
         for (int i = start; i < end; i++)
@@ -697,7 +701,7 @@ struct ACollectionOfParticles_2d : public ::testing::Test
                 part.v[1]   = -1.;
                 part.v[2]   = +1.;
             }
-        interpolator(makeIndexRange(particles), rho, v, layout);
+        interpolator(makeIndexRange(particles), rho, rho_c, v, layout);
     }
 };
 TYPED_TEST_SUITE_P(ACollectionOfParticles_2d);
diff --git a/tests/core/numerics/ion_updater/test_updater.cpp b/tests/core/numerics/ion_updater/test_updater.cpp
index 2d939cf7d..2ab71bfb4 100644
--- a/tests/core/numerics/ion_updater/test_updater.cpp
+++ b/tests/core/numerics/ion_updater/test_updater.cpp
@@ -211,10 +211,12 @@ struct IonsBuffers
     using ParticleArray              = typename PHARETypes::ParticleArray_t;
     using ParticleInitializerFactory = typename PHARETypes::ParticleInitializerFactory;
 
-    Grid ionDensity;
+    Grid ionChargeDensity;
     Grid ionMassDensity;
-    Grid protonDensity;
-    Grid alphaDensity;
+    Grid protonParticleDensity;
+    Grid protonChargeDensity;
+    Grid alphaParticleDensity;
+    Grid alphaChargeDensity;
 
     UsableVecFieldND protonF, alphaF, Vi;
     UsableTensorField<dim> M, alpha_M, protons_M;
@@ -237,14 +239,18 @@ struct IonsBuffers
     ParticlesPack<ParticleArray> alphaPack;
 
     IonsBuffers(GridLayout const& layout)
-        : ionDensity{"rho", HybridQuantity::Scalar::rho,
-                     layout.allocSize(HybridQuantity::Scalar::rho)}
+        : ionChargeDensity{"chargeDensity", HybridQuantity::Scalar::rho,
+                           layout.allocSize(HybridQuantity::Scalar::rho)}
         , ionMassDensity{"massDensity", HybridQuantity::Scalar::rho,
                          layout.allocSize(HybridQuantity::Scalar::rho)}
-        , protonDensity{"protons_rho", HybridQuantity::Scalar::rho,
-                        layout.allocSize(HybridQuantity::Scalar::rho)}
-        , alphaDensity{"alpha_rho", HybridQuantity::Scalar::rho,
-                       layout.allocSize(HybridQuantity::Scalar::rho)}
+        , protonParticleDensity{"protons_particleDensity", HybridQuantity::Scalar::rho,
+                                layout.allocSize(HybridQuantity::Scalar::rho)}
+        , protonChargeDensity{"protons_chargeDensity", HybridQuantity::Scalar::rho,
+                              layout.allocSize(HybridQuantity::Scalar::rho)}
+        , alphaParticleDensity{"alpha_particleDensity", HybridQuantity::Scalar::rho,
+                               layout.allocSize(HybridQuantity::Scalar::rho)}
+        , alphaChargeDensity{"alpha_chargeDensity", HybridQuantity::Scalar::rho,
+                             layout.allocSize(HybridQuantity::Scalar::rho)}
         , protonF{"protons_flux", layout, HybridQuantity::Vector::V}
         , alphaF{"alpha_flux", layout, HybridQuantity::Vector::V}
         , Vi{"bulkVel", layout, HybridQuantity::Vector::V}
@@ -270,14 +276,18 @@ struct IonsBuffers
 
 
     IonsBuffers(IonsBuffers const& source, GridLayout const& layout)
-        : ionDensity{"rho", HybridQuantity::Scalar::rho,
-                     layout.allocSize(HybridQuantity::Scalar::rho)}
+        : ionChargeDensity{"chargeDensity", HybridQuantity::Scalar::rho,
+                           layout.allocSize(HybridQuantity::Scalar::rho)}
         , ionMassDensity{"massDensity", HybridQuantity::Scalar::rho,
                          layout.allocSize(HybridQuantity::Scalar::rho)}
-        , protonDensity{"protons_rho", HybridQuantity::Scalar::rho,
-                        layout.allocSize(HybridQuantity::Scalar::rho)}
-        , alphaDensity{"alpha_rho", HybridQuantity::Scalar::rho,
-                       layout.allocSize(HybridQuantity::Scalar::rho)}
+        , protonParticleDensity{"protons_particleDensity", HybridQuantity::Scalar::rho,
+                                layout.allocSize(HybridQuantity::Scalar::rho)}
+        , protonChargeDensity{"protons_chargeDensity", HybridQuantity::Scalar::rho,
+                              layout.allocSize(HybridQuantity::Scalar::rho)}
+        , alphaParticleDensity{"alpha_particleDensity", HybridQuantity::Scalar::rho,
+                               layout.allocSize(HybridQuantity::Scalar::rho)}
+        , alphaChargeDensity{"alpha_chargeDensity", HybridQuantity::Scalar::rho,
+                             layout.allocSize(HybridQuantity::Scalar::rho)}
         , protonF{"protons_flux", layout, HybridQuantity::Vector::V}
         , alphaF{"alpha_flux", layout, HybridQuantity::Vector::V}
         , Vi{"bulkVel", layout, HybridQuantity::Vector::V}
@@ -300,10 +310,12 @@ struct IonsBuffers
                     &alphaLevelGhost, &alphaLevelGhostOld, &alphaLevelGhostNew}
 
     {
-        ionDensity.copyData(source.ionDensity);
+        ionChargeDensity.copyData(source.ionChargeDensity);
         ionMassDensity.copyData(source.ionMassDensity);
-        protonDensity.copyData(source.protonDensity);
-        alphaDensity.copyData(source.alphaDensity);
+        protonParticleDensity.copyData(source.protonParticleDensity);
+        protonChargeDensity.copyData(source.protonChargeDensity);
+        alphaParticleDensity.copyData(source.alphaParticleDensity);
+        alphaChargeDensity.copyData(source.alphaChargeDensity);
 
         protonF.copyData(source.protonF);
         alphaF.copyData(source.alphaF);
@@ -313,25 +325,27 @@ struct IonsBuffers
     void setBuffers(Ions& ions)
     {
         {
-            auto const& [V, m, d, md] = ions.getCompileTimeResourcesViewList();
+            auto const& [V, m, cd, md] = ions.getCompileTimeResourcesViewList();
             Vi.set_on(V);
             M.set_on(m);
-            d.setBuffer(&ionDensity);
+            cd.setBuffer(&ionChargeDensity);
             md.setBuffer(&ionMassDensity);
         }
 
         auto& pops = ions.getRunTimeResourcesViewList();
         {
-            auto const& [F, M, d, particles] = pops[0].getCompileTimeResourcesViewList();
-            d.setBuffer(&protonDensity);
+            auto const& [F, M, d, c, particles] = pops[0].getCompileTimeResourcesViewList();
+            d.setBuffer(&protonParticleDensity);
+            c.setBuffer(&protonChargeDensity);
             protons_M.set_on(M);
             protonF.set_on(F);
             particles.setBuffer(&protonPack);
         }
 
         {
-            auto const& [F, M, d, particles] = pops[1].getCompileTimeResourcesViewList();
-            d.setBuffer(&alphaDensity);
+            auto const& [F, M, d, c, particles] = pops[1].getCompileTimeResourcesViewList();
+            d.setBuffer(&alphaParticleDensity);
+            c.setBuffer(&alphaChargeDensity);
             alpha_M.set_on(M);
             alphaF.set_on(F);
             particles.setBuffer(&alphaPack);
@@ -515,7 +529,7 @@ struct IonUpdaterTest : public ::testing::Test
                 }
 
             } // end 1D
-        }     // end pop loop
+        } // end pop loop
         PHARE::core::depositParticles(ions, layout, Interpolator<dim, interp_order>{},
                                       PHARE::core::DomainDeposit{});
 
@@ -526,7 +540,7 @@ struct IonUpdaterTest : public ::testing::Test
                                       PHARE::core::LevelGhostDeposit{});
 
 
-        ions.computeDensity();
+        ions.computeChargeDensity();
         ions.computeBulkVelocity();
     } // end Ctor
 
@@ -539,17 +553,17 @@ struct IonUpdaterTest : public ::testing::Test
 
         for (auto& pop : this->ions)
         {
-            interpolate(makeIndexRange(pop.patchGhostParticles()), pop.density(), pop.flux(),
-                        layout);
+            interpolate(makeIndexRange(pop.patchGhostParticles()), pop.particleDensity(),
+                        pop.chargeDensity(), pop.flux(), layout);
 
             double alpha = 0.5;
-            interpolate(makeIndexRange(pop.levelGhostParticlesNew()), pop.density(), pop.flux(),
-                        layout,
+            interpolate(makeIndexRange(pop.levelGhostParticlesNew()), pop.particleDensity(),
+                        pop.chargeDensity(), pop.flux(), layout,
                         /*coef = */ alpha);
 
 
-            interpolate(makeIndexRange(pop.levelGhostParticlesOld()), pop.density(), pop.flux(),
-                        layout,
+            interpolate(makeIndexRange(pop.levelGhostParticlesOld()), pop.particleDensity(),
+                        pop.chargeDensity(), pop.flux(), layout,
                         /*coef = */ (1. - alpha));
         }
     }
@@ -560,15 +574,17 @@ struct IonUpdaterTest : public ::testing::Test
     {
         auto& populations = this->ions.getRunTimeResourcesViewList();
 
-        auto& protonDensity = populations[0].density();
-        auto& protonFx      = populations[0].flux().getComponent(Component::X);
-        auto& protonFy      = populations[0].flux().getComponent(Component::Y);
-        auto& protonFz      = populations[0].flux().getComponent(Component::Z);
+        auto& protonParticleDensity = populations[0].particleDensity();
+        auto& protonChargeDensity   = populations[0].chargeDensity();
+        auto& protonFx              = populations[0].flux().getComponent(Component::X);
+        auto& protonFy              = populations[0].flux().getComponent(Component::Y);
+        auto& protonFz              = populations[0].flux().getComponent(Component::Z);
 
-        auto& alphaDensity = populations[1].density();
-        auto& alphaFx      = populations[1].flux().getComponent(Component::X);
-        auto& alphaFy      = populations[1].flux().getComponent(Component::Y);
-        auto& alphaFz      = populations[1].flux().getComponent(Component::Z);
+        auto& alphaParticleDensity = populations[1].particleDensity();
+        auto& alphaChargeDensity   = populations[1].chargeDensity();
+        auto& alphaFx              = populations[1].flux().getComponent(Component::X);
+        auto& alphaFy              = populations[1].flux().getComponent(Component::Y);
+        auto& alphaFz              = populations[1].flux().getComponent(Component::Z);
 
         auto ix0 = this->layout.physicalStartIndex(QtyCentering::primal, Direction::X);
         auto ix1 = this->layout.physicalEndIndex(QtyCentering::primal, Direction::X);
@@ -597,19 +613,18 @@ struct IonUpdaterTest : public ::testing::Test
             }
         };
 
-        check(protonDensity, ionsBufferCpy.protonDensity);
-
-
+        check(protonParticleDensity, ionsBufferCpy.protonParticleDensity);
+        check(protonChargeDensity, ionsBufferCpy.protonChargeDensity);
         check(protonFx, ionsBufferCpy.protonF(Component::X));
         check(protonFy, ionsBufferCpy.protonF(Component::Y));
         check(protonFz, ionsBufferCpy.protonF(Component::Z));
 
-        check(alphaDensity, ionsBufferCpy.alphaDensity);
+        check(alphaParticleDensity, ionsBufferCpy.alphaParticleDensity);
+        check(alphaChargeDensity, ionsBufferCpy.alphaChargeDensity);
         check(alphaFx, ionsBufferCpy.alphaF(Component::X));
         check(alphaFy, ionsBufferCpy.alphaF(Component::Y));
         check(alphaFz, ionsBufferCpy.alphaF(Component::Z));
 
-        check(ions.density(), ionsBufferCpy.ionDensity);
         check(ions.velocity().getComponent(Component::X), ionsBufferCpy.Vi(Component::X));
         check(ions.velocity().getComponent(Component::Y), ionsBufferCpy.Vi(Component::Y));
         check(ions.velocity().getComponent(Component::Z), ionsBufferCpy.Vi(Component::Z));
@@ -650,12 +665,12 @@ struct IonUpdaterTest : public ::testing::Test
             }
         };
 
-        auto& populations   = this->ions.getRunTimeResourcesViewList();
-        auto& protonDensity = populations[0].density();
-        auto& alphaDensity  = populations[1].density();
+        auto& populations           = this->ions.getRunTimeResourcesViewList();
+        auto& protonParticleDensity = populations[0].particleDensity();
+        auto& alphaParticleDensity  = populations[1].particleDensity();
 
-        check(protonDensity, density);
-        check(alphaDensity, density);
+        check(protonParticleDensity, density);
+        check(alphaParticleDensity, density);
     }
 };
 
@@ -909,7 +924,8 @@ TYPED_TEST(IonUpdaterTest, thatNoNaNsExistOnPhysicalNodesMoments)
     {
         for (auto ix = ix0; ix <= ix1; ++ix)
         {
-            auto& density = pop.density();
+            auto& density = pop.particleDensity();
+            auto& chargeDensity = pop.chargeDensity();
             auto& flux    = pop.flux();
 
             auto& fx = flux.getComponent(Component::X);
@@ -917,6 +933,7 @@ TYPED_TEST(IonUpdaterTest, thatNoNaNsExistOnPhysicalNodesMoments)
             auto& fz = flux.getComponent(Component::Z);
 
             EXPECT_FALSE(std::isnan(density(ix)));
+            EXPECT_FALSE(std::isnan(chargeDensity(ix)));
             EXPECT_FALSE(std::isnan(fx(ix)));
             EXPECT_FALSE(std::isnan(fy(ix)));
             EXPECT_FALSE(std::isnan(fz(ix)));
diff --git a/tests/diagnostic/__init__.py b/tests/diagnostic/__init__.py
index fb7c52d96..370bfe7a3 100644
--- a/tests/diagnostic/__init__.py
+++ b/tests/diagnostic/__init__.py
@@ -28,7 +28,12 @@ def dump_all_diags(pops=[], flush_every=100, timestamps=None):
     #     write_timestamps=timestamps
     # )
 
-    for quantity in ["density", "bulkVelocity", "pressure_tensor"]:
+    for quantity in [
+        "charge_density",
+        "mass_density",
+        "bulkVelocity",
+        "pressure_tensor",
+    ]:
         ph.FluidDiagnostics(
             quantity=quantity,
             write_timestamps=timestamps,
@@ -36,7 +41,7 @@ def dump_all_diags(pops=[], flush_every=100, timestamps=None):
         )
 
     for pop in pops:
-        for quantity in ["density", "flux", "pressure_tensor"]:
+        for quantity in ["density", "charge_density", "flux", "pressure_tensor"]:
             ph.FluidDiagnostics(
                 quantity=quantity,
                 write_timestamps=timestamps,
diff --git a/tests/diagnostic/test_diagnostics.hpp b/tests/diagnostic/test_diagnostics.hpp
index f397148a8..02b24b052 100644
--- a/tests/diagnostic/test_diagnostics.hpp
+++ b/tests/diagnostic/test_diagnostics.hpp
@@ -122,10 +122,10 @@ void validateFluidDump(Simulator& sim, Hi5Diagnostic& hi5)
         auto& ions = hi5.modelView.getIons();
         for (auto& pop : ions)
         {
-            checkF(layout, path, "/ions/pop/" + pop.name(), "/density"s, pop.density());
+            checkF(layout, path, "/ions/pop/" + pop.name(), "/density"s, pop.chargeDensity());
             checkVF(layout, path, "/ions/pop/" + pop.name(), "/flux"s, pop.flux());
         }
-        checkF(layout, path, "/ions"s, "/density"s, ions.density());
+        checkF(layout, path, "/ions"s, "/charge_density"s, ions.chargeDensity());
 
         std::string tree{"/ions"}, var{"/bulkVelocity"};
         auto hifile = hi5.writer.makeFile(hi5.writer.fileString(tree + var), hi5.flags_);
@@ -226,7 +226,7 @@ void validateAttributes(Simulator& sim, Hi5Diagnostic& hi5)
     using GridLayout                           = typename Simulator::PHARETypes::GridLayout_t;
     constexpr auto dimension                   = Simulator::dimension;
     constexpr std::size_t expectedPopNbr       = 2;
-    constexpr std::size_t expectedPopAttrFiles = 5;
+    constexpr std::size_t expectedPopAttrFiles = 6;
 
     std::string const ionsPopPath = "/ions/pop/";
 
@@ -236,7 +236,8 @@ void validateAttributes(Simulator& sim, Hi5Diagnostic& hi5)
     auto nbrPop = dict["simulation"]["ions"]["nbrPopulations"].template to<std::size_t>();
     EXPECT_EQ(nbrPop, expectedPopNbr);
 
-    std::vector<std::string> h5FileTypes{"/EM_B", "/EM_E", "/ions/density", "/ions/bulkVelocity"};
+    std::vector<std::string> h5FileTypes{"/EM_B", "/EM_E", "/ions/charge_density",
+                                         "/ions/mass_density", "/ions/bulkVelocity"};
 
     for (std::size_t i = 0; i < nbrPop; ++i)
     {
@@ -247,6 +248,7 @@ void validateAttributes(Simulator& sim, Hi5Diagnostic& hi5)
         h5FileTypes.emplace_back(ionsPopPath + popName + "/levelGhost");
         h5FileTypes.emplace_back(ionsPopPath + popName + "/patchGhost");
         h5FileTypes.emplace_back(ionsPopPath + popName + "/density");
+        h5FileTypes.emplace_back(ionsPopPath + popName + "/charge_density");
         h5FileTypes.emplace_back(ionsPopPath + popName + "/flux");
     }
 
diff --git a/tests/diagnostic/test_diagnostics.ipp b/tests/diagnostic/test_diagnostics.ipp
index b6813b165..e108c465b 100644
--- a/tests/diagnostic/test_diagnostics.ipp
+++ b/tests/diagnostic/test_diagnostics.ipp
@@ -14,7 +14,7 @@ void fluid_test(Simulator&& sim, std::string out_dir)
 
     { // Scoped to destruct after dump
         Hi5Diagnostic<Hierarchy, HybridModel> hi5{hierarchy, hybridModel, out_dir, NEW_HI5_FILE};
-        hi5.dMan.addDiagDict(hi5.fluid("/ions/density"))
+        hi5.dMan.addDiagDict(hi5.fluid("/ions/charge_density"))
             .addDiagDict(hi5.fluid("/ions/bulkVelocity"))
             .addDiagDict(hi5.fluid("/ions/momentum_tensor"))
             .addDiagDict(hi5.fluid("/ions/pop/alpha/density"))
diff --git a/tests/functional/alfven_wave/alfven_wave1d.py b/tests/functional/alfven_wave/alfven_wave1d.py
index dc830ce13..1fe582bc1 100644
--- a/tests/functional/alfven_wave/alfven_wave1d.py
+++ b/tests/functional/alfven_wave/alfven_wave1d.py
@@ -96,7 +96,7 @@ def vthz(x):
     for quantity in ["E", "B"]:
         ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
-    for quantity in ["density", "bulkVelocity"]:
+    for quantity in ["charge_density", "bulkVelocity"]:
         ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
     return sim
diff --git a/tests/functional/dispersion/dispersion.py b/tests/functional/dispersion/dispersion.py
index ad6619984..67ab8d831 100644
--- a/tests/functional/dispersion/dispersion.py
+++ b/tests/functional/dispersion/dispersion.py
@@ -171,7 +171,7 @@ def vthz(x):
             quantity=quantity, write_timestamps=timestamps, flush_every=0
         )
 
-    for quantity in ["density", "bulkVelocity"]:
+    for quantity in ["charge_density", "bulkVelocity"]:
         ph.FluidDiagnostics(
             quantity=quantity, write_timestamps=timestamps, flush_every=0
         )
diff --git a/tests/functional/harris/harris_2d_lb.py b/tests/functional/harris/harris_2d_lb.py
index 1153b6e8b..5b553ffb6 100644
--- a/tests/functional/harris/harris_2d_lb.py
+++ b/tests/functional/harris/harris_2d_lb.py
@@ -134,11 +134,11 @@ def vthxyz(x, y):
 
     for quantity in ["E", "B"]:
         ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps)
-    for quantity in ["density", "bulkVelocity"]:
+    for quantity in ["charge_density", "bulkVelocity"]:
         ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
     ph.FluidDiagnostics(
-        quantity="density", write_timestamps=timestamps, population_name="protons"
+        quantity="charge_density", write_timestamps=timestamps, population_name="protons"
     )
     ph.InfoDiagnostics(quantity="particle_count")
 
diff --git a/tests/functional/shock/shock.py b/tests/functional/shock/shock.py
index 7fd0a6ec2..6cc63f297 100644
--- a/tests/functional/shock/shock.py
+++ b/tests/functional/shock/shock.py
@@ -107,7 +107,7 @@ def vthz(x):
     for quantity in ["E", "B"]:
         ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
-    for quantity in ["density", "bulkVelocity"]:
+    for quantity in ["charge_density", "bulkVelocity"]:
         ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
     return sim
diff --git a/tests/functional/td/td1d.py b/tests/functional/td/td1d.py
index 7a715a182..e21d9f28d 100644
--- a/tests/functional/td/td1d.py
+++ b/tests/functional/td/td1d.py
@@ -100,7 +100,7 @@ def vthz(x):
     for quantity in ["E", "B"]:
         ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
-    for quantity in ["density", "bulkVelocity"]:
+    for quantity in ["charge_density", "bulkVelocity"]:
         ph.FluidDiagnostics(
             quantity=quantity,
             write_timestamps=timestamps,
diff --git a/tests/functional/tdtagged/td1dtagged.py b/tests/functional/tdtagged/td1dtagged.py
index 2d30e59f5..35ab9ec36 100644
--- a/tests/functional/tdtagged/td1dtagged.py
+++ b/tests/functional/tdtagged/td1dtagged.py
@@ -120,7 +120,7 @@ def config(**options):
 
     for quantity in ["E", "B"]:
         ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps)
-    for quantity in ["density", "bulkVelocity"]:
+    for quantity in ["charge_density", "bulkVelocity"]:
         ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
     for pop in sim.model.populations:
diff --git a/tests/functional/translation/translat1d.py b/tests/functional/translation/translat1d.py
index 53eb7daa0..6a5321698 100644
--- a/tests/functional/translation/translat1d.py
+++ b/tests/functional/translation/translat1d.py
@@ -87,7 +87,7 @@ def vthz(x):
     for quantity in ["E", "B"]:
         ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
-    for quantity in ["density", "bulkVelocity"]:
+    for quantity in ["charge_density", "bulkVelocity"]:
         ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
     return sim
@@ -182,7 +182,7 @@ def vthz(x):
     for quantity in ["E", "B"]:
         ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
-    for quantity in ["density", "bulkVelocity"]:
+    for quantity in ["charge_density", "bulkVelocity"]:
         ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
     return sim
diff --git a/tests/simulator/initialize/CMakeLists.txt b/tests/simulator/initialize/CMakeLists.txt
index cbaf44d02..496c45ab4 100644
--- a/tests/simulator/initialize/CMakeLists.txt
+++ b/tests/simulator/initialize/CMakeLists.txt
@@ -20,4 +20,6 @@ if(HighFive)
 
   endif()
 
+  phare_mpi_python3_exec(9 ${PHARE_MPI_PROCS} init-densities test_density_init.py  ${CMAKE_CURRENT_BINARY_DIR})
+
 endif()
diff --git a/tests/simulator/initialize/test_density_init.py b/tests/simulator/initialize/test_density_init.py
new file mode 100644
index 000000000..9ed30faae
--- /dev/null
+++ b/tests/simulator/initialize/test_density_init.py
@@ -0,0 +1,395 @@
+import os
+
+from pyphare.simulator.simulator import Simulator
+from pyphare.pharesee.run import Run
+
+import pyphare.pharein as ph
+
+from pyphare.pharein import global_vars
+from tests.diagnostic import all_timestamps
+
+from pyphare.pharein.diagnostics import FluidDiagnostics
+
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+import numpy as np
+
+from pyphare.pharesee.hierarchy import hierarchy_from
+
+mpl.use("Agg")
+
+
+def ions_mass_density_func1d(x, **kwargs):
+    masses = kwargs["masses"]  # list of float : the ion pop masses
+    densities = kwargs["densities"]  # list of callable : the ion pop density profiles
+    if len(masses) != len(densities):
+        raise ValueError("Length of masses and densities must be equal")
+
+    funcs = np.zeros((x.size, len(masses)))
+
+    for i, (mass, density) in enumerate(zip(masses, densities)):
+        funcs[:, i] = mass * density(x)
+
+    return funcs.sum(axis=1)
+
+
+def ions_charge_density_func1d(x, **kwargs):
+    charges = kwargs["charges"]  # list of float : the ion pop charges
+    densities = kwargs["densities"]  # list of callable : the ion pop density profiles
+    if len(charges) != len(densities):
+        raise ValueError("Length of charges and densities must be equal")
+
+    funcs = np.zeros((x.size, len(charges)))
+
+    for i, (charge, density) in enumerate(zip(charges, densities)):
+        funcs[:, i] = charge * density(x)
+
+    return funcs.sum(axis=1)
+
+
+def ions_mass_density_func2d(x, y, **kwargs):
+    masses = kwargs["masses"]  # list of float : the ion pop masses
+    densities = kwargs["densities"]  # list of callable : the ion pop density profiles
+    if len(masses) != len(densities):
+        raise ValueError("Length of masses and densities must be equal")
+
+    yv, xv = np.meshgrid(y, x)
+
+    funcs = np.zeros((x.size, y.size, len(masses)))
+
+    for i, (mass, density) in enumerate(zip(masses, densities)):
+        funcs[:, :, i] = mass * density(xv, yv)
+
+    return funcs.sum(axis=2)
+
+
+def ions_charge_density_func2d(x, y, **kwargs):
+    charges = kwargs["charges"]  # list of float : the ion pop charges
+    densities = kwargs["densities"]  # list of callable : the ion pop density profiles
+    if len(charges) != len(densities):
+        raise ValueError("Length of charges and densities must be equal")
+
+    yv, xv = np.meshgrid(y, x)
+
+    funcs = np.zeros((x.size, y.size, len(charges)))
+
+    for i, (charge, density) in enumerate(zip(charges, densities)):
+        funcs[:, :, i] = charge * density(xv, yv)
+
+    return funcs.sum(axis=2)
+
+
+ncell = 100
+dl = 0.2
+L = ncell * dl
+ts = 0.01
+masses = (2, 3)
+charges = (1, 2)
+
+
+def densityMain_1d(x):
+    return 1.0
+
+
+def densityBeam_1d(x):
+    u = x / L - 0.5
+    return np.exp(-(u**2))
+
+
+def bx_1d(x):
+    return 1.0
+
+
+def by_1d(x):
+    return 0.0
+
+
+def bz_1d(x):
+    return 0.0
+
+
+def v0_1d(x):
+    return 0.0
+
+
+def vth_1d(x):
+    return np.sqrt(1.0)
+
+
+def config_1d():
+
+    sim = ph.Simulation(
+        smallest_patch_size=20,
+        largest_patch_size=60,
+        time_step=ts,
+        time_step_nbr=1,
+        boundary_types="periodic",
+        cells=ncell,
+        dl=dl,
+        diag_options={
+            "format": "phareh5",
+            "options": {"dir": "test_densities_1d", "mode": "overwrite"},
+        },
+    )
+
+    v_pop = {
+        "vbulkx": v0_1d,
+        "vbulky": v0_1d,
+        "vbulkz": v0_1d,
+        "vthx": vth_1d,
+        "vthy": vth_1d,
+        "vthz": vth_1d,
+    }
+
+    ph.MaxwellianFluidModel(
+        bx=bx_1d,
+        by=by_1d,
+        bz=bz_1d,
+        main={
+            "mass": masses[0],
+            "charge": charges[0],
+            "density": densityMain_1d,
+            "nbr_part_per_cell": 100,
+            **v_pop,
+        },
+        beam={
+            "mass": masses[1],
+            "charge": charges[1],
+            "density": densityBeam_1d,
+            "nbr_part_per_cell": 100,
+            **v_pop,
+        },
+    )
+
+    ph.ElectronModel(closure="isothermal", Te=0.0)
+
+    timestamps = all_timestamps(global_vars.sim)
+
+    for quantity in ["charge_density", "mass_density"]:
+        FluidDiagnostics(quantity=quantity, write_timestamps=timestamps)
+
+    poplist = ["main", "beam"]
+    for pop in poplist:
+        for quantity in ["density", "charge_density"]:
+            FluidDiagnostics(
+                quantity=quantity,
+                write_timestamps=timestamps,
+                population_name=pop,
+            )
+
+    return sim
+
+
+def densityMain_2d(x, y):
+    assert len(x) == len(y)
+    return 1.0 * np.ones_like(x)
+
+
+def densityBeam_2d(x, y):
+    assert len(x) == len(y)
+    u = x / L - 0.5
+    v = y / L - 0.5
+    return np.exp(-(u**2) - v**2)
+
+
+def bx_2d(x, y):
+    return 1.0
+
+
+def by_2d(x, y):
+    return 0.0
+
+
+def bz_2d(x, y):
+    return 0.0
+
+
+def v0_2d(x, y):
+    return 0.0
+
+
+def vth_2d(x, y):
+    return np.sqrt(1.0)
+
+
+def config_2d():
+
+    sim = ph.Simulation(
+        smallest_patch_size=20,
+        largest_patch_size=60,
+        time_step=ts,
+        time_step_nbr=1,
+        boundary_types=("periodic", "periodic"),
+        cells=(ncell, ncell),
+        dl=(dl, dl),
+        diag_options={
+            "format": "phareh5",
+            "options": {"dir": "test_densities_2d", "mode": "overwrite"},
+        },
+    )
+
+    v_pop = {
+        "vbulkx": v0_2d,
+        "vbulky": v0_2d,
+        "vbulkz": v0_2d,
+        "vthx": vth_2d,
+        "vthy": vth_2d,
+        "vthz": vth_2d,
+    }
+
+    ph.MaxwellianFluidModel(
+        bx=bx_2d,
+        by=by_2d,
+        bz=bz_2d,
+        main={
+            "mass": masses[0],
+            "charge": charges[0],
+            "density": densityMain_2d,
+            "nbr_part_per_cell": 100,
+            **v_pop,
+        },
+        beam={
+            "mass": masses[1],
+            "charge": charges[1],
+            "density": densityBeam_2d,
+            "nbr_part_per_cell": 100,
+            **v_pop,
+        },
+    )
+
+    ph.ElectronModel(closure="isothermal", Te=0.0)
+
+    timestamps = all_timestamps(global_vars.sim)
+
+    for quantity in ["charge_density", "mass_density"]:
+        FluidDiagnostics(quantity=quantity, write_timestamps=timestamps)
+
+    poplist = ["main", "beam"]
+    for pop in poplist:
+        for quantity in ["density", "charge_density"]:
+            FluidDiagnostics(
+                quantity=quantity,
+                write_timestamps=timestamps,
+                population_name=pop,
+            )
+
+    return sim
+
+
+def main():
+
+    Simulator(config_1d()).run().reset()
+    ph.global_vars.sim = None
+    Simulator(config_2d()).run().reset()
+
+    def noise_level(h, H):
+        noises = list()
+        for lvl_h, lvl_H in zip(h.levels(time).values(), H.levels(time).values()):
+            for patch_h, patch_H in zip(lvl_h.patches, lvl_H.patches):
+                pd_h = patch_h.patch_datas["value"]
+                pd_H = patch_H.patch_datas["value"]
+
+                dset_h = pd_h[pd_h.box]
+                dset_H = pd_H[pd_H.box]
+
+                noise = np.std(dset_h - dset_H)
+                # print(f"noise = ", noise)
+                noises.append(noise)
+        return noises
+
+    fig, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(3, 2, figsize=(6, 8))
+
+    # 1d stuffs
+    run_path = os.path.join(os.curdir, "test_densities_1d")
+    time = 0.0
+    r = Run(run_path)
+
+    h1 = r.GetMassDensity(time)
+    h2 = r.GetNi(time)
+
+    H1 = hierarchy_from(
+        hier=h1,
+        func=ions_mass_density_func1d,
+        masses=masses,
+        densities=(densityMain_1d, densityBeam_1d),
+    )
+    H2 = hierarchy_from(
+        hier=h2,
+        func=ions_charge_density_func1d,
+        charges=charges,
+        densities=(densityMain_1d, densityBeam_1d),
+    )
+
+    assert np.mean(noise_level(h1, H1)) < 0.18
+    assert np.mean(noise_level(h2, H2)) < 0.12
+
+    cycle = plt.rcParams["axes.prop_cycle"].by_key()["color"]
+
+    h1.plot(ax=ax1, ls="-", lw=2.0, color=cycle[0])
+    H1.plot(ax=ax1, ls="-", lw=2.0, color=cycle[1])
+
+    h2.plot(ax=ax2, ls="-", lw=2.0, color=cycle[0])
+    H2.plot(ax=ax2, ls="-", lw=2.0, color=cycle[1])
+
+    ax1.set_title("mass density : 1d")
+    ax2.set_title("charge density : 1d")
+
+    # 2d stuffs
+    run_path = os.path.join(os.curdir, "test_densities_2d")
+    time = 0.0
+    r = Run(run_path)
+
+    h1 = r.GetMassDensity(time)
+    h2 = r.GetNi(time)
+
+    H1 = hierarchy_from(
+        hier=h1,
+        func=ions_mass_density_func2d,
+        masses=masses,
+        densities=(densityMain_2d, densityBeam_2d),
+    )
+    H2 = hierarchy_from(
+        hier=h2,
+        func=ions_charge_density_func2d,
+        charges=charges,
+        densities=(densityMain_2d, densityBeam_2d),
+    )
+
+    assert np.mean(noise_level(h1, H1)) < 0.18
+    assert np.mean(noise_level(h2, H2)) < 0.12
+
+    cmap = mpl.colormaps["viridis"]
+
+    h1.plot(ax=ax3, vmin=3.75, vmax=5, cmap=cmap, title="computed mass density : 2d")
+    H1.plot(ax=ax4, vmin=3.75, vmax=5, cmap=cmap, title="expected mass density : 2d")
+    h2.plot(ax=ax5, vmin=2.2, vmax=3, cmap=cmap, title="computed charge density : 2d")
+    H2.plot(ax=ax6, vmin=2.2, vmax=3, cmap=cmap, title="expected charge density : 2d")
+
+    plt.tight_layout()
+    plt.savefig("nCheck.pdf", dpi=300)
+
+
+from tests.simulator import SimulatorTest
+
+
+class DensityInitTest(SimulatorTest):
+    def __init__(self, *args, **kwargs):
+        super(DensityInitTest, self).__init__(*args, **kwargs)
+        self.simulator = None
+
+    def tearDown(self):
+        super(DensityInitTest, self).tearDown()
+        if self.simulator is not None:
+            self.simulator.reset()
+        self.simulator = None
+        ph.global_vars.sim = None
+
+    def test_run(self):
+        self.register_diag_dir_for_cleanup("test_densities_1d")
+        self.register_diag_dir_for_cleanup("test_densities_2d")
+        main()
+        return self
+
+
+if __name__ == "__main__":
+    DensityInitTest().test_run().tearDown()
diff --git a/tests/simulator/refinement/test_2d_10_core.py b/tests/simulator/refinement/test_2d_10_core.py
index dfb88703a..114b7296b 100644
--- a/tests/simulator/refinement/test_2d_10_core.py
+++ b/tests/simulator/refinement/test_2d_10_core.py
@@ -101,7 +101,7 @@ def vthz(x, y):
     for quantity in ["E", "B"]:
         ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
-    for quantity in ["density", "bulkVelocity"]:
+    for quantity in ["charge_density", "bulkVelocity"]:
         ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
     return sim
diff --git a/tests/simulator/refinement/test_2d_2_core.py b/tests/simulator/refinement/test_2d_2_core.py
index af017f60c..c64d488fc 100644
--- a/tests/simulator/refinement/test_2d_2_core.py
+++ b/tests/simulator/refinement/test_2d_2_core.py
@@ -101,7 +101,7 @@ def vthz(x, y):
     for quantity in ["E", "B"]:
         ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
-    for quantity in ["density", "bulkVelocity"]:
+    for quantity in ["charge_density", "bulkVelocity"]:
         ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
     return sim
@@ -167,12 +167,7 @@ def post_advance_1(new_time):
         errors = test.base_test_overlaped_fields_are_equal(L0L1_datahier, new_time)
         #        errors = test.base_test_field_level_ghosts_via_subcycles_and_coarser_interpolation(L0_datahier, L0L1_datahier)
         if isinstance(errors, list):
-            extra_collections += [
-                {
-                    "boxes": errors,
-                    "facecolor": "black",
-                }
-            ]
+            extra_collections += [{"boxes": errors, "facecolor": "black"}]
         make_fig(L0L1_datahier, L0L1_diags.split("/")[-1], 1, extra_collections)
 
 
diff --git a/tests/simulator/test_advance.py b/tests/simulator/test_advance.py
index 569f2c873..a5b40bcc6 100644
--- a/tests/simulator/test_advance.py
+++ b/tests/simulator/test_advance.py
@@ -152,7 +152,7 @@ def vthz(*xyz):
         for quantity in ["E", "B"]:
             ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
-        for quantity in ["density", "bulkVelocity"]:
+        for quantity in ["charge_density", "bulkVelocity"]:
             FluidDiagnostics(
                 quantity=quantity,
                 write_timestamps=timestamps,
@@ -222,7 +222,7 @@ def vthz(*xyz):
                 h5_filename=diag_outputs + "/EM_B.h5", hier=eb_hier
             )
             mom_hier = hierarchy_from(
-                h5_filename=diag_outputs + "/ions_density.h5", hier=eb_hier
+                h5_filename=diag_outputs + "/ions_charge_density.h5", hier=eb_hier
             )
             mom_hier = hierarchy_from(
                 h5_filename=diag_outputs + "/ions_bulkVelocity.h5", hier=mom_hier
@@ -231,7 +231,7 @@ def vthz(*xyz):
 
         if qty == "moments" or qty == "fields":
             mom_hier = hierarchy_from(
-                h5_filename=diag_outputs + "/ions_density.h5", hier=eb_hier
+                h5_filename=diag_outputs + "/ions_charge_density.h5", hier=eb_hier
             )
             mom_hier = hierarchy_from(
                 h5_filename=diag_outputs + "/ions_bulkVelocity.h5", hier=mom_hier
@@ -284,12 +284,10 @@ def base_test_overlaped_fields_are_equal(self, datahier, coarsest_time):
                         ]
 
                     try:
-                        # empirical max absolute observed 5.2e-15
-                        # https://hephaistos.lpp.polytechnique.fr/teamcity/buildConfiguration/Phare_Phare_BuildGithubPrClang/78544
-                        # seems correct considering ghosts are filled with schedules
-                        # involving linear/spatial interpolations and so on where
-                        # rounding errors may occur.... setting atol to 5.5e-15
-                        assert_fp_any_all_close(slice1, slice2, atol=5.5e-15, rtol=0)
+                        # empirical max absolute observed 2.5e-14
+                        # set after introducing the new charge density calculation
+                        # different from the mass density one.
+                        assert_fp_any_all_close(slice1, slice2, atol=2.5e-14, rtol=0)
                         checks += 1
                     except AssertionError as e:
                         print("AssertionError", pd1.name, e)
diff --git a/tests/simulator/test_initialization.py b/tests/simulator/test_initialization.py
index 2ea135909..de00a9ab3 100644
--- a/tests/simulator/test_initialization.py
+++ b/tests/simulator/test_initialization.py
@@ -170,7 +170,7 @@ def vthz(*xyz):
                 quantity=quantity, write_timestamps=np.zeros(time_step_nbr)
             )
 
-        for quantity in ["density", "bulkVelocity"]:
+        for quantity in ["charge_density", "bulkVelocity"]:
             FluidDiagnostics(
                 quantity=quantity, write_timestamps=np.zeros(time_step_nbr)
             )
@@ -232,7 +232,9 @@ def vthz(*xyz):
             return particle_hier
 
         if qty == "moments":
-            mom_hier = hierarchy_from(h5_filename=diag_outputs + "/ions_density.h5")
+            mom_hier = hierarchy_from(
+                h5_filename=diag_outputs + "/ions_charge_density.h5"
+            )
             mom_hier = hierarchy_from(
                 h5_filename=diag_outputs + "/ions_bulkVelocity.h5", hier=mom_hier
             )
diff --git a/tests/simulator/test_run.py b/tests/simulator/test_run.py
index f935e9768..0b8accc37 100644
--- a/tests/simulator/test_run.py
+++ b/tests/simulator/test_run.py
@@ -142,7 +142,7 @@ def vthz(x, y):
 
     for quantity in ["E", "B"]:
         ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps)
-    for quantity in ["density", "bulkVelocity"]:
+    for quantity in ["charge_density", "bulkVelocity"]:
         ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
     pop = "protons"