diff --git a/src/Evolution/Executables/RadiationTransport/CMakeLists.txt b/src/Evolution/Executables/RadiationTransport/CMakeLists.txt
index 9a5cb7938cb7..060886a54225 100644
--- a/src/Evolution/Executables/RadiationTransport/CMakeLists.txt
+++ b/src/Evolution/Executables/RadiationTransport/CMakeLists.txt
@@ -2,3 +2,4 @@
 # See LICENSE.txt for details.
 
 add_subdirectory(M1Grey)
+add_subdirectory(MonteCarlo)
diff --git a/src/Evolution/Executables/RadiationTransport/MonteCarlo/CMakeLists.txt b/src/Evolution/Executables/RadiationTransport/MonteCarlo/CMakeLists.txt
new file mode 100644
index 000000000000..8680009bc4e2
--- /dev/null
+++ b/src/Evolution/Executables/RadiationTransport/MonteCarlo/CMakeLists.txt
@@ -0,0 +1,45 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+set(EXECUTABLE EvolveMonteCarlo)
+
+add_spectre_executable(
+  ${EXECUTABLE}
+  EXCLUDE_FROM_ALL
+  EvolveMonteCarlo.cpp
+  )
+
+target_link_libraries(
+  ${EXECUTABLE}
+  PRIVATE
+  Actions
+  Charmxx::main
+  CoordinateMaps
+  DataStructures
+  DgSubcell
+  DiscontinuousGalerkin
+  DomainCreators
+  Events
+  EventsAndDenseTriggers
+  EventsAndTriggers
+  Evolution
+  GeneralRelativity
+  GeneralRelativitySolutions
+  GrMhdAnalyticData
+  GrMhdSolutions
+  H5
+  Hydro
+  Informer
+  LinearOperators
+  MathFunctions
+  MonteCarlo
+  Observer
+  Options
+  Parallel
+  PhaseControl
+  Serialization
+  Time
+  Utilities
+  ValenciaDivClean
+  )
+
diff --git a/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.cpp b/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.cpp
new file mode 100644
index 000000000000..db4a08f04d8b
--- /dev/null
+++ b/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.cpp
@@ -0,0 +1,33 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.hpp"
+
+#include <vector>
+
+#include "Domain/Creators/RegisterDerivedWithCharm.hpp"
+#include "Domain/Creators/TimeDependence/RegisterDerivedWithCharm.hpp"
+#include "Domain/FunctionsOfTime/RegisterDerivedWithCharm.hpp"
+#include "Parallel/CharmMain.tpp"
+#include "PointwiseFunctions/Hydro/EquationsOfState/RegisterDerivedWithCharm.hpp"
+#include "Utilities/Serialization/RegisterDerivedClassesWithCharm.hpp"
+
+void register_neutrino_tables() {
+  register_classes_with_charm(
+      tmpl::list<Particles::MonteCarlo::NeutrinoInteractionTable<2, 2>,
+                 Particles::MonteCarlo::NeutrinoInteractionTable<2, 3>,
+                 Particles::MonteCarlo::NeutrinoInteractionTable<4, 3>,
+                 Particles::MonteCarlo::NeutrinoInteractionTable<16, 3>>{});
+}
+
+extern "C" void CkRegisterMainModule() {
+  Parallel::charmxx::register_main_module<EvolutionMetavars<4, 3>>();
+  Parallel::charmxx::register_init_node_and_proc(
+      {&domain::creators::register_derived_with_charm,
+       &domain::creators::time_dependence::register_derived_with_charm,
+       &domain::FunctionsOfTime::register_derived_with_charm,
+       &EquationsOfState::register_derived_with_charm,
+       &register_factory_classes_with_charm<EvolutionMetavars<4, 3>>,
+       &register_neutrino_tables},
+      {});
+}
diff --git a/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.hpp b/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.hpp
new file mode 100644
index 000000000000..e1bdde9784fd
--- /dev/null
+++ b/src/Evolution/Executables/RadiationTransport/MonteCarlo/EvolveMonteCarlo.hpp
@@ -0,0 +1,254 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstdint>
+#include <vector>
+
+#include "Domain/Creators/Factory3D.hpp"
+#include "Domain/Tags.hpp"
+#include "Evolution/Actions/RunEventsAndDenseTriggers.hpp"
+#include "Evolution/Actions/RunEventsAndTriggers.hpp"
+#include "Evolution/ComputeTags.hpp"
+#include "Evolution/DgSubcell/Actions/Initialize.hpp"
+#include "Evolution/DgSubcell/BackgroundGrVars.hpp"
+#include "Evolution/DiscontinuousGalerkin/Actions/ApplyBoundaryCorrections.hpp"
+#include "Evolution/DiscontinuousGalerkin/Actions/ComputeTimeDerivative.hpp"
+#include "Evolution/DiscontinuousGalerkin/Actions/VolumeTermsImpl.tpp"
+#include "Evolution/DiscontinuousGalerkin/BackgroundGrVars.hpp"
+#include "Evolution/DiscontinuousGalerkin/DgElementArray.hpp"
+#include "Evolution/DiscontinuousGalerkin/Initialization/Mortars.hpp"
+#include "Evolution/DiscontinuousGalerkin/Initialization/QuadratureTag.hpp"
+#include "Evolution/DiscontinuousGalerkin/Limiters/Minmod.hpp"
+#include "Evolution/DiscontinuousGalerkin/Limiters/Tags.hpp"
+#include "Evolution/Initialization/ConservativeSystem.hpp"
+#include "Evolution/Initialization/DgDomain.hpp"
+#include "Evolution/Initialization/Evolution.hpp"
+#include "Evolution/Initialization/Limiter.hpp"
+#include "Evolution/Initialization/SetVariables.hpp"
+#include "Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp"
+#include "Evolution/Particles/MonteCarlo/Actions/TimeStepActions.hpp"
+#include "Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp"
+#include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationTags.hpp"
+#include "Evolution/Particles/MonteCarlo/System.hpp"
+#include "Evolution/Particles/MonteCarlo/Tags.hpp"
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/AllSolutions.hpp"
+#include "Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/SwapGrTags.hpp"
+#include "IO/Observer/Actions/RegisterEvents.hpp"
+#include "IO/Observer/Helpers.hpp"
+#include "IO/Observer/ObserverComponent.hpp"
+#include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
+#include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
+#include "Options/Protocols/FactoryCreation.hpp"
+#include "Options/String.hpp"
+#include "Parallel/Local.hpp"
+#include "Parallel/Phase.hpp"
+#include "Parallel/PhaseControl/CheckpointAndExitAfterWallclock.hpp"
+#include "Parallel/PhaseControl/ExecutePhaseChange.hpp"
+#include "Parallel/PhaseControl/Factory.hpp"
+#include "Parallel/PhaseControl/VisitAndReturn.hpp"
+#include "Parallel/PhaseDependentActionList.hpp"
+#include "Parallel/Protocols/RegistrationMetavariables.hpp"
+#include "ParallelAlgorithms/Actions/AddComputeTags.hpp"
+#include "ParallelAlgorithms/Actions/AddSimpleTags.hpp"
+#include "ParallelAlgorithms/Actions/InitializeItems.hpp"
+#include "ParallelAlgorithms/Actions/LimiterActions.hpp"
+#include "ParallelAlgorithms/Actions/MutateApply.hpp"
+#include "ParallelAlgorithms/Actions/TerminatePhase.hpp"
+#include "ParallelAlgorithms/Events/Factory.hpp"
+#include "ParallelAlgorithms/EventsAndDenseTriggers/DenseTrigger.hpp"
+#include "ParallelAlgorithms/EventsAndDenseTriggers/DenseTriggers/Factory.hpp"
+#include "ParallelAlgorithms/EventsAndTriggers/Completion.hpp"
+#include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
+#include "ParallelAlgorithms/EventsAndTriggers/EventsAndTriggers.hpp"
+#include "ParallelAlgorithms/EventsAndTriggers/LogicalTriggers.hpp"
+#include "ParallelAlgorithms/EventsAndTriggers/Trigger.hpp"
+#include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
+#include "PointwiseFunctions/AnalyticData/GrMhd/InitialMagneticFields/InitialMagneticField.hpp"
+#include "PointwiseFunctions/AnalyticData/Tags.hpp"
+#include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
+#include "PointwiseFunctions/AnalyticSolutions/RadiationTransport/M1Grey/ConstantM1.hpp"
+#include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
+#include "PointwiseFunctions/Hydro/LowerSpatialFourVelocity.hpp"
+#include "PointwiseFunctions/Hydro/Tags.hpp"
+#include "PointwiseFunctions/InitialDataUtilities/Tags/InitialData.hpp"
+#include "Time/Actions/AdvanceTime.hpp"
+#include "Time/Actions/CleanHistory.hpp"
+#include "Time/Actions/RecordTimeStepperData.hpp"
+#include "Time/Actions/SelfStartActions.hpp"
+#include "Time/Actions/UpdateU.hpp"
+#include "Time/ChangeSlabSize/Action.hpp"
+#include "Time/StepChoosers/Factory.hpp"
+#include "Time/StepChoosers/StepChooser.hpp"
+#include "Time/Tags/Time.hpp"
+#include "Time/Tags/TimeStepId.hpp"
+#include "Time/TimeSequence.hpp"
+#include "Time/TimeSteppers/Factory.hpp"
+#include "Time/TimeSteppers/LtsTimeStepper.hpp"
+#include "Time/TimeSteppers/TimeStepper.hpp"
+#include "Time/Triggers/TimeTriggers.hpp"
+#include "Utilities/Functional.hpp"
+#include "Utilities/ProtocolHelpers.hpp"
+#include "Utilities/TMPL.hpp"
+
+/// \cond
+namespace Frame {
+struct Inertial;
+}  // namespace Frame
+namespace PUP {
+class er;
+}  // namespace PUP
+namespace Parallel {
+template <typename Metavariables>
+class CProxy_GlobalCache;
+}  // namespace Parallel
+/// \endcond
+
+// NEED:
+// Initial data
+
+template <size_t EnergyBins, size_t NeutrinoSpecies>
+struct EvolutionMetavars {
+  static constexpr size_t volume_dim = 3;
+
+  using system = Particles::MonteCarlo::System;
+  using temporal_id = Tags::TimeStepId;
+  using TimeStepperBase = TimeStepper;
+  static constexpr bool use_dg_subcell = true;
+
+  using initial_data_list =
+      grmhd::ValenciaDivClean::InitialData::initial_data_list;
+  using equation_of_state_tag = hydro::Tags::GrmhdEquationOfState;
+
+  struct SubcellOptions {
+    using evolved_vars_tags = typename system::variables_tag::tags_list;
+
+    static constexpr bool subcell_enabled = use_dg_subcell;
+    static constexpr bool subcell_enabled_at_external_boundary = true;
+  };
+
+  using observe_fields =
+      tmpl::list<domain::Tags::Coordinates<volume_dim, Frame::Grid>,
+                 domain::Tags::Coordinates<volume_dim, Frame::Inertial>>;
+  using non_tensor_compute_tags =
+      tmpl::list<::Events::Tags::ObserverMeshCompute<volume_dim>,
+                 ::Events::Tags::ObserverDetInvJacobianCompute<
+                     Frame::ElementLogical, Frame::Inertial>>;
+
+  using analytic_variables_tags = typename system::variables_tag::tags_list;
+  using analytic_compute = evolution::Tags::AnalyticSolutionsCompute<
+      volume_dim, analytic_variables_tags, false, initial_data_list>;
+
+  struct factory_creation
+      : tt::ConformsTo<Options::protocols::FactoryCreation> {
+    using factory_classes = tmpl::map<
+        tmpl::pair<DenseTrigger, DenseTriggers::standard_dense_triggers>,
+        tmpl::pair<DomainCreator<volume_dim>, domain_creators<volume_dim>>,
+        tmpl::pair<Event, tmpl::flatten<tmpl::list<Events::Completion>>>,
+        tmpl::pair<evolution::initial_data::InitialData, initial_data_list>,
+        tmpl::pair<
+            grmhd::AnalyticData::InitialMagneticFields::InitialMagneticField,
+            grmhd::AnalyticData::InitialMagneticFields::
+                initial_magnetic_fields>,
+        tmpl::pair<PhaseChange,
+                   tmpl::list<PhaseControl::CheckpointAndExitAfterWallclock>>,
+        tmpl::pair<StepChooser<StepChooserUse::Slab>,
+                   StepChoosers::standard_slab_choosers<system, false, false>>,
+        tmpl::pair<TimeSequence<double>,
+                   TimeSequences::all_time_sequences<double>>,
+        tmpl::pair<TimeSequence<std::uint64_t>,
+                   TimeSequences::all_time_sequences<std::uint64_t>>,
+        tmpl::pair<TimeStepper, TimeSteppers::time_steppers>,
+        tmpl::pair<Trigger, Triggers::time_triggers>>;
+  };
+
+  using observed_reduction_data_tags =
+      observers::collect_reduction_data_tags<tmpl::flatten<tmpl::list<
+          tmpl::at<typename factory_creation::factory_classes, Event>>>>;
+
+  using dg_registration_list =
+      tmpl::list<observers::Actions::RegisterEventsWithObservers>;
+
+  using initialization_actions = tmpl::list<
+      Initialization::Actions::InitializeItems<
+          Initialization::TimeStepping<EvolutionMetavars, TimeStepperBase>,
+          evolution::dg::Initialization::Domain<volume_dim>  //,
+          >,
+      Initialization::Actions::AddSimpleTags<
+          evolution::dg::BackgroundGrVars<system, EvolutionMetavars, true>>,
+      evolution::dg::subcell::Actions::SetSubcellGrid<volume_dim, system,
+                                                      false>,
+      Initialization::Actions::AddSimpleTags<
+          evolution::dg::subcell::BackgroundGrVars<system, EvolutionMetavars,
+                                                   true, false>>,
+      Actions::MutateApply<grmhd::ValenciaDivClean::subcell::SwapGrTags>,
+      Initialization::Actions::InitializeMCTags<system, EnergyBins,
+                                                NeutrinoSpecies>,
+      Initialization::Actions::AddComputeTags<tmpl::list<
+          hydro::Tags::LowerSpatialFourVelocityCompute,
+          Particles::MonteCarlo::InverseJacobianInertialToFluidCompute,
+          domain::Tags::JacobianCompute<4, Frame::Inertial, Frame::Fluid>>>,
+      evolution::Actions::InitializeRunEventsAndDenseTriggers,
+      Parallel::Actions::TerminatePhase>;
+
+  using dg_element_array = DgElementArray<
+      EvolutionMetavars,
+      tmpl::list<
+          Parallel::PhaseActions<Parallel::Phase::Initialization,
+                                 initialization_actions>,
+          Parallel::PhaseActions<
+              Parallel::Phase::Evolve,
+              tmpl::list<
+                  Actions::AdvanceTime,
+                  evolution::Actions::RunEventsAndTriggers<false>,
+                  Actions::AdvanceTime,
+                  evolution::Actions::RunEventsAndTriggers<false>,
+                  Actions::AdvanceTime,
+                  evolution::Actions::RunEventsAndTriggers<false>,
+                  Particles::MonteCarlo::Actions::SendDataForMcCommunication<
+                      volume_dim,
+                      // No local time stepping
+                      false, Particles::MonteCarlo::CommunicationStep::PreStep>,
+                  Particles::MonteCarlo::Actions::ReceiveDataForMcCommunication<
+                      volume_dim,
+                      Particles::MonteCarlo::CommunicationStep::PreStep>,
+                  Particles::MonteCarlo::Actions::TakeTimeStep<EnergyBins,
+                                                               NeutrinoSpecies>,
+                  Particles::MonteCarlo::Actions::SendDataForMcCommunication<
+                      volume_dim,
+                      // No local time stepping
+                      false,
+                      Particles::MonteCarlo::CommunicationStep::PostStep>,
+                  Particles::MonteCarlo::Actions::ReceiveDataForMcCommunication<
+                      volume_dim,
+                      Particles::MonteCarlo::CommunicationStep::PostStep>,
+                  PhaseControl::Actions::ExecutePhaseChange>>>>;
+
+  struct registration
+      : tt::ConformsTo<Parallel::protocols::RegistrationMetavariables> {
+    using element_registrars =
+        tmpl::map<tmpl::pair<dg_element_array, dg_registration_list>>;
+  };
+
+  using component_list =
+      tmpl::list<observers::Observer<EvolutionMetavars>,
+                 observers::ObserverWriter<EvolutionMetavars>,
+                 dg_element_array>;
+
+  using const_global_cache_tags = tmpl::list<
+      equation_of_state_tag, evolution::initial_data::Tags::InitialData,
+      Particles::MonteCarlo::Tags::InteractionRatesTable<EnergyBins,
+                                                         NeutrinoSpecies>>;
+
+  static constexpr Options::String help{
+      "Evolve Monte Carlo transport (without coupling to hydro).\n\n"};
+
+  static constexpr std::array<Parallel::Phase, 5> default_phase_order{
+      {Parallel::Phase::Initialization,
+       Parallel::Phase::InitializeTimeStepperHistory, Parallel::Phase::Register,
+       Parallel::Phase::Evolve, Parallel::Phase::Exit}};
+
+  // NOLINTNEXTLINE(google-runtime-references)
+  void pup(PUP::er& /*p*/) {}
+};
diff --git a/src/Evolution/Initialization/SetVariables.hpp b/src/Evolution/Initialization/SetVariables.hpp
index f8a28f5930c1..53eef9d8ee31 100644
--- a/src/Evolution/Initialization/SetVariables.hpp
+++ b/src/Evolution/Initialization/SetVariables.hpp
@@ -155,15 +155,18 @@ struct SetVariables {
       using variables_tag = typename system::variables_tag;
 
       // Set initial data from analytic solution
-      using Vars = typename variables_tag::type;
-      db::mutate<variables_tag>(
-          [&initial_time, &inertial_coords,
-           &solution_or_data](const gsl::not_null<Vars*> vars) {
-            vars->assign_subset(evolution::Initialization::initial_data(
-                solution_or_data, inertial_coords, initial_time,
-                typename Vars::tags_list{}));
-          },
-          box);
+      if constexpr (not std::is_same_v<typename variables_tag::tags_list,
+                                       tmpl::list<>>) {
+        using Vars = typename variables_tag::type;
+        db::mutate<variables_tag>(
+            [&initial_time, &inertial_coords,
+             &solution_or_data](const gsl::not_null<Vars*> vars) {
+              vars->assign_subset(evolution::Initialization::initial_data(
+                  solution_or_data, inertial_coords, initial_time,
+                  typename Vars::tags_list{}));
+            },
+            box);
+      }
     }
   }
 };
diff --git a/src/Evolution/Particles/MonteCarlo/Actions/CMakeLists.txt b/src/Evolution/Particles/MonteCarlo/Actions/CMakeLists.txt
index 2ab80257813e..d9e1b37ef3f3 100644
--- a/src/Evolution/Particles/MonteCarlo/Actions/CMakeLists.txt
+++ b/src/Evolution/Particles/MonteCarlo/Actions/CMakeLists.txt
@@ -5,5 +5,6 @@ spectre_target_headers(
   ${LIBRARY}
   INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
   HEADERS
+  InitializeMonteCarlo.hpp
   TimeStepActions.hpp
   )
diff --git a/src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp b/src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp
new file mode 100644
index 000000000000..d36e4b2dda05
--- /dev/null
+++ b/src/Evolution/Particles/MonteCarlo/Actions/InitializeMonteCarlo.hpp
@@ -0,0 +1,167 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <random>
+#include <tuple>
+
+#include "DataStructures/DataBox/DataBox.hpp"
+#include "Domain/Structure/Element.hpp"
+#include "Domain/Tags.hpp"
+#include "Evolution/DgSubcell/Tags/Coordinates.hpp"
+#include "Evolution/DgSubcell/Tags/Mesh.hpp"
+#include "Evolution/Initialization/InitialData.hpp"
+#include "Evolution/Particles/MonteCarlo/GhostZoneCommunicationTags.hpp"
+#include "Evolution/Particles/MonteCarlo/MortarData.hpp"
+#include "Evolution/Particles/MonteCarlo/Packet.hpp"
+#include "Evolution/Particles/MonteCarlo/Tags.hpp"
+#include "Parallel/AlgorithmExecution.hpp"
+#include "Parallel/GlobalCache.hpp"
+#include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
+#include "PointwiseFunctions/Hydro/Tags.hpp"
+#include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
+#include "Utilities/TMPL.hpp"
+
+/// \cond
+namespace evolution::initial_data::Tags {
+struct InitialData;
+}  // namespace evolution::initial_data::Tags
+
+namespace tuples {
+template <typename...>
+class TaggedTuple;
+}  // namespace tuples
+
+namespace Parallel {
+template <typename Metavariables>
+class GlobalCache;
+}  // namespace Parallel
+/// \endcond
+
+namespace Initialization::Actions {
+
+/// \ingroup InitializationGroup
+/// \brief Allocate variables needed for evolution of Monte Carlo transport
+///
+/// Uses:
+/// - evolution::dg::subcell::Tags::Mesh<dim>
+/// - evolution::dg::subcell::Tags::Coordinates<dim, Frame::Inertial>
+/// - evolution::initial_data::Tags::InitialData
+///
+/// DataBox changes:
+/// - Adds:
+///   * Particles::MonteCarlo::Tags::PacketsOnElement
+///   * Particles::MonteCarlo::Tags::RandomNumberGenerator
+///   * Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission<
+///                                  NeutrinoSpecies>
+///   * Particles::MonteCarlo::Tags::CellLightCrossingTime<DataVector>
+///   * Background hydro variables
+///   * Particles::MonteCarlo::Tags::MortarDataTag<dim>
+///   * Particles::MonteCarlo::Tags::McGhostZoneDataTag<dim>
+///
+/// - Removes: nothing
+/// - Modifies: nothing
+template <typename System, size_t EnergyBins, size_t NeutrinoSpecies>
+struct InitializeMCTags {
+ public:
+  using hydro_variables_tag = typename System::hydro_variables_tag;
+
+  static constexpr size_t dim = System::volume_dim;
+  using simple_tags =
+      tmpl::list<Particles::MonteCarlo::Tags::PacketsOnElement,
+                 Particles::MonteCarlo::Tags::RandomNumberGenerator,
+                 Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission<
+                     NeutrinoSpecies>,
+                 Particles::MonteCarlo::Tags::CellLightCrossingTime<DataVector>,
+                 hydro_variables_tag,
+                 Particles::MonteCarlo::Tags::MortarDataTag<dim>,
+                 Particles::MonteCarlo::Tags::McGhostZoneDataTag<dim>>;
+
+  using compute_tags = tmpl::list<>;
+
+  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
+            typename ArrayIndex, typename ActionList,
+            typename ParallelComponent>
+  static Parallel::iterable_action_return_t apply(
+      db::DataBox<DbTagsList>& box,
+      const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
+      const Parallel::GlobalCache<Metavariables>& /*cache*/,
+      const ArrayIndex& /*array_index*/, ActionList /*meta*/,
+      const ParallelComponent* const /*meta*/) {
+    const size_t num_grid_points =
+        db::get<evolution::dg::subcell::Tags::Mesh<dim>>(box)
+            .number_of_grid_points();
+    using derived_classes =
+        tmpl::at<typename Metavariables::factory_creation::factory_classes,
+                 evolution::initial_data::InitialData>;
+    using HydroVars = typename hydro_variables_tag::type;
+    call_with_dynamic_type<void, derived_classes>(
+        &db::get<evolution::initial_data::Tags::InitialData>(box),
+        [&box](const auto* const data_or_solution) {
+          static constexpr size_t dim = System::volume_dim;
+          const double initial_time = db::get<::Tags::Time>(box);
+          const size_t num_grid_points =
+              db::get<evolution::dg::subcell::Tags::Mesh<dim>>(box)
+                  .number_of_grid_points();
+          const auto& inertial_coords = db::get<
+              evolution::dg::subcell::Tags::Coordinates<dim, Frame::Inertial>>(
+              box);
+          // Get hydro variables
+          HydroVars hydro_variables{num_grid_points};
+          hydro_variables.assign_subset(evolution::Initialization::initial_data(
+              *data_or_solution, inertial_coords, initial_time,
+              typename hydro_variables_tag::tags_list{}));
+          Initialization::mutate_assign<tmpl::list<hydro_variables_tag>>(
+              make_not_null(&box), std::move(hydro_variables));
+        });
+
+    typename Particles::MonteCarlo::Tags::PacketsOnElement::type all_packets;
+    Initialization::mutate_assign<
+        tmpl::list<Particles::MonteCarlo::Tags::PacketsOnElement>>(
+        make_not_null(&box), std::move(all_packets));
+
+    const unsigned long seed =
+        std::random_device{}();  // static_cast<unsigned long>(time(NULL));
+    typename Particles::MonteCarlo::Tags::RandomNumberGenerator::type rng(seed);
+
+    Initialization::mutate_assign<
+        tmpl::list<Particles::MonteCarlo::Tags::RandomNumberGenerator>>(
+        make_not_null(&box), std::move(rng));
+
+    // Currently hard-code energy at emission; should be set by option
+    typename Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission<
+        NeutrinoSpecies>::type packet_energy_at_emission =
+        make_with_value<std::array<DataVector, NeutrinoSpecies>>(
+            DataVector{num_grid_points}, 1.e-12);
+    Initialization::mutate_assign<
+        tmpl::list<Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission<
+            NeutrinoSpecies>>>(make_not_null(&box),
+                               std::move(packet_energy_at_emission));
+
+    typename Particles::MonteCarlo::Tags::CellLightCrossingTime<
+        DataVector>::type cell_light_crossing_time(num_grid_points, 1.0);
+    Initialization::mutate_assign<tmpl::list<
+        Particles::MonteCarlo::Tags::CellLightCrossingTime<DataVector>>>(
+        make_not_null(&box), std::move(cell_light_crossing_time));
+
+    // Initialize empty mortar data; do we need more at initialization stage?
+    using MortarData =
+        typename Particles::MonteCarlo::Tags::MortarDataTag<dim>::type;
+    MortarData mortar_data;
+    Initialization::mutate_assign<
+        tmpl::list<Particles::MonteCarlo::Tags::MortarDataTag<dim>>>(
+        make_not_null(&box), std::move(mortar_data));
+
+    using GhostZoneData =
+        typename Particles::MonteCarlo::Tags::McGhostZoneDataTag<dim>::type;
+    GhostZoneData ghost_zone_data{};
+    Initialization::mutate_assign<
+        tmpl::list<Particles::MonteCarlo::Tags::McGhostZoneDataTag<dim>>>(
+        make_not_null(&box), std::move(ghost_zone_data));
+
+    return {Parallel::AlgorithmExecution::Continue, std::nullopt};
+  }
+};
+
+}  // namespace Initialization::Actions
diff --git a/src/Evolution/Particles/MonteCarlo/Actions/TimeStepActions.hpp b/src/Evolution/Particles/MonteCarlo/Actions/TimeStepActions.hpp
index f31cca715fee..42bf52adaad2 100644
--- a/src/Evolution/Particles/MonteCarlo/Actions/TimeStepActions.hpp
+++ b/src/Evolution/Particles/MonteCarlo/Actions/TimeStepActions.hpp
@@ -41,7 +41,8 @@ struct TimeStepMutator {
   using return_tags =
       tmpl::list<Particles::MonteCarlo::Tags::PacketsOnElement,
                  Particles::MonteCarlo::Tags::RandomNumberGenerator,
-                 Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission<3>>;
+                 Particles::MonteCarlo::Tags::DesiredPacketEnergyAtEmission<
+                     NeutrinoSpecies>>;
   // To do : check carefully DG vs Subcell quantities... everything should
   // be on the Subcell grid!
   using argument_tags = tmpl::list<
@@ -60,11 +61,11 @@ struct TimeStepMutator {
                     Frame::Inertial>,
       ::Tags::deriv<gr::Tags::Shift<DataVector, Dim>, tmpl::size_t<Dim>,
                     Frame::Inertial>,
-      ::Tags::deriv<gr::Tags::InverseSpatialMetric<DataVector, Dim>,
-                    tmpl::size_t<Dim>, Frame::Inertial>,
+      ::Tags::deriv<gr::Tags::SpatialMetric<DataVector, Dim>, tmpl::size_t<Dim>,
+                    Frame::Inertial>,
       gr::Tags::SpatialMetric<DataVector, Dim, Frame::Inertial>,
       gr::Tags::InverseSpatialMetric<DataVector, Dim, Frame::Inertial>,
-      gr::Tags::DetSpatialMetric<DataVector>,
+      gr::Tags::SqrtDetSpatialMetric<DataVector>,
       Particles::MonteCarlo::Tags::CellLightCrossingTime<DataVector>,
       evolution::dg::subcell::Tags::Mesh<Dim>,
       evolution::dg::subcell::Tags::Coordinates<Dim, Frame::ElementLogical>,
@@ -96,10 +97,10 @@ struct TimeStepMutator {
 
       const tnsr::i<DataVector, Dim, Frame::Inertial>& d_lapse,
       const tnsr::iJ<DataVector, Dim, Frame::Inertial>& d_shift,
-      const tnsr::iJJ<DataVector, Dim, Frame::Inertial>& d_inv_spatial_metric,
+      const tnsr::ijj<DataVector, Dim, Frame::Inertial>& d_spatial_metric,
       const tnsr::ii<DataVector, Dim, Frame::Inertial>& spatial_metric,
       const tnsr::II<DataVector, Dim, Frame::Inertial>& inv_spatial_metric,
-      const Scalar<DataVector>& determinant_spatial_metric,
+      const Scalar<DataVector>& sqrt_determinant_spatial_metric,
       const Scalar<DataVector>& cell_light_crossing_time, const Mesh<Dim>& mesh,
       const tnsr::I<DataVector, Dim, Frame::ElementLogical>& mesh_coordinates,
       const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
@@ -131,6 +132,22 @@ struct TimeStepMutator {
     const DirectionalIdMap<Dim, std::optional<DataVector>>&
         cell_light_crossing_time_ghost = mortar_data.cell_light_crossing_time;
 
+    tnsr::iJJ<DataVector, 3, Frame::Inertial> d_inv_spatial_metric =
+        make_with_value<tnsr::iJJ<DataVector, 3, Frame::Inertial>>(lapse, 0.0);
+    for (size_t i = 0; i < 3; i++) {
+      for (size_t j = i; j < 3; j++) {
+        for (size_t k = 0; k < 3; k++) {
+          for (size_t l = 0; l < 3; l++) {
+            for (size_t m = 0; m < 3; m++) {
+              d_inv_spatial_metric.get(k, i, j) -=
+                  inv_spatial_metric.get(i, l) * inv_spatial_metric.get(j, m) *
+                  d_spatial_metric.get(k, l, m);
+            }
+          }
+        }
+      }
+    }
+
     TemplatedLocalFunctions<EnergyBins, NeutrinoSpecies> templated_functions;
     templated_functions.take_time_step_on_element(
         packets, random_number_generator, single_packet_energy, start_time,
@@ -138,7 +155,7 @@ struct TimeStepMutator {
         rest_mass_density, temperature, lorentz_factor,
         lower_spatial_four_velocity, lapse, shift, d_lapse, d_shift,
         d_inv_spatial_metric, spatial_metric, inv_spatial_metric,
-        determinant_spatial_metric, cell_light_crossing_time, mesh,
+        sqrt_determinant_spatial_metric, cell_light_crossing_time, mesh,
         mesh_coordinates, num_ghost_zones, mesh_velocity,
         inverse_jacobian_logical_to_inertial, det_jacobian_logical_to_inertial,
         inertial_to_fluid_jacobian, inertial_to_fluid_inverse_jacobian,
diff --git a/src/Evolution/Particles/MonteCarlo/CMakeLists.txt b/src/Evolution/Particles/MonteCarlo/CMakeLists.txt
index 0f8c3f106098..daca3ec8939d 100644
--- a/src/Evolution/Particles/MonteCarlo/CMakeLists.txt
+++ b/src/Evolution/Particles/MonteCarlo/CMakeLists.txt
@@ -38,6 +38,7 @@ spectre_target_headers(
   NeutrinoInteractionTable.hpp
   Packet.hpp
   Scattering.hpp
+  System.hpp
   Tags.hpp
   TakeTimeStep.tpp
   TemplatedLocalFunctions.hpp
diff --git a/src/Evolution/Particles/MonteCarlo/CellVolume.cpp b/src/Evolution/Particles/MonteCarlo/CellVolume.cpp
index 89f2a71267ee..4f416722fff6 100644
--- a/src/Evolution/Particles/MonteCarlo/CellVolume.cpp
+++ b/src/Evolution/Particles/MonteCarlo/CellVolume.cpp
@@ -12,13 +12,13 @@ namespace Particles::MonteCarlo {
 void cell_proper_four_volume_finite_difference(
     const gsl::not_null<Scalar<DataVector>*> cell_proper_four_volume,
     const Scalar<DataVector>& lapse,
-    const Scalar<DataVector>& determinant_spatial_metric,
+    const Scalar<DataVector>& sqrt_determinant_spatial_metric,
     const double time_step, const Mesh<3>& mesh,
     const Scalar<DataVector>& det_jacobian_logical_to_inertial) {
   const double cell_logical_volume =
     8.0 / static_cast<double>(mesh.number_of_grid_points());
   cell_proper_four_volume->get() =
-      get(lapse) * get(determinant_spatial_metric) * time_step *
+      get(lapse) * get(sqrt_determinant_spatial_metric) * time_step *
       cell_logical_volume * get(det_jacobian_logical_to_inertial);
 }
 
diff --git a/src/Evolution/Particles/MonteCarlo/CellVolume.hpp b/src/Evolution/Particles/MonteCarlo/CellVolume.hpp
index f813836ef9a9..52db8d797b98 100644
--- a/src/Evolution/Particles/MonteCarlo/CellVolume.hpp
+++ b/src/Evolution/Particles/MonteCarlo/CellVolume.hpp
@@ -25,12 +25,11 @@ namespace Particles::MonteCarlo {
 /// in inertial coordinates, hence the need for
 /// det_jacobian_logical_to_inertial
 void cell_proper_four_volume_finite_difference(
-  gsl::not_null<Scalar<DataVector>* > cell_proper_four_volume,
-  const Scalar<DataVector>& lapse,
-  const Scalar<DataVector>& determinant_spatial_metric,
-  double time_step,
-  const Mesh<3>& mesh,
-  const Scalar<DataVector>& det_jacobian_logical_to_inertial);
+    gsl::not_null<Scalar<DataVector>*> cell_proper_four_volume,
+    const Scalar<DataVector>& lapse,
+    const Scalar<DataVector>& sqrt_determinant_spatial_metric, double time_step,
+    const Mesh<3>& mesh,
+    const Scalar<DataVector>& det_jacobian_logical_to_inertial);
 
 /// 3-volume of a cell in inertial coordinate. Note that this is
 /// the coordinate volume, not the proper volume. This quantity
diff --git a/src/Evolution/Particles/MonteCarlo/EvolvePacketsInElement.tpp b/src/Evolution/Particles/MonteCarlo/EvolvePacketsInElement.tpp
index 3abfe9d8d072..effa473212c1 100644
--- a/src/Evolution/Particles/MonteCarlo/EvolvePacketsInElement.tpp
+++ b/src/Evolution/Particles/MonteCarlo/EvolvePacketsInElement.tpp
@@ -162,6 +162,22 @@ void TemplatedLocalFunctions<EnergyBins, NeutrinoSpecies>::evolve_packets(
     initial_time = packet.time;
     dt_end_step = final_time - initial_time;
 
+    // Find closest grid point to packet at current time, using
+    // extents for live points only.
+    {
+      std::array<size_t, 3> closest_point_index_3d{0, 0, 0};
+      for (size_t d = 0; d < 3; d++) {
+        gsl::at(closest_point_index_3d, d) =
+          std::floor((packet.coordinates[d] - gsl::at(bottom_coord_mesh, d)) /
+                         gsl::at(dx_mesh, d) +
+                     0.5);
+      }
+      packet.index_of_closest_grid_point =
+          closest_point_index_3d[0] +
+          extents[0] * (closest_point_index_3d[1] +
+                        extents[1] * closest_point_index_3d[2]);
+    }
+
     // Get quantities that we do NOT update if the packet
     // changes cell.
     // local_idx is the index on the mesh without ghost zones
diff --git a/src/Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp b/src/Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp
index 2f9422a5848a..0ced96571ea5 100644
--- a/src/Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp
+++ b/src/Evolution/Particles/MonteCarlo/GhostZoneCommunication.hpp
@@ -136,6 +136,16 @@ struct GhostDataMcPackets {
         const size_t d = max_distance_direction.value().dimension();
         const Side side = max_distance_direction.value().side();
         packet.coordinates.get(d) += (side == Side::Lower) ? (2.0) : (-2.0);
+        // Corner/edge treatment; move packet in a live point for now
+        // Definitely needs improvement...
+        for (size_t dd = 0; dd < Dim; dd++) {
+          if (dd != d && packet.coordinates.get(dd) < -1.0) {
+            packet.coordinates.get(dd) = -2.0 - packet.coordinates.get(dd);
+          }
+          if (dd != d && packet.coordinates.get(dd) > 1.0) {
+            packet.coordinates.get(dd) = 2.0 - packet.coordinates.get(dd);
+          }
+        }
         if (output.contains(max_distance_direction.value())) {
           output[max_distance_direction.value()].push_back(packet);
         }
@@ -311,8 +321,8 @@ struct ReceiveDataForMcCommunication {
                 const DataVector& received_data_direction =
                     received_data[directional_element_id]
                         .ghost_zone_hydro_variables;
-                REQUIRE(received_data[directional_element_id]
-                            .packets_entering_this_element == std::nullopt);
+                // REQUIRE(received_data[directional_element_id]
+                //             .packets_entering_this_element == std::nullopt);
                 if (mortar_data->rest_mass_density[directional_element_id] ==
                     std::nullopt) {
                   continue;
@@ -368,8 +378,8 @@ struct ReceiveDataForMcCommunication {
                         received_data[directional_element_id]
                             .packets_entering_this_element;
                 // Temporary: currently no data for coupling to the fluid
-                REQUIRE(received_data[directional_element_id]
-                            .ghost_zone_hydro_variables.size() == 0);
+                // REQUIRE(received_data[directional_element_id]
+                //            .ghost_zone_hydro_variables.size() == 0);
                 if (received_data_packets == std::nullopt) {
                   continue;
                 } else {
diff --git a/src/Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.cpp b/src/Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.cpp
index cbbbf8b9e251..e5894408c2bf 100644
--- a/src/Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.cpp
+++ b/src/Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.cpp
@@ -306,6 +306,7 @@ void NeutrinoInteractionTable<EnergyBins, NeutrinoSpecies>::
   const size_t n_ye_points = table_electron_fraction.size();
   double temperature_correction_factor = 0.0;
 
+  const double dye = table_electron_fraction[1] - table_electron_fraction[0];
   double ye = 0.0;
   double log_rho = 0.0;
   double log_temp = 0.0;
@@ -323,8 +324,8 @@ void NeutrinoInteractionTable<EnergyBins, NeutrinoSpecies>::
             ? 1.0
             : get(temperature)[p] / minimum_temperature;
 
-    ye = std::clamp(ye, table_electron_fraction[0],
-                    table_electron_fraction[n_ye_points - 1]);
+    ye = std::clamp(ye, table_electron_fraction[0] + 1.e-8 * dye,
+                    table_electron_fraction[n_ye_points - 1] - 1.e-8 * dye);
     log_rho = std::clamp(log_rho, table_log_density[0],
                          table_log_density[n_rho_points - 1]);
     log_temp = std::clamp(log_temp, table_log_temperature[0],
diff --git a/src/Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.hpp b/src/Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.hpp
index 69c057c841ae..db027ad2567f 100644
--- a/src/Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.hpp
+++ b/src/Evolution/Particles/MonteCarlo/NeutrinoInteractionTable.hpp
@@ -11,6 +11,7 @@
 #include "DataStructures/BoostMultiArray.hpp"
 #include "DataStructures/Tensor/TypeAliases.hpp"
 #include "NumericalAlgorithms/Interpolation/MultiLinearSpanInterpolation.hpp"
+#include "Options/String.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/Serialization/CharmPupable.hpp"
 
@@ -28,6 +29,19 @@ class NeutrinoInteractionTable : public PUP::able {
   /// Read table from disk and stores interaction rates.
   explicit NeutrinoInteractionTable(const std::string& filename);
 
+  static constexpr Options::String help = {
+      "Tabulated neutrino-matter interactions in NuLib format.\n"
+      "Emissivity, absorption and scattering opacities are    \n"
+      "tabulated as a function of density, eletron fraction   \n"
+      "and temperature."};
+
+  struct TableFilename {
+    using type = std::string;
+    static constexpr Options::String help{"File name of the NuLib table"};
+  };
+
+  using options = tmpl::list<TableFilename>;
+
   /// Explicit instantiation from table values, for tests
   NeutrinoInteractionTable(
       std::vector<double> table_data_,
diff --git a/src/Evolution/Particles/MonteCarlo/System.hpp b/src/Evolution/Particles/MonteCarlo/System.hpp
new file mode 100644
index 000000000000..c638ed2b45b8
--- /dev/null
+++ b/src/Evolution/Particles/MonteCarlo/System.hpp
@@ -0,0 +1,41 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstddef>
+
+#include "DataStructures/VariablesTag.hpp"
+#include "Evolution/Particles/MonteCarlo/Tags.hpp"
+#include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
+#include "PointwiseFunctions/Hydro/Tags.hpp"
+#include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
+#include "Utilities/TMPL.hpp"
+
+namespace Particles::MonteCarlo {
+
+struct System {
+  static constexpr bool is_in_flux_conservative_form = false;
+  static constexpr bool has_primitive_and_conservative_vars = false;
+  static constexpr size_t volume_dim = 3;
+  // The EoS is used within the MC code itself, even if we provide
+  // the fluid variables as a background.
+  static constexpr size_t thermodynamic_dim = 3;
+
+  using mc_variables_tag = ::Tags::Variables<
+      tmpl::list<Particles::MonteCarlo::Tags::PacketsOnElement>>;
+  using variables_tag = ::Tags::Variables<tmpl::list<>>;  // mc_variables_tag;
+  using flux_variables = tmpl::list<>;
+  using gradient_variables = tmpl::list<>;
+  // GR tags needed for background metric
+  using spacetime_variables_tag =
+      ::Tags::Variables<gr::tags_for_hydro<volume_dim, DataVector>>;
+  using flux_spacetime_variables_tag = ::Tags::Variables<tmpl::list<>>;
+  // Hydro tags needed for background metric
+  using hydro_variables_tag = ::Tags::Variables<hydro::grmhd_tags<DataVector>>;
+  using primitive_variables_tag = hydro_variables_tag;
+
+  using inverse_spatial_metric_tag =
+      gr::Tags::InverseSpatialMetric<DataVector, 3>;
+};
+}  // namespace Particles::MonteCarlo
diff --git a/src/Evolution/Particles/MonteCarlo/Tags.hpp b/src/Evolution/Particles/MonteCarlo/Tags.hpp
index ed0fb3080df8..e0274fc1ae89 100644
--- a/src/Evolution/Particles/MonteCarlo/Tags.hpp
+++ b/src/Evolution/Particles/MonteCarlo/Tags.hpp
@@ -51,6 +51,18 @@ template <size_t EnergyBins, size_t NeutrinoSpecies>
 struct InteractionRatesTable : db::SimpleTag {
   using type =
       std::unique_ptr<NeutrinoInteractionTable<EnergyBins, NeutrinoSpecies>>;
+  static constexpr bool pass_metavariables = false;
+  using option_tags =
+      typename NeutrinoInteractionTable<EnergyBins, NeutrinoSpecies>::options;
+  static type create_from_options(const std::string filename) {
+    std::unique_ptr<Particles::MonteCarlo::NeutrinoInteractionTable<
+        EnergyBins, NeutrinoSpecies>>
+        interaction_table_ptr =
+            std::make_unique<Particles::MonteCarlo::NeutrinoInteractionTable<
+                EnergyBins, NeutrinoSpecies>>(filename);
+    return interaction_table_ptr;
+    ;
+  }
 };
 
 }  // namespace Particles::MonteCarlo::Tags
diff --git a/src/Evolution/Particles/MonteCarlo/TakeTimeStep.tpp b/src/Evolution/Particles/MonteCarlo/TakeTimeStep.tpp
index 00838bf3f48b..06c715cbd48f 100644
--- a/src/Evolution/Particles/MonteCarlo/TakeTimeStep.tpp
+++ b/src/Evolution/Particles/MonteCarlo/TakeTimeStep.tpp
@@ -98,7 +98,7 @@ void TemplatedLocalFunctions<EnergyBins, NeutrinoSpecies>::
         const tnsr::iJJ<DataVector, 3, Frame::Inertial>& d_inv_spatial_metric,
         const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
         const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
-        const Scalar<DataVector>& determinant_spatial_metric,
+        const Scalar<DataVector>& sqrt_determinant_spatial_metric,
         const Scalar<DataVector>& cell_light_crossing_time,
 
         const Mesh<3>& mesh,
@@ -137,8 +137,8 @@ void TemplatedLocalFunctions<EnergyBins, NeutrinoSpecies>::
   Scalar<DataVector> cell_inertial_three_volume =
       make_with_value<Scalar<DataVector>>(lapse, 0.0);
   cell_proper_four_volume_finite_difference(
-      &cell_proper_four_volume, lapse, determinant_spatial_metric, time_step,
-      mesh, det_jacobian_logical_to_inertial);
+      &cell_proper_four_volume, lapse, sqrt_determinant_spatial_metric,
+      time_step, mesh, det_jacobian_logical_to_inertial);
   cell_inertial_coordinate_three_volume_finite_difference(
       &cell_inertial_three_volume, mesh, det_jacobian_logical_to_inertial);
 
diff --git a/src/Evolution/Particles/MonteCarlo/TemplatedLocalFunctions.hpp b/src/Evolution/Particles/MonteCarlo/TemplatedLocalFunctions.hpp
index b88054158b99..e62fbc5b67a8 100644
--- a/src/Evolution/Particles/MonteCarlo/TemplatedLocalFunctions.hpp
+++ b/src/Evolution/Particles/MonteCarlo/TemplatedLocalFunctions.hpp
@@ -85,7 +85,7 @@ struct TemplatedLocalFunctions {
       const tnsr::iJJ<DataVector, 3, Frame::Inertial>& d_inv_spatial_metric,
       const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
       const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
-      const Scalar<DataVector>& determinant_spatial_metric,
+      const Scalar<DataVector>& sqrt_determinant_spatial_metric,
       const Scalar<DataVector>& cell_light_crossing_time, const Mesh<3>& mesh,
       const tnsr::I<DataVector, 3, Frame::ElementLogical>& mesh_coordinates,
       size_t num_ghost_zones,
diff --git a/src/NumericalAlgorithms/Interpolation/MultiLinearSpanInterpolation.hpp b/src/NumericalAlgorithms/Interpolation/MultiLinearSpanInterpolation.hpp
index ad29e279b095..7de83e9076a3 100644
--- a/src/NumericalAlgorithms/Interpolation/MultiLinearSpanInterpolation.hpp
+++ b/src/NumericalAlgorithms/Interpolation/MultiLinearSpanInterpolation.hpp
@@ -234,19 +234,32 @@ size_t MultiLinearSpanInterpolation<
              << which_dimension << "\ncurrent_index: " << current_index
              << "\nnumber of points: " << number_of_points_[which_dimension]
              << "\ntarget point: " << target_points
-             << "\nrelative coordinate: " << relative_coordinate);
+             << "\nrelative coordinate: " << relative_coordinate
+             << "\nspacing: " << inverse_spacing_[which_dimension]);
+
+  // Do not trigger errors for roundoff errors in index calculation
+  if (UNLIKELY(current_index + 1 == number_of_points_[which_dimension])) {
+    if (relative_coordinate * inverse_spacing_[which_dimension] -
+            static_cast<double>(current_index) <
+        1.e-13) {
+      current_index -= 1;
+    }
+  }
 
   // We are exceeding the table bounds:
   // Use linear extrapolation based of the highest
   // two points in the table
 
-  ASSERT(allow_extrapolation_abov_data_[which_dimension] or
-             UNLIKELY(current_index + 1 < number_of_points_[which_dimension]),
-         "Interpolation exceeds upper table bounds.\nwhich_dimension: "
-             << which_dimension << "\ncurrent_index: " << current_index
-             << "\nnumber of points: " << number_of_points_[which_dimension]
-             << "\ntarget point: " << target_points
-             << "\nrelative coordinate: " << relative_coordinate);
+  ASSERT(
+      allow_extrapolation_abov_data_[which_dimension] or
+          UNLIKELY(current_index + 1 < number_of_points_[which_dimension]),
+      "Interpolation exceeds upper table bounds.\nwhich_dimension: "
+          << which_dimension << "\ncurrent_index: " << current_index
+          << "\nnumber of points: " << number_of_points_[which_dimension]
+          << "\ntarget point: " << target_points
+          << "\nrelative coordinate: " << relative_coordinate << "\nposition: "
+          << relative_coordinate * inverse_spacing_[which_dimension] -
+                 static_cast<double>(number_of_points_[which_dimension] - 1));
 
   // Enforce index ranges
   current_index = std::min(number_of_points_[which_dimension] - 2,
diff --git a/tests/Unit/Evolution/Particles/MonteCarlo/Test_CellVolume.cpp b/tests/Unit/Evolution/Particles/MonteCarlo/Test_CellVolume.cpp
index 7e440d0a20b6..61ac032d2d0b 100644
--- a/tests/Unit/Evolution/Particles/MonteCarlo/Test_CellVolume.cpp
+++ b/tests/Unit/Evolution/Particles/MonteCarlo/Test_CellVolume.cpp
@@ -15,18 +15,20 @@ SPECTRE_TEST_CASE("Unit.Evolution.Particles.MonteCarloCellVolume",
                      Spectral::Quadrature::CellCentered);
 
   const size_t dv_size = 27;
-  DataVector zero_dv(dv_size, 0.0);
-  Scalar<DataVector> lapse{DataVector(dv_size, 1.2)};
-  Scalar<DataVector> determinant_spatial_metric{DataVector(dv_size, 0.9)};
-  Scalar<DataVector> det_jacobian_logical_to_inertial{DataVector(dv_size, 1.1)};
+  const DataVector zero_dv(dv_size, 0.0);
+  const Scalar<DataVector> lapse{DataVector(dv_size, 1.2)};
+  const Scalar<DataVector> sqrt_determinant_spatial_metric
+    {DataVector(dv_size, 0.9)};
+  const Scalar<DataVector> det_jacobian_logical_to_inertial
+    {DataVector(dv_size, 1.1)};
   const double time_step = 0.6;
 
   Scalar<DataVector> cell_proper_four_volume{DataVector(dv_size, 0.0)};
   Scalar<DataVector> expected_cell_proper_four_volume{
       DataVector(dv_size, 8.0 / 27.0 * 1.2 * 0.9 * 0.6 * 1.1)};
   Particles::MonteCarlo::cell_proper_four_volume_finite_difference(
-      &cell_proper_four_volume, lapse, determinant_spatial_metric, time_step,
-      mesh, det_jacobian_logical_to_inertial);
+      &cell_proper_four_volume, lapse, sqrt_determinant_spatial_metric,
+      time_step, mesh, det_jacobian_logical_to_inertial);
 
   Scalar<DataVector> cell_inertial_three_volume{DataVector(dv_size, 0.0)};
   Scalar<DataVector> expected_cell_inertial_three_volume{
diff --git a/tests/Unit/Evolution/Particles/MonteCarlo/Test_CommunicationTags.cpp b/tests/Unit/Evolution/Particles/MonteCarlo/Test_CommunicationTags.cpp
index a91b6f4c8f5e..46aaa87fdf8a 100644
--- a/tests/Unit/Evolution/Particles/MonteCarlo/Test_CommunicationTags.cpp
+++ b/tests/Unit/Evolution/Particles/MonteCarlo/Test_CommunicationTags.cpp
@@ -524,6 +524,10 @@ void test_send_receive_actions() {
   // topological coordinate of their new element.
   packet_east.coordinates[0] -= 2.0;
   packet_south.coordinates[1] += 2.0;
+  // Current hack for edges/corners
+  if constexpr (Dim > 2){
+    packet_south.coordinates[2] = -2.0 - packet_south.coordinates[2];
+  }
   {
     const auto& east_data = ActionTesting::get_inbox_tag<
         comp, Particles::MonteCarlo::McGhostZoneDataInboxTag<
diff --git a/tests/Unit/Evolution/Particles/MonteCarlo/Test_TakeTimeStep.cpp b/tests/Unit/Evolution/Particles/MonteCarlo/Test_TakeTimeStep.cpp
index f15870ef41f6..15ebadbc052f 100644
--- a/tests/Unit/Evolution/Particles/MonteCarlo/Test_TakeTimeStep.cpp
+++ b/tests/Unit/Evolution/Particles/MonteCarlo/Test_TakeTimeStep.cpp
@@ -33,15 +33,15 @@ void test_flat_space_time_step() {
   const size_t num_ghost_zones = 1;
 
   const size_t dv_size = cube(size_1d);
-  DataVector zero_dv(dv_size, 0.0);
+  const DataVector zero_dv(dv_size, 0.0);
   const size_t dv_size_with_ghost = cube(size_1d + 2 * num_ghost_zones);
-  DataVector zero_dv_with_ghost(dv_size_with_ghost, 0.0);
+  const DataVector zero_dv_with_ghost(dv_size_with_ghost, 0.0);
   const size_t dv_size_in_ghost = square(size_1d) * num_ghost_zones;
-  DataVector zero_dv_in_ghost(dv_size_in_ghost, 0.0);
-  DataVector one_dv_in_ghost(dv_size_in_ghost, 1.0);
+  const DataVector zero_dv_in_ghost(dv_size_in_ghost, 0.0);
+  const DataVector one_dv_in_ghost(dv_size_in_ghost, 1.0);
 
   // Minkowski metric
-  Scalar<DataVector> lapse{DataVector(dv_size, 1.0)};
+  const Scalar<DataVector> lapse{DataVector(dv_size, 1.0)};
   tnsr::II<DataVector, 3, Frame::Inertial> inv_spatial_metric =
       make_with_value<tnsr::II<DataVector, 3, Frame::Inertial>>(lapse, 0.0);
   inv_spatial_metric.get(0, 0) = 1.0;
@@ -52,14 +52,14 @@ void test_flat_space_time_step() {
   spatial_metric.get(0, 0) = 1.0;
   spatial_metric.get(1, 1) = 1.0;
   spatial_metric.get(2, 2) = 1.0;
-  Scalar<DataVector> determinant_spatial_metric(dv_size, 1.0);
-  tnsr::I<DataVector, 3, Frame::Inertial> shift =
+  const Scalar<DataVector> sqrt_determinant_spatial_metric(dv_size, 1.0);
+  const tnsr::I<DataVector, 3, Frame::Inertial> shift =
       make_with_value<tnsr::I<DataVector, 3, Frame::Inertial>>(lapse, 0.0);
-  tnsr::i<DataVector, 3, Frame::Inertial> d_lapse =
+  const tnsr::i<DataVector, 3, Frame::Inertial> d_lapse =
       make_with_value<tnsr::i<DataVector, 3, Frame::Inertial>>(lapse, 0.0);
-  tnsr::iJ<DataVector, 3, Frame::Inertial> d_shift =
+  const tnsr::iJ<DataVector, 3, Frame::Inertial> d_shift =
       make_with_value<tnsr::iJ<DataVector, 3, Frame::Inertial>>(lapse, 0.0);
-  tnsr::iJJ<DataVector, 3, Frame::Inertial> d_inv_spatial_metric =
+  const tnsr::iJJ<DataVector, 3, Frame::Inertial> d_inv_spatial_metric =
       make_with_value<tnsr::iJJ<DataVector, 3, Frame::Inertial>>(lapse, 0.0);
 
   // Mesh velocity set to std::null for now
@@ -90,7 +90,7 @@ void test_flat_space_time_step() {
   inverse_jacobian_logical_to_inertial.get(0, 0) = 1.0;
   inverse_jacobian_logical_to_inertial.get(1, 1) = 1.0;
   inverse_jacobian_logical_to_inertial.get(2, 2) = 1.0;
-  Scalar<DataVector> det_jacobian_logical_to_inertial(dv_size, 1.0);
+  const Scalar<DataVector> det_jacobian_logical_to_inertial(dv_size, 1.0);
 
   // Coordinates
   tnsr::I<DataVector, 3, Frame::ElementLogical> mesh_coordinates =
@@ -238,7 +238,7 @@ void test_flat_space_time_step() {
         electron_fraction, baryon_density, temperature, lorentz_factor,
         lower_spatial_four_velocity, lapse, shift, d_lapse, d_shift,
         d_inv_spatial_metric, spatial_metric, inv_spatial_metric,
-        determinant_spatial_metric, cell_light_crossing_time, mesh,
+        sqrt_determinant_spatial_metric, cell_light_crossing_time, mesh,
         mesh_coordinates, num_ghost_zones, mesh_velocity,
         inverse_jacobian_logical_to_inertial, det_jacobian_logical_to_inertial,
         jacobian_inertial_to_fluid, inverse_jacobian_inertial_to_fluid,
diff --git a/tests/Unit/Evolution/Particles/MonteCarlo/Test_TimeStepAction.cpp b/tests/Unit/Evolution/Particles/MonteCarlo/Test_TimeStepAction.cpp
index 2f2aec53b353..d7dca4d480ee 100644
--- a/tests/Unit/Evolution/Particles/MonteCarlo/Test_TimeStepAction.cpp
+++ b/tests/Unit/Evolution/Particles/MonteCarlo/Test_TimeStepAction.cpp
@@ -88,11 +88,11 @@ struct component {
                   Frame::Inertial>,
       Tags::deriv<gr::Tags::Shift<DataVector, Dim>, tmpl::size_t<Dim>,
                   Frame::Inertial>,
-      Tags::deriv<gr::Tags::InverseSpatialMetric<DataVector, Dim>,
+      Tags::deriv<gr::Tags::SpatialMetric<DataVector, Dim>,
                   tmpl::size_t<Dim>, Frame::Inertial>,
       gr::Tags::SpatialMetric<DataVector, Dim, Frame::Inertial>,
       gr::Tags::InverseSpatialMetric<DataVector, Dim, Frame::Inertial>,
-      gr::Tags::DetSpatialMetric<DataVector>,
+      gr::Tags::SqrtDetSpatialMetric<DataVector>,
       evolution::dg::subcell::Tags::Coordinates<Dim, Frame::ElementLogical>,
       domain::Tags::MeshVelocity<Dim>,
       evolution::dg::subcell::fd::Tags::InverseJacobianLogicalToInertial<Dim>,
@@ -185,15 +185,15 @@ void test_advance_packets() {
   spatial_metric.get(0, 0) = 1.0;
   spatial_metric.get(1, 1) = 1.0;
   spatial_metric.get(2, 2) = 1.0;
-  Scalar<DataVector> determinant_spatial_metric(n_pts, 1.0);
+  Scalar<DataVector> sqrt_determinant_spatial_metric(n_pts, 1.0);
   tnsr::I<DataVector, 3, Frame::Inertial> shift =
       make_with_value<tnsr::I<DataVector, 3, Frame::Inertial>>(lapse, 0.0);
   tnsr::i<DataVector, 3, Frame::Inertial> d_lapse =
       make_with_value<tnsr::i<DataVector, 3, Frame::Inertial>>(lapse, 0.0);
   tnsr::iJ<DataVector, 3, Frame::Inertial> d_shift =
       make_with_value<tnsr::iJ<DataVector, 3, Frame::Inertial>>(lapse, 0.0);
-  tnsr::iJJ<DataVector, 3, Frame::Inertial> d_inv_spatial_metric =
-      make_with_value<tnsr::iJJ<DataVector, 3, Frame::Inertial>>(lapse, 0.0);
+  tnsr::ijj<DataVector, 3, Frame::Inertial> d_spatial_metric =
+      make_with_value<tnsr::ijj<DataVector, 3, Frame::Inertial>>(lapse, 0.0);
 
   // Fluid variables
   Scalar<DataVector> rest_mass_density(zero_dv);
@@ -313,10 +313,10 @@ void test_advance_packets() {
        shift,
        d_lapse,
        d_shift,
-       d_inv_spatial_metric,
+       d_spatial_metric,
        spatial_metric,
        inv_spatial_metric,
-       determinant_spatial_metric,
+       sqrt_determinant_spatial_metric,
        mesh_coordinates,
        mesh_velocity,
        inverse_jacobian_logical_to_inertial,