diff --git a/README.md b/README.md
index 41a5fe47e..d6f4597f5 100644
--- a/README.md
+++ b/README.md
@@ -12,6 +12,8 @@ Our [detailed documentation](https://entity-toolkit.github.io/) includes everyth
 
 ## Core developers (alphabetical)
 
+👀 __Yangyang Cai__ {[@StaticObserver](https://github.com/StaticObserver): GRPIC}
+
 💁‍♂️ __Alexander Chernoglazov__ {[@SChernoglazov](https://github.com/SChernoglazov): PIC}
 
 🍵 __Benjamin Crinquand__ {[@bcrinquand](https://github.com/bcrinquand): GRPIC, cubed-sphere}
diff --git a/dev/runners/Dockerfile.runner.cuda b/dev/runners/Dockerfile.runner.cuda
index 4ff132990..6e0b0755c 100644
--- a/dev/runners/Dockerfile.runner.cuda
+++ b/dev/runners/Dockerfile.runner.cuda
@@ -59,14 +59,17 @@ ARG HOME=/home/$USER
 WORKDIR $HOME
 
 # gh runner
+ARG RUNNER_VERSION=2.317.0
 RUN mkdir actions-runner
 WORKDIR $HOME/actions-runner
 
-RUN --mount=type=secret,id=ghtoken \
-    curl -o actions-runner-linux-x64-2.317.0.tar.gz \
-         -L https://github.com/actions/runner/releases/download/v2.317.0/actions-runner-linux-x64-2.317.0.tar.gz && \
-    tar xzf ./actions-runner-linux-x64-2.317.0.tar.gz && \
-    sudo ./bin/installdependencies.sh && \
-    ./config.sh --url https://github.com/entity-toolkit/entity --token "$(sudo cat /run/secrets/ghtoken)" --labels nvidia-gpu
+RUN curl -o actions-runner-linux-x64-${RUNNER_VERSION}.tar.gz \
+         -L https://github.com/actions/runner/releases/download/v${RUNNER_VERSION}/actions-runner-linux-x64-${RUNNER_VERSION}.tar.gz && \
+    tar xzf ./actions-runner-linux-x64-${RUNNER_VERSION}.tar.gz && \
+    sudo ./bin/installdependencies.sh
 
-ENTRYPOINT ["./run.sh"]
+ADD start.sh start.sh
+RUN sudo chown $USER:$USER start.sh && \
+    sudo chmod +x start.sh
+
+ENTRYPOINT ["./start.sh"]
diff --git a/legacy/src/framework/utils/particle_injectors.hpp b/legacy/src/framework/utils/particle_injectors.hpp
index 21e3dad72..c275f170a 100644
--- a/legacy/src/framework/utils/particle_injectors.hpp
+++ b/legacy/src/framework/utils/particle_injectors.hpp
@@ -165,7 +165,7 @@ namespace ntt {
    * @brief Volumetrically uniform particle injector parallelized over particles.
    * @tparam D dimension.
    * @tparam S simulation engine.
-   * @tparam EnDist energy distribution [default = ColdDist].
+   * @tparam EnDist energy distribution [default = Cold].
    *
    * @param params simulation parameters.
    * @param mblock meshblock.
@@ -174,7 +174,7 @@ namespace ntt {
    * @param region region to inject particles as a list of coordinates [optional].
    * @param time current time [optional].
    */
-  template <Dimension D, SimulationEngine S, template <Dimension, SimulationEngine> class EnDist = ColdDist>
+  template <Dimension D, SimulationEngine S, template <Dimension, SimulationEngine> class EnDist = Cold>
   inline void InjectUniform(const SimulationParams& params,
                             Meshblock<D, S>&        mblock,
                             const std::vector<int>& species,
@@ -613,8 +613,8 @@ namespace ntt {
    * @brief Particle injector parallelized by cells in a volume.
    * @tparam D dimension.
    * @tparam S simulation engine.
-   * @tparam EnDist energy distribution [default = ColdDist].
-   * @tparam SpDist spatial distribution [default = UniformDist].
+   * @tparam EnDist energy distribution [default = Cold].
+   * @tparam SpDist spatial distribution [default = Uniform].
    * @tparam InjCrit injection criterion [default = NoCriterion].
    *
    * @param params simulation parameters.
@@ -626,8 +626,8 @@ namespace ntt {
    */
   template <Dimension        D,
             SimulationEngine S,
-            template <Dimension, SimulationEngine> class EnDist  = ColdDist,
-            template <Dimension, SimulationEngine> class SpDist  = UniformDist,
+            template <Dimension, SimulationEngine> class EnDist  = Cold,
+            template <Dimension, SimulationEngine> class SpDist  = Uniform,
             template <Dimension, SimulationEngine> class InjCrit = NoCriterion>
   inline void InjectInVolume(const SimulationParams& params,
                              Meshblock<D, S>&        mblock,
@@ -928,7 +928,7 @@ namespace ntt {
    * @brief ... up to certain number density.
    * @tparam D dimension.
    * @tparam S simulation engine.
-   * @tparam EnDist energy distribution [default = ColdDist].
+   * @tparam EnDist energy distribution [default = Cold].
    * @tparam InjCrit injection criterion [default = NoCriterion].
    *
    * @param params simulation parameters.
@@ -940,7 +940,7 @@ namespace ntt {
    */
   template <Dimension        D,
             SimulationEngine S,
-            template <Dimension, SimulationEngine> class EnDist  = ColdDist,
+            template <Dimension, SimulationEngine> class EnDist  = Cold,
             template <Dimension, SimulationEngine> class InjCrit = NoCriterion>
   inline void InjectNonUniform(const SimulationParams&      params,
                                Meshblock<D, S>&             mblock,
diff --git a/setups/srpic/turbulence/pgen.hpp b/setups/srpic/turbulence/pgen.hpp
index bd8b0ce41..bbd61cc3a 100644
--- a/setups/srpic/turbulence/pgen.hpp
+++ b/setups/srpic/turbulence/pgen.hpp
@@ -185,7 +185,7 @@ namespace user {
         const auto injector = arch::UniformInjector<S, M, arch::Maxwellian>(
           energy_dist,
           { 1, 2 });
-        const real_t ndens = 1.0;
+        const real_t ndens = 0.9;
         arch::InjectUniform<S, M, decltype(injector)>(params,
                                                       local_domain,
                                                       injector,
@@ -193,22 +193,17 @@ namespace user {
       }
 
       {
-        // const auto energy_dist = arch::Maxwellian<S, M>(local_domain.mesh.metric,
-        //                                                 local_domain.random_pool,
-        //                                                 temperature*10);        
-        // // const auto energy_dist = arch::Maxwellian<S, M>(local_domain.mesh.metric,
-        // //                                                 local_domain.random_pool,
-        // //                                                 temperature * 2,
-        // //                                                 10.0,
-        // //                                                 1);
-        // const auto injector = arch::UniformInjector<S, M, arch::Maxwellian>(
-        //   energy_dist,
-        //   { 1, 2 });
-        // const real_t ndens = 0.01;
-        // arch::InjectUniform<S, M, decltype(injector)>(params,
-        //                                               local_domain,
-        //                                               injector,
-        //                                               ndens);
+        const auto energy_dist = arch::PowerlawDist<S, M>(local_domain.mesh.metric,
+                                                        local_domain.random_pool,
+                                                        0.1, 100.0, -3.0);
+        const auto injector = arch::UniformInjector<S, M, arch::PowerlawDist>(
+          energy_dist,
+          { 1, 2 });
+        const real_t ndens = 0.1;
+        arch::InjectUniform<S, M, decltype(injector)>(params,
+                                                      local_domain,
+                                                      injector,
+                                                      ndens);
       }
     }
 
diff --git a/src/archetypes/energy_dist.h b/src/archetypes/energy_dist.h
index c5b9d8a87..e9bc9051a 100644
--- a/src/archetypes/energy_dist.h
+++ b/src/archetypes/energy_dist.h
@@ -3,7 +3,8 @@
  * @brief Defines an archetype for energy distributions
  * @implements
  *   - arch::EnergyDistribution<>
- *   - arch::ColdDist<> : arch::EnergyDistribution<>
+ *   - arch::Cold<> : arch::EnergyDistribution<>
+ *   - arch::Powerlaw<> : arch::EnergyDistribution<>
  *   - arch::Maxwellian<> : arch::EnergyDistribution<>
  * @namespaces:
  *   - arch::
@@ -55,8 +56,8 @@ namespace arch {
   };
 
   template <SimEngine::type S, class M>
-  struct ColdDist : public EnergyDistribution<S, M> {
-    ColdDist(const M& metric) : EnergyDistribution<S, M> { metric } {}
+  struct Cold : public EnergyDistribution<S, M> {
+    Cold(const M& metric) : EnergyDistribution<S, M> { metric } {}
 
     Inline void operator()(const coord_t<M::Dim>&,
                            vec_t<Dim::_3D>& v,
@@ -67,6 +68,63 @@ namespace arch {
     }
   };
 
+  template <SimEngine::type S, class M>
+  struct Powerlaw : public EnergyDistribution<S, M> {
+    using EnergyDistribution<S, M>::metric;
+
+    Powerlaw(const M&              metric,
+             random_number_pool_t& pool,
+             real_t                g_min,
+             real_t                g_max,
+             real_t                pl_ind)
+      : EnergyDistribution<S, M> { metric }
+      , pool { pool }
+      , g_min { g_min }
+      , g_max { g_max }
+      , pl_ind { pl_ind } {}
+
+    Inline void operator()(const coord_t<M::Dim>& x_Code,
+                           vec_t<Dim::_3D>&       v,
+                           unsigned short         sp = 0) const override {
+      auto rand_gen = pool.get_state();
+      auto rand_X1  = Random<real_t>(rand_gen);
+      auto rand_gam = ONE;
+
+      // Power-law distribution from uniform (see https://mathworld.wolfram.com/RandomNumber.html)
+      if (pl_ind != -ONE) {
+        rand_gam += math::pow(
+          math::pow(g_min, ONE + pl_ind) +
+            (-math::pow(g_min, ONE + pl_ind) + math::pow(g_max, ONE + pl_ind)) *
+              rand_X1,
+          ONE / (ONE + pl_ind));
+      } else {
+        rand_gam += math::pow(g_min, ONE - rand_X1) * math::pow(g_max, rand_X1);
+      }
+      auto rand_u  = math::sqrt(SQR(rand_gam) - ONE);
+      auto rand_X2 = Random<real_t>(rand_gen);
+      auto rand_X3 = Random<real_t>(rand_gen);
+      v[0]         = rand_u * (TWO * rand_X2 - ONE);
+      v[2]         = TWO * rand_u * math::sqrt(rand_X2 * (ONE - rand_X2));
+      v[1]         = v[2] * math::cos(constant::TWO_PI * rand_X3);
+      v[2]         = v[2] * math::sin(constant::TWO_PI * rand_X3);
+
+      if constexpr (S == SimEngine::GRPIC) {
+        // convert from the tetrad basis to covariant
+        vec_t<Dim::_3D> v_Hat;
+        v_Hat[0] = v[0];
+        v_Hat[1] = v[1];
+        v_Hat[2] = v[2];
+        metric.template transform<Idx::T, Idx::D>(x_Code, v_Hat, v);
+      }
+
+      pool.free_state(rand_gen);
+    }
+
+  private:
+    const real_t         g_min, g_max, pl_ind;
+    random_number_pool_t pool;
+  };
+
   template <SimEngine::type S, class M>
   struct Maxwellian : public EnergyDistribution<S, M> {
     using EnergyDistribution<S, M>::metric;
diff --git a/src/archetypes/particle_injector.h b/src/archetypes/particle_injector.h
index 4e75003bb..cbcbbd389 100644
--- a/src/archetypes/particle_injector.h
+++ b/src/archetypes/particle_injector.h
@@ -142,7 +142,7 @@ namespace arch {
     };
 
     using energy_dist_t  = Maxwellian<S, M>;
-    using spatial_dist_t = ReplenishDist<S, M, TargetDensityProfile>;
+    using spatial_dist_t = Replenish<S, M, TargetDensityProfile>;
     static_assert(M::is_metric, "M must be a metric class");
     static constexpr bool      is_nonuniform_injector { true };
     static constexpr Dimension D { M::Dim };
diff --git a/src/archetypes/spatial_dist.h b/src/archetypes/spatial_dist.h
index 6c19d44d0..be2836da2 100644
--- a/src/archetypes/spatial_dist.h
+++ b/src/archetypes/spatial_dist.h
@@ -3,8 +3,8 @@
  * @brief Spatial distribution class passed to injectors
  * @implements
  *   - arch::SpatialDistribution<>
- *   - arch::UniformDist<> : arch::SpatialDistribution<>
- *   - arch::ReplenishDist<> : arch::SpatialDistribution<>
+ *   - arch::Uniform<> : arch::SpatialDistribution<>
+ *   - arch::Replenish<> : arch::SpatialDistribution<>
  * @namespace
  *   - arch::
  * @note
@@ -41,8 +41,8 @@ namespace arch {
   };
 
   template <SimEngine::type S, class M>
-  struct UniformDist : public SpatialDistribution<S, M> {
-    UniformDist(const M& metric) : SpatialDistribution<S, M> { metric } {}
+  struct Uniform : public SpatialDistribution<S, M> {
+    Uniform(const M& metric) : SpatialDistribution<S, M> { metric } {}
 
     Inline auto operator()(const coord_t<M::Dim>&) const -> real_t override {
       return ONE;
@@ -50,7 +50,7 @@ namespace arch {
   };
 
   template <SimEngine::type S, class M, class T>
-  struct ReplenishDist : public SpatialDistribution<S, M> {
+  struct Replenish : public SpatialDistribution<S, M> {
     using SpatialDistribution<S, M>::metric;
     const ndfield_t<M::Dim, 6> density;
     const unsigned short       idx;
@@ -58,11 +58,11 @@ namespace arch {
     const T      target_density;
     const real_t target_max_density;
 
-    ReplenishDist(const M&                    metric,
-                  const ndfield_t<M::Dim, 6>& density,
-                  unsigned short              idx,
-                  const T&                    target_density,
-                  real_t                      target_max_density)
+    Replenish(const M&                    metric,
+              const ndfield_t<M::Dim, 6>& density,
+              unsigned short              idx,
+              const T&                    target_density,
+              real_t                      target_max_density)
       : SpatialDistribution<S, M> { metric }
       , density { density }
       , idx { idx }
diff --git a/src/archetypes/tests/CMakeLists.txt b/src/archetypes/tests/CMakeLists.txt
index ceee1edc9..4ffc35322 100644
--- a/src/archetypes/tests/CMakeLists.txt
+++ b/src/archetypes/tests/CMakeLists.txt
@@ -22,4 +22,5 @@ endfunction()
 
 gen_test(energy_dist)
 gen_test(spatial_dist)
-gen_test(field_setter)
\ No newline at end of file
+gen_test(field_setter)
+gen_test(powerlaw)
diff --git a/src/archetypes/tests/powerlaw.cpp b/src/archetypes/tests/powerlaw.cpp
new file mode 100644
index 000000000..dfcb6b247
--- /dev/null
+++ b/src/archetypes/tests/powerlaw.cpp
@@ -0,0 +1,171 @@
+#include "enums.h"
+#include "global.h"
+
+#include "utils/error.h"
+
+#include "metrics/kerr_schild.h"
+#include "metrics/kerr_schild_0.h"
+#include "metrics/minkowski.h"
+#include "metrics/qkerr_schild.h"
+#include "metrics/qspherical.h"
+#include "metrics/spherical.h"
+
+#include "archetypes/energy_dist.h"
+
+#include <Kokkos_Core.hpp>
+
+#include <iostream>
+
+using namespace ntt;
+using namespace metric;
+using namespace arch;
+
+template <class EnrgDist, SimEngine::type S, class M>
+struct Caller {
+  static constexpr auto D = M::Dim;
+
+  Caller(const M& metric, const EnrgDist& dist)
+    : metric { metric }
+    , dist { dist } {}
+
+  Inline void operator()(index_t) const {
+    vec_t<Dim::_3D> vp { ZERO };
+    coord_t<D>      xp { ZERO };
+    for (unsigned short d = 0; d < D; ++d) {
+      xp[d] = 2.0;
+    }
+    dist(xp, vp);
+    if (not Kokkos::isfinite(vp[0]) or not Kokkos::isfinite(vp[1]) or
+        not Kokkos::isfinite(vp[2])) {
+      raise::KernelError(HERE, "Non-finite velocity generated");
+    }
+    if constexpr (S == SimEngine::SRPIC) {
+      const auto gamma = math::sqrt(ONE + SQR(vp[0]) + SQR(vp[1]) + SQR(vp[2]));
+      if (gamma < 10 or gamma > 1000) {
+        raise::KernelError(HERE, "Gamma out of bounds");
+      }
+    } else {
+      vec_t<Dim::_3D> vup { ZERO };
+      metric.template transform<Idx::D, Idx::U>(xp, vp, vup);
+      const auto gamma = math::sqrt(
+        ONE + vup[0] * vp[0] + vup[1] * vp[1] + vup[2] * vp[2]);
+      if (gamma < 10 or gamma > 1000) {
+        raise::KernelError(HERE, "Gamma out of bounds");
+      }
+    }
+  }
+
+private:
+  M        metric;
+  EnrgDist dist;
+};
+
+template <SimEngine::type S, typename M>
+void testEnergyDist(const std::vector<std::size_t>&      res,
+                    const boundaries_t<real_t>&          ext,
+                    const std::map<std::string, real_t>& params = {}) {
+  raise::ErrorIf(res.size() != M::Dim, "res.size() != M::Dim", HERE);
+
+  boundaries_t<real_t> extent;
+  if constexpr (M::CoordType == Coord::Cart) {
+    extent = ext;
+  } else {
+    if constexpr (M::Dim == Dim::_2D) {
+      extent = {
+        ext[0],
+        {ZERO, constant::PI}
+      };
+    } else if constexpr (M::Dim == Dim::_3D) {
+      extent = {
+        ext[0],
+        {ZERO,     constant::PI},
+        {ZERO, constant::TWO_PI}
+      };
+    }
+  }
+  raise::ErrorIf(extent.size() != M::Dim, "extent.size() != M::Dim", HERE);
+
+  M metric { res, extent, params };
+
+  random_number_pool_t pool { constant::RandomSeed };
+  Powerlaw<S, M>       plaw { metric,
+                        pool,
+                        static_cast<real_t>(10),
+                        static_cast<real_t>(1000),
+                        static_cast<real_t>(-2.5) };
+  Kokkos::parallel_for("Powerlaw", 100, Caller<Powerlaw<S, M>, S, M>(metric, plaw));
+}
+
+auto main(int argc, char* argv[]) -> int {
+  Kokkos::initialize(argc, argv);
+
+  try {
+    using namespace ntt;
+    testEnergyDist<SimEngine::SRPIC, Minkowski<Dim::_1D>>(
+      {
+        10
+    },
+      { { 0.0, 55.0 } });
+
+    testEnergyDist<SimEngine::SRPIC, Minkowski<Dim::_2D>>(
+      {
+        10,
+        10
+    },
+      { { 0.0, 55.0 }, { 0.0, 55.0 } });
+
+    testEnergyDist<SimEngine::SRPIC, Minkowski<Dim::_3D>>(
+      {
+        10,
+        10,
+        10
+    },
+      { { 0.0, 55.0 }, { 0.0, 55.0 }, { 0.0, 55.0 } });
+
+    testEnergyDist<SimEngine::SRPIC, Spherical<Dim::_2D>>(
+      {
+        10,
+        10
+    },
+      { { 1.0, 100.0 } });
+
+    testEnergyDist<SimEngine::SRPIC, QSpherical<Dim::_2D>>(
+      {
+        10,
+        10
+    },
+      { { 1.0, 100.0 } },
+      { { "r0", 0.0 }, { "h", 0.25 } });
+
+    testEnergyDist<SimEngine::GRPIC, KerrSchild<Dim::_2D>>(
+      {
+        10,
+        10
+    },
+      { { 1.0, 100.0 } },
+      { { "a", 0.9 } });
+
+    testEnergyDist<SimEngine::GRPIC, QKerrSchild<Dim::_2D>>(
+      {
+        10,
+        10
+    },
+      { { 1.0, 100.0 } },
+      { { "r0", 0.0 }, { "h", 0.25 }, { "a", 0.9 } });
+
+    testEnergyDist<SimEngine::GRPIC, KerrSchild0<Dim::_2D>>(
+      {
+        10,
+        10
+    },
+      { { 1.0, 100.0 } },
+      { { "a", 0.9 } });
+
+  } catch (std::exception& e) {
+    std::cerr << e.what() << std::endl;
+    Kokkos::finalize();
+    return 1;
+  }
+  Kokkos::finalize();
+  return 0;
+}
diff --git a/src/framework/domain/output.cpp b/src/framework/domain/output.cpp
index cabb37093..be154ce16 100644
--- a/src/framework/domain/output.cpp
+++ b/src/framework/domain/output.cpp
@@ -26,6 +26,10 @@
 #include <Kokkos_ScatterView.hpp>
 #include <Kokkos_StdAlgorithms.hpp>
 
+#if defined(MPI_ENABLED)
+  #include <mpi.h>
+#endif // MPI_ENABLED
+
 #include <algorithm>
 #include <iterator>
 #include <vector>
@@ -458,34 +462,55 @@ namespace ntt {
                       ((D == Dim::_2D) and (M::CoordType != Coord::Cart))) {
           buff_x3 = array_t<real_t*> { "x3", nout };
         }
-        // clang-format off
-        Kokkos::parallel_for(
-          "PrtlToPhys",
-          nout,
-          kernel::PrtlToPhys_kernel<S, M>(prtl_stride,
-                                          buff_x1, buff_x2, buff_x3,
-                                          buff_ux1, buff_ux2, buff_ux3,
-                                          buff_wei,
-                                          species.i1, species.i2, species.i3,
-                                          species.dx1, species.dx2, species.dx3,
-                                          species.ux1, species.ux2, species.ux3,
-                                          species.phi, species.weight,
-                                          local_domain->mesh.metric));
-        // clang-format on
-        g_writer.writeParticleQuantity(buff_wei, prtl.name("W", 0));
-        g_writer.writeParticleQuantity(buff_ux1, prtl.name("U", 1));
-        g_writer.writeParticleQuantity(buff_ux2, prtl.name("U", 2));
-        g_writer.writeParticleQuantity(buff_ux3, prtl.name("U", 3));
+        if (nout > 0) {
+          // clang-format off
+          Kokkos::parallel_for(
+            "PrtlToPhys",
+            nout,
+            kernel::PrtlToPhys_kernel<S, M>(prtl_stride,
+                                            buff_x1, buff_x2, buff_x3,
+                                            buff_ux1, buff_ux2, buff_ux3,
+                                            buff_wei,
+                                            species.i1, species.i2, species.i3,
+                                            species.dx1, species.dx2, species.dx3,
+                                            species.ux1, species.ux2, species.ux3,
+                                            species.phi, species.weight,
+                                            local_domain->mesh.metric));
+          // clang-format on
+        }
+        std::size_t offset   = 0;
+        std::size_t glob_tot = nout;
+#if defined(MPI_ENABLED)
+        auto glob_nout = std::vector<std::size_t>(g_ndomains);
+        MPI_Allgather(&nout,
+                      1,
+                      mpi::get_type<std::size_t>(),
+                      glob_nout.data(),
+                      1,
+                      mpi::get_type<std::size_t>(),
+                      MPI_COMM_WORLD);
+        glob_tot = 0;
+        for (auto r = 0; r < g_mpi_size; ++r) {
+          if (r < g_mpi_rank) {
+            offset += glob_nout[r];
+          }
+          glob_tot += glob_nout[r];
+        }
+#endif // MPI_ENABLED
+        g_writer.writeParticleQuantity(buff_wei, glob_tot, offset, prtl.name("W", 0));
+        g_writer.writeParticleQuantity(buff_ux1, glob_tot, offset, prtl.name("U", 1));
+        g_writer.writeParticleQuantity(buff_ux2, glob_tot, offset, prtl.name("U", 2));
+        g_writer.writeParticleQuantity(buff_ux3, glob_tot, offset, prtl.name("U", 3));
         if constexpr (M::Dim == Dim::_1D or M::Dim == Dim::_2D or
                       M::Dim == Dim::_3D) {
-          g_writer.writeParticleQuantity(buff_x1, prtl.name("X", 1));
+          g_writer.writeParticleQuantity(buff_x1, glob_tot, offset, prtl.name("X", 1));
         }
         if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) {
-          g_writer.writeParticleQuantity(buff_x2, prtl.name("X", 2));
+          g_writer.writeParticleQuantity(buff_x2, glob_tot, offset, prtl.name("X", 2));
         }
         if constexpr (M::Dim == Dim::_3D or
                       ((D == Dim::_2D) and (M::CoordType != Coord::Cart))) {
-          g_writer.writeParticleQuantity(buff_x3, prtl.name("X", 3));
+          g_writer.writeParticleQuantity(buff_x3, glob_tot, offset, prtl.name("X", 3));
         }
       }
     } // end shouldWrite("particles", step, time)
diff --git a/src/framework/parameters.cpp b/src/framework/parameters.cpp
index 6e581a6b7..20c7e83d8 100644
--- a/src/framework/parameters.cpp
+++ b/src/framework/parameters.cpp
@@ -408,11 +408,16 @@ namespace ntt {
         toml::find_or(raw_data, "output", "fields", "stride", defaults::output::flds_stride));
 
     // particles
-    const auto prtl_out = toml::find_or(raw_data,
-                                        "output",
-                                        "particles",
-                                        "species",
-                                        std::vector<unsigned short> {});
+    auto prtl_out = toml::find_or(raw_data,
+                                  "output",
+                                  "particles",
+                                  "species",
+                                  std::vector<unsigned short> {});
+    if (prtl_out.size() == 0) {
+      for (unsigned short i = 0; i < species.size(); ++i) {
+        prtl_out.push_back(i + 1);
+      }
+    }
     set("output.particles.species", prtl_out);
     set("output.particles.stride",
         toml::find_or(raw_data,
diff --git a/src/output/writer.cpp b/src/output/writer.cpp
index 182705efa..91bf596a8 100644
--- a/src/output/writer.cpp
+++ b/src/output/writer.cpp
@@ -130,17 +130,20 @@ namespace out {
     for (const auto& prtl : m_prtl_writers) {
       for (auto d { 0u }; d < dim; ++d) {
         m_io.DefineVariable<real_t>(prtl.name("X", d + 1),
-                                    {},
-                                    {},
+                                    { adios2::UnknownDim },
+                                    { adios2::UnknownDim },
                                     { adios2::UnknownDim });
       }
       for (auto d { 0u }; d < Dim::_3D; ++d) {
         m_io.DefineVariable<real_t>(prtl.name("U", d + 1),
-                                    {},
-                                    {},
+                                    { adios2::UnknownDim },
+                                    { adios2::UnknownDim },
                                     { adios2::UnknownDim });
       }
-      m_io.DefineVariable<real_t>(prtl.name("W", 0), {}, {}, { adios2::UnknownDim });
+      m_io.DefineVariable<real_t>(prtl.name("W", 0),
+                                  { adios2::UnknownDim },
+                                  { adios2::UnknownDim },
+                                  { adios2::UnknownDim });
     }
   }
 
@@ -216,9 +219,13 @@ namespace out {
   }
 
   void Writer::writeParticleQuantity(const array_t<real_t*>& array,
+                                     std::size_t             glob_total,
+                                     std::size_t             loc_offset,
                                      const std::string&      varname) {
     auto var = m_io.InquireVariable<real_t>(varname);
-    var.SetSelection(adios2::Box<adios2::Dims>({}, { array.extent(0) }));
+    var.SetShape({ glob_total });
+    var.SetSelection(
+      adios2::Box<adios2::Dims>({ loc_offset }, { array.extent(0) }));
     auto array_h = Kokkos::create_mirror_view(array);
     Kokkos::deep_copy(array_h, array);
     m_writer.Put<real_t>(var, array_h);
diff --git a/src/output/writer.h b/src/output/writer.h
index d6a089095..517a1655d 100644
--- a/src/output/writer.h
+++ b/src/output/writer.h
@@ -114,7 +114,10 @@ namespace out {
                     const ndfield_t<D, N>&,
                     const std::vector<std::size_t>&);
 
-    void writeParticleQuantity(const array_t<real_t*>&, const std::string&);
+    void writeParticleQuantity(const array_t<real_t*>&,
+                               std::size_t,
+                               std::size_t,
+                               const std::string&);
     void writeSpectrum(const array_t<real_t*>&, const std::string&);
     void writeSpectrumBins(const array_t<real_t*>&, const std::string&);