entity-toolkit · haykh · May 21, 2026 · May 18, 2026 · May 18, 2026 · May 18, 2026
diff --git a/CITATION b/CITATION
@@ -41,3 +41,23 @@ archivePrefix = {arXiv},
       adsnote = {Provided by the SAO/NASA Astrophysics Data System}
 }
 ```
+
+For the higher-order shape functions and generalized field stencils, please cite the following paper:
+
+```latex
+@ARTICLE{EntityHighOrder_2026,
+       author = {{B{\"o}ss}, Ludwig M. and {Vanthieghem}, Arno and {Hakobyan}, Hayk and {Gorbunov}, Evgeny A. and {Caprioli}, Damiano},
+        title = "{Entity -- Hardware-agnostic Particle-in-Cell Code for Plasma Astrophysics. III: Higher-order shape functions \& generalized field stencils}",
+      journal = {arXiv e-prints},
+     keywords = {High Energy Astrophysical Phenomena},
+         year = 2026,
+        month = may,
+          eid = {arXiv:2605.15260},
+        pages = {arXiv:2605.15260},
+archivePrefix = {arXiv},
+       eprint = {2605.15260},
+ primaryClass = {astro-ph.HE},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2026arXiv260515260B},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+```
diff --git a/README.md b/README.md
@@ -12,12 +12,16 @@ Our [detailed documentation](https://entity-toolkit.github.io/) includes everyth
 
 ## News
 
+- [May 2026]: our **method paper** for higher-order shape functions and generalized field stencils is [online](https://ui.adsabs.harvard.edu/abs/2026arXiv260515260B/abstract)
+- [Apr 2026]: user-defined **custom particle update** functionality [PR #198](https://github.com/entity-toolkit/entity/pull/198)
+- [Apr 2026]: **moving window** [PR #196](https://github.com/entity-toolkit/entity/pull/196)
+- [Apr 2026]: **built-in piston** [PR #192](https://github.com/entity-toolkit/entity/pull/192)
 - [Mar 2026]: **single-particle emission** [PR #174](https://github.com/entity-toolkit/entity/pull/174) and [PR #188](https://github.com/entity-toolkit/entity/pull/188)
 - [Mar 2026]: **spatial sorting of particles** [PR #181](https://github.com/entity-toolkit/entity/pull/181)
 - [Mar 2026]: **external EM & force fields** [PR #183](https://github.com/entity-toolkit/entity/pull/183)
 - [Mar 2026]: **examples** of most commonly used custom patterns (see `pgens/examples`)
-- [Dec 2025]: **high-order** shape functions [PR #109](https://github.com/entity-toolkit/entity/pull/109) and advanced field stencils [PR #103](https://github.com/entity-toolkit/entity/pull/103) are now supported.
-- [Dec 2025]: **particle tracking** is now fully supported via [PR #144](https://github.com/entity-toolkit/entity/pull/144).
+- [Dec 2025]: **high-order** shape functions [PR #109](https://github.com/entity-toolkit/entity/pull/109) and advanced field stencils [PR #103](https://github.com/entity-toolkit/entity/pull/103) are now supported
+- [Dec 2025]: **particle tracking** is now fully supported via [PR #144](https://github.com/entity-toolkit/entity/pull/144)
 - [Nov 2025]: our **method papers** are online: [Special relativistic module](https://ui.adsabs.harvard.edu/abs/2025arXiv251117710H/abstract), [GR module](https://ui.adsabs.harvard.edu/abs/2025arXiv251117701G/abstract)!
 
 ## Citation

diff --git a/input.example.toml b/input.example.toml
@@ -514,10 +514,6 @@
     #   @type: array<string>
     #   @default: []
     custom = ""
-    # Smoothing window for the output of moments ("Rho", "Charge", "T", ...)
-    #   @type: ushort
-    #   @default: 0
-    mom_smooth = ""
     # Number of timesteps between field outputs
     #   @type: uint
     #   @default: 0
@@ -537,6 +533,19 @@
     #   @note: If a scalar is given, it is applied to all directions
     downsampling = ""
 
+    [output.fields.smoothing]
+      # Smoothing order for the output of moments ("Rho", "Charge", "T", ...)
+      #   @type: ushort
+      #   @default: 0
+      order = ""
+      # Smoothing algorithm 
+      #   @type: string
+      #   @enum: "const", "spline"
+      #   @default: "spline"
+      #   @note: When using "spline", `order` corresponds to the order of the polynomial used
+      #   @note: When using "const", the smoothing window is `ceil(order / 2)` in both directions
+      method = ""
+
   [output.particles]
     # Toggle for the particles output
     #   @type: bool

diff --git a/pgens/accretion/pgen.hpp b/pgens/accretion/pgen.hpp
@@ -11,9 +11,9 @@
 
 #include "archetypes/energy_dist.h"
 #include "archetypes/particle_injector.h"
+#include "archetypes/utils.h"
 #include "framework/domain/metadomain.h"
 #include "framework/parameters/parameters.h"
-#include "kernels/particle_moments.hpp"
 
 namespace user {
   using namespace ntt;
@@ -117,39 +117,18 @@ namespace user {
       std::copy(xi_min.begin(), xi_min.end(), x_min);
       std::copy(xi_max.begin(), xi_max.end(), x_max);
 
-      std::vector<unsigned short> specs {};
+      std::vector<spidx_t> specs {};
       for (auto& sp : domain_ptr->species) {
         if (sp.mass() > 0) {
           specs.push_back(sp.index());
         }
       }
 
-      Kokkos::deep_copy(density, ZERO);
-      auto  scatter_buff = Kokkos::Experimental::create_scatter_view(density);
-      // some parameters
-      auto& mesh         = domain_ptr->mesh;
-      const auto use_weights = params.template get<bool>(
-        "particles.use_weights");
-      const auto ni2 = mesh.n_active(in::x2);
-
-      for (const auto& sp : specs) {
-        auto& prtl_spec = domain_ptr->species[sp - 1];
-        // clang-format off
-        Kokkos::parallel_for(
-          "ComputeMoments",
-          prtl_spec.rangeActiveParticles(),
-          kernel::ParticleMoments_kernel<SimEngine::GRPIC, M, FldsID::Rho, 3>({}, scatter_buff, 0u,
-                                                               prtl_spec.i1, prtl_spec.i2, prtl_spec.i3,
-                                                               prtl_spec.dx1, prtl_spec.dx2, prtl_spec.dx3,
-                                                               prtl_spec.ux1, prtl_spec.ux2, prtl_spec.ux3,
-                                                               prtl_spec.phi, prtl_spec.weight, prtl_spec.tag,
-                                                               prtl_spec.mass(), prtl_spec.charge(),
-                                                               use_weights,
-                                                               metric, mesh.flds_bc(),
-                                                               ni2, inv_n0, ZERO));
-        // clang-format on
-      }
-      Kokkos::Experimental::contribute(density, scatter_buff);
+      arch::ComputeMomentWithSpecies<SimEngine::GRPIC, M, FldsID::Rho, 3>(
+        params,
+        *domain_ptr,
+        specs,
+        density);
     }
 
     Inline auto sigma_crit(const coord_t<M::Dim>& x_Ph) const -> bool {

diff --git a/pgens/reconnection/pgen.hpp b/pgens/reconnection/pgen.hpp
@@ -13,7 +13,6 @@
 #include "archetypes/spatial_dist.h"
 #include "archetypes/utils.h"
 #include "framework/domain/metadomain.h"
-#include "kernels/particle_moments.hpp"
 
 namespace user {
   using namespace ntt;
@@ -247,35 +246,11 @@ namespace user {
         inj_box_down.push_back(Range::All);
       }
 
-      {
-        // compute density of species #1 and #2
-        const auto use_weights = params.template get<bool>(
-          "particles.use_weights");
-        const auto ni2    = domain.mesh.n_active(in::x2);
-        const auto inv_n0 = ONE / params.template get<real_t>("scales.n0");
-
-        auto scatter_buff = Kokkos::Experimental::create_scatter_view(
-          domain.fields.buff);
-        Kokkos::deep_copy(domain.fields.buff, ZERO);
-        for (const auto sp : std::vector<spidx_t> { 1, 2 }) {
-          const auto& prtl_spec = domain.species[sp - 1];
-          // clang-format off
-          Kokkos::parallel_for(
-            "ComputeMoments",
-            prtl_spec.rangeActiveParticles(),
-            kernel::ParticleMoments_kernel<S, M, FldsID::N, 3>({}, scatter_buff, 0u,
-                                                              prtl_spec.i1, prtl_spec.i2, prtl_spec.i3,
-                                                              prtl_spec.dx1, prtl_spec.dx2, prtl_spec.dx3,
-                                                              prtl_spec.ux1, prtl_spec.ux2, prtl_spec.ux3,
-                                                              prtl_spec.phi, prtl_spec.weight, prtl_spec.tag,
-                                                              prtl_spec.mass(), prtl_spec.charge(),
-                                                              use_weights,
-                                                              domain.mesh.metric, domain.mesh.flds_bc(),
-                                                              ni2, inv_n0, 0u));
-          // clang-format on
-        }
-        Kokkos::Experimental::contribute(domain.fields.buff, scatter_buff);
-      }
+      // compute density of species #1 and #2
+      arch::ComputeMomentWithSpecies<S, M, FldsID::Rho, 3>(params,
+                                                           domain,
+                                                           { 1, 2 },
+                                                           domain.fields.buff);
 
       const auto replenish_sdist = arch::spatial_dist::ReplenishUniform<M, 3>(
         domain.mesh.metric,

diff --git a/src/archetypes/moving_window.h b/src/archetypes/moving_window.h
@@ -28,10 +28,10 @@ namespace arch {
     const BCTags    tags;
     const int       window_shift;
 
-    FieldShift_kernel(ndfield_t<D, 6> Fld,
-                      ndfield_t<D, 6> backup_Fld,
-                      const int       window_shift,
-                      BCTags          tags)
+    FieldShift_kernel(ndfield_t<D, 6>&       Fld,
+                      const ndfield_t<D, 6>& backup_Fld,
+                      int                    window_shift,
+                      BCTags                 tags)
       : Fld { Fld }
       , backup_Fld { backup_Fld }
       , window_shift { window_shift }

diff --git a/src/archetypes/utils.h b/src/archetypes/utils.h
@@ -124,21 +124,24 @@ namespace arch {
    * @param components Vector of field components to compute (e.g. {} for N, {0}
    * {1} {2} for V, {0, 1} for T, etc., default: empty, i.e. scalar)
    * @param buffer_idx Index of the field component in the buffer to save the result to
-   * @param window Window size for smoothing (in number of cells, default: 0, i.e. no smoothing)
+   * @param smoothing_order Smoothing order
+   * @param smoothing_method Smoothing algorithm (e.g. SPLINE, CONST, etc.)
    *
    * @tparam S Simulation engine type
    * @tparam M Metric type
    * @tparam F Field ID for the moment to compute (e.g. FldsID::N, FldsID::T, etc.)
    * @tparam N Last dimension of the buffer (e.g. 3 or 6)
    */
   template <SimEngine::type S, MetricClass M, FldsID::type F, int N>
-  inline void ComputeMomentWithSpecies(const SimulationParams&     params,
-                                       Domain<S, M>&               domain,
-                                       const std::vector<spidx_t>& species,
-                                       ndfield_t<M::Dim, N>&       buffer,
-                                       const std::vector<uint8_t>& components = {},
-                                       idx_t          buffer_idx = 0u,
-                                       unsigned short window     = 0u) {
+  inline void ComputeMomentWithSpecies(
+    const SimulationParams&     params,
+    Domain<S, M>&               domain,
+    const std::vector<spidx_t>& species,
+    ndfield_t<M::Dim, N>&       buffer,
+    const std::vector<uint8_t>& components      = {},
+    idx_t                       buffer_idx      = 0u,
+    uint8_t                     smoothing_order = 0u,
+    OutputSmoothingTypeFlag smoothing_method    = OutputSmoothingType::SPLINE) {
     const auto ni2         = domain.mesh.n_active(in::x2);
     const auto inv_n0      = ONE / params.template get<real_t>("scales.n0");
     const auto use_weights = params.template get<bool>("particles.use_weights");
@@ -153,26 +156,14 @@ namespace arch {
         kernel::ParticleMoments_kernel<S, M, F, N>(components,
                                                    scatter_buff,
                                                    buffer_idx,
-                                                   prtl_spec.i1,
-                                                   prtl_spec.i2,
-                                                   prtl_spec.i3,
-                                                   prtl_spec.dx1,
-                                                   prtl_spec.dx2,
-                                                   prtl_spec.dx3,
-                                                   prtl_spec.ux1,
-                                                   prtl_spec.ux2,
-                                                   prtl_spec.ux3,
-                                                   prtl_spec.phi,
-                                                   prtl_spec.weight,
-                                                   prtl_spec.tag,
-                                                   prtl_spec.mass(),
-                                                   prtl_spec.charge(),
+                                                   prtl_spec,
                                                    use_weights,
                                                    domain.mesh.metric,
                                                    domain.mesh.flds_bc(),
                                                    ni2,
                                                    inv_n0,
-                                                   window));
+                                                   smoothing_order,
+                                                   smoothing_method));
     }
     Kokkos::Experimental::contribute(buffer, scatter_buff);
   }

diff --git a/src/engines/grpic/fields_bcs.h b/src/engines/grpic/fields_bcs.h
@@ -39,13 +39,13 @@ namespace ntt {
     };
 
     template <GRMetricClass M, PGenClass<SimEngine::GRPIC, M> PG>
-    void MatchFieldsIn(dir::direction_t<M::Dim>     direction,
-                       Domain<SimEngine::GRPIC, M>& domain,
-                       const Grid<M::Dim>&          global_grid,
-                       const PG&                    pgen,
-                       const SimulationParams&      params,
-                       BCTags                       tags,
-                       const gr_bc&                 g) {
+    void MatchFieldsIn(const dir::direction_t<M::Dim>& direction,
+                       Domain<SimEngine::GRPIC, M>&    domain,
+                       const Grid<M::Dim>&             global_grid,
+                       const PG&                       pgen,
+                       const SimulationParams&         params,
+                       BCTags                          tags,
+                       const gr_bc&                    g) {
       /**
        * match boundaries
        */
@@ -138,11 +138,11 @@ namespace ntt {
     }
 
     template <GRMetricClass M>
-    void HorizonFieldsIn(dir::direction_t<M::Dim>     direction,
-                         Domain<SimEngine::GRPIC, M>& domain,
-                         const SimulationParams&      params,
-                         BCTags                       tags,
-                         const gr_bc&                 g) {
+    void HorizonFieldsIn(const dir::direction_t<M::Dim>& direction,
+                         Domain<SimEngine::GRPIC, M>&    domain,
+                         const SimulationParams&         params,
+                         BCTags                          tags,
+                         const gr_bc&                    g) {
       /**
        * open boundaries
        */
@@ -176,9 +176,9 @@ namespace ntt {
     }
 
     template <GRMetricClass M>
-    void AxisFieldsIn(dir::direction_t<M::Dim>     direction,
-                      Domain<SimEngine::GRPIC, M>& domain,
-                      BCTags                       tags) {
+    void AxisFieldsIn(const dir::direction_t<M::Dim>& direction,
+                      Domain<SimEngine::GRPIC, M>&    domain,
+                      BCTags                          tags) {
       /**
        * axis boundaries
        */
@@ -220,10 +220,10 @@ namespace ntt {
     }
 
     template <GRMetricClass M>
-    void CustomFieldsIn(dir::direction_t<M::Dim>     direction,
-                        Domain<SimEngine::GRPIC, M>& domain,
-                        BCTags                       tags,
-                        const gr_bc&                 g) {
+    void CustomFieldsIn(const dir::direction_t<M::Dim>& direction,
+                        Domain<SimEngine::GRPIC, M>&    domain,
+                        BCTags                          tags,
+                        const gr_bc&                    g) {
       (void)direction;
       (void)domain;
       (void)tags;

diff --git a/src/engines/srpic/particles_bcs.h b/src/engines/srpic/particles_bcs.h
@@ -88,21 +88,20 @@ namespace ntt {
           if (prtl_spec.npart() == 0) {
             continue;
           }
-          // clang-format off
           Kokkos::parallel_for(
             "ComputeMoments",
             prtl_spec.rangeActiveParticles(),
             kernel::ParticleMoments_kernel<SimEngine::SRPIC, M, FldsID::Rho, 6>(
-              {}, scatter_bckp, 0,
-              prtl_spec.i1, prtl_spec.i2, prtl_spec.i3,
-              prtl_spec.dx1, prtl_spec.dx2, prtl_spec.dx3,
-              prtl_spec.ux1, prtl_spec.ux2, prtl_spec.ux3,
-              prtl_spec.phi, prtl_spec.weight, prtl_spec.tag,
-              prtl_spec.mass(), prtl_spec.charge(),
+              {},
+              scatter_bckp,
+              0,
+              prtl_spec,
               use_weights,
-              domain.mesh.metric, domain.mesh.flds_bc(),
-              ni2, inv_n0, 0));
-          // clang-format on
+              domain.mesh.metric,
+              domain.mesh.flds_bc(),
+              ni2,
+              inv_n0,
+              0u));
           prtl_spec.set_unsorted();
         }
         Kokkos::Experimental::contribute(domain.fields.bckp, scatter_bckp);

diff --git a/src/framework/domain/metadomain_io.cpp b/src/framework/domain/metadomain_io.cpp
@@ -1,3 +1,4 @@
+#include "defaults.h"
 #include "enums.h"
 #include "global.h"
 
@@ -146,25 +147,27 @@ namespace ntt {
     const auto use_weights = params.template get<bool>("particles.use_weights");
     const auto ni2         = mesh.n_active(in::x2);
     const auto inv_n0      = ONE / params.template get<real_t>("scales.n0");
-    const auto window      = params.template get<unsigned short>(
-      "output.fields.mom_smooth");
+    const auto smooth_order = params.template get<unsigned short>(
+      "output.fields.smoothing.order");
+    const auto smooth_method = OutputSmoothingType::from_string(
+      params.template get<std::string>("output.fields.smoothing.method"));
 
     for (const auto& sp : specs) {
       auto& prtl_spec = prtl_species[sp - 1];
-      // clang-format off
       Kokkos::parallel_for(
         "ComputeMoments",
         prtl_spec.rangeActiveParticles(),
-        kernel::ParticleMoments_kernel<S, M, F, 6>(components, scatter_buff, buff_idx,
-                                                   prtl_spec.i1, prtl_spec.i2, prtl_spec.i3,
-                                                   prtl_spec.dx1, prtl_spec.dx2, prtl_spec.dx3,
-                                                   prtl_spec.ux1, prtl_spec.ux2, prtl_spec.ux3,
-                                                   prtl_spec.phi, prtl_spec.weight, prtl_spec.tag,
-                                                   prtl_spec.mass(), prtl_spec.charge(),
+        kernel::ParticleMoments_kernel<S, M, F, 6>(components,
+                                                   scatter_buff,
+                                                   buff_idx,
+                                                   prtl_spec,
                                                    use_weights,
-                                                   mesh.metric, mesh.flds_bc(),
-                                                   ni2, inv_n0, window));
-      // clang-format on
+                                                   mesh.metric,
+                                                   mesh.flds_bc(),
+                                                   ni2,
+                                                   inv_n0,
+                                                   smooth_order,
+                                                   smooth_method));
     }
     Kokkos::Experimental::contribute(buffer, scatter_buff);
   }