diff --git a/CI/physmon/reference/performance_ambi_ttbar.root b/CI/physmon/reference/performance_ambi_ttbar.root
index 447854c935b9..3113b5bf5597 100644
Binary files a/CI/physmon/reference/performance_ambi_ttbar.root and b/CI/physmon/reference/performance_ambi_ttbar.root differ
diff --git a/CI/physmon/reference/performance_amvf_gridseeder_ttbar_hist.root b/CI/physmon/reference/performance_amvf_gridseeder_ttbar_hist.root
index 22f53f42a82d..e2b42ad1e84e 100644
Binary files a/CI/physmon/reference/performance_amvf_gridseeder_ttbar_hist.root and b/CI/physmon/reference/performance_amvf_gridseeder_ttbar_hist.root differ
diff --git a/CI/physmon/reference/performance_amvf_ttbar_hist.root b/CI/physmon/reference/performance_amvf_ttbar_hist.root
index 8e78161a3628..dd394d86bc77 100644
Binary files a/CI/physmon/reference/performance_amvf_ttbar_hist.root and b/CI/physmon/reference/performance_amvf_ttbar_hist.root differ
diff --git a/CI/physmon/reference/performance_ckf_ttbar.root b/CI/physmon/reference/performance_ckf_ttbar.root
index 7a8bb020d163..4cf17f1902c8 100644
Binary files a/CI/physmon/reference/performance_ckf_ttbar.root and b/CI/physmon/reference/performance_ckf_ttbar.root differ
diff --git a/CI/physmon/reference/performance_gsf.root b/CI/physmon/reference/performance_gsf.root
index 725398e90a7c..faa6c70a7df5 100644
Binary files a/CI/physmon/reference/performance_gsf.root and b/CI/physmon/reference/performance_gsf.root differ
diff --git a/CI/physmon/reference/performance_seeding_ttbar.root b/CI/physmon/reference/performance_seeding_ttbar.root
index e99cd72f2f45..700dfed6e7d2 100644
Binary files a/CI/physmon/reference/performance_seeding_ttbar.root and b/CI/physmon/reference/performance_seeding_ttbar.root differ
diff --git a/CI/physmon/reference/tracksummary_ckf_ttbar_hist.root b/CI/physmon/reference/tracksummary_ckf_ttbar_hist.root
index 467ed85119dd..f71e399df040 100644
Binary files a/CI/physmon/reference/tracksummary_ckf_ttbar_hist.root and b/CI/physmon/reference/tracksummary_ckf_ttbar_hist.root differ
diff --git a/Core/include/Acts/Geometry/Layer.hpp b/Core/include/Acts/Geometry/Layer.hpp
index 2195020c78f1..494126ccaaa9 100644
--- a/Core/include/Acts/Geometry/Layer.hpp
+++ b/Core/include/Acts/Geometry/Layer.hpp
@@ -35,7 +35,7 @@ class VolumeBounds;
 class TrackingVolume;
 class ApproachDescriptor;
 class IMaterialDecorator;
-template <typename T>
+template <typename object_t>
 struct NavigationOptions;
 
 // Simple surface intersection
diff --git a/Core/include/Acts/Propagator/EigenStepper.hpp b/Core/include/Acts/Propagator/EigenStepper.hpp
index fc5cf28f64a3..c17e45e8bb61 100644
--- a/Core/include/Acts/Propagator/EigenStepper.hpp
+++ b/Core/include/Acts/Propagator/EigenStepper.hpp
@@ -162,8 +162,7 @@ class EigenStepper {
   };
 
   /// Constructor requires knowledge of the detector's magnetic field
-  EigenStepper(std::shared_ptr<const MagneticFieldProvider> bField,
-               double overstepLimit = 100 * UnitConstants::um);
+  EigenStepper(std::shared_ptr<const MagneticFieldProvider> bField);
 
   State makeState(std::reference_wrapper<const GeometryContext> gctx,
                   std::reference_wrapper<const MagneticFieldContext> mctx,
@@ -317,12 +316,6 @@ class EigenStepper {
     return state.stepSize.toString();
   }
 
-  /// Overstep limit
-  double overstepLimit(const State& /*state*/) const {
-    // A dynamic overstep limit could sit here
-    return -m_overstepLimit;
-  }
-
   /// Create and return the bound state at the current position
   ///
   /// @brief This transports (if necessary) the covariance
@@ -435,9 +428,6 @@ class EigenStepper {
  protected:
   /// Magnetic field inside of the detector
   std::shared_ptr<const MagneticFieldProvider> m_bField;
-
-  /// Overstep limit
-  double m_overstepLimit{};
 };
 
 template <typename navigator_t>
diff --git a/Core/include/Acts/Propagator/EigenStepper.ipp b/Core/include/Acts/Propagator/EigenStepper.ipp
index e41a3fb67c13..4ecad2c31ebb 100644
--- a/Core/include/Acts/Propagator/EigenStepper.ipp
+++ b/Core/include/Acts/Propagator/EigenStepper.ipp
@@ -15,8 +15,8 @@
 
 template <typename E, typename A>
 Acts::EigenStepper<E, A>::EigenStepper(
-    std::shared_ptr<const MagneticFieldProvider> bField, double overstepLimit)
-    : m_bField(std::move(bField)), m_overstepLimit(overstepLimit) {}
+    std::shared_ptr<const MagneticFieldProvider> bField)
+    : m_bField(std::move(bField)) {}
 
 template <typename E, typename A>
 auto Acts::EigenStepper<E, A>::makeState(
diff --git a/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp b/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp
index 30f60177e710..4637a4acb34c 100644
--- a/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp
+++ b/Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp
@@ -726,14 +726,6 @@ class MultiEigenStepperLoop
     return ss.str();
   }
 
-  /// Overstep limit
-  ///
-  /// @param state [in] The stepping state (thread-local cache)
-  double overstepLimit(const State& state) const {
-    // A dynamic overstep limit could sit here
-    return SingleStepper::overstepLimit(state.components.front().state);
-  }
-
   /// Create and return the bound state at the current position
   ///
   /// @brief This transports (if necessary) the covariance
diff --git a/Core/include/Acts/Propagator/Navigator.hpp b/Core/include/Acts/Propagator/Navigator.hpp
index c56d21248055..f01e33a60af5 100644
--- a/Core/include/Acts/Propagator/Navigator.hpp
+++ b/Core/include/Acts/Propagator/Navigator.hpp
@@ -9,7 +9,6 @@
 #pragma once
 
 #include "Acts/Definitions/Units.hpp"
-#include "Acts/Geometry/BoundarySurfaceT.hpp"
 #include "Acts/Geometry/GeometryIdentifier.hpp"
 #include "Acts/Geometry/Layer.hpp"
 #include "Acts/Geometry/TrackingGeometry.hpp"
@@ -19,8 +18,6 @@
 #include "Acts/Utilities/Logger.hpp"
 #include "Acts/Utilities/StringHelpers.hpp"
 
-#include <iomanip>
-#include <iterator>
 #include <sstream>
 #include <string>
 
@@ -31,7 +28,6 @@ namespace Acts {
 /// @brief struct for the Navigation options that are forwarded to
 ///        the geometry
 ///
-/// @tparam propagator_state_t Type of the object for navigation state
 /// @tparam object_t Type of the object for navigation to check against
 template <typename object_t>
 struct NavigationOptions {
@@ -46,9 +42,9 @@ struct NavigationOptions {
   /// always look for passive
   bool resolvePassive = false;
 
-  /// object to check against: at start
+  /// Hint for start object
   const object_t* startObject = nullptr;
-  /// object to check against: at end
+  /// Hint for end object
   const object_t* endObject = nullptr;
 
   /// External surface identifier for which the boundary check is ignored
@@ -856,7 +852,13 @@ class Navigator {
               state.geoContext, stepper.position(state.stepping),
               state.options.direction * stepper.direction(state.stepping),
               navOpts, logger());
-      // The number of boundary candidates
+      std::sort(state.navigation.navBoundaries.begin(),
+                state.navigation.navBoundaries.end(),
+                [](const auto& a, const auto& b) {
+                  return SurfaceIntersection::pathLengthOrder(a.first, b.first);
+                });
+
+      // Print boundary information
       if (logger().doPrint(Logging::VERBOSE)) {
         std::ostringstream os;
         os << state.navigation.navBoundaries.size();
@@ -866,6 +868,7 @@ class Navigator {
         }
         logger().log(Logging::VERBOSE, os.str());
       }
+
       // Set the begin index
       state.navigation.navBoundaryIndex = 0;
       if (!state.navigation.navBoundaries.empty()) {
@@ -1049,18 +1052,23 @@ class Navigator {
     state.navigation.navSurfaces = currentLayer->compatibleSurfaces(
         state.geoContext, stepper.position(state.stepping),
         state.options.direction * stepper.direction(state.stepping), navOpts);
-    // the number of layer candidates
-    if (!state.navigation.navSurfaces.empty()) {
-      if (logger().doPrint(Logging::VERBOSE)) {
-        std::ostringstream os;
-        os << state.navigation.navSurfaces.size();
-        os << " surface candidates found at path(s): ";
-        for (auto& sfc : state.navigation.navSurfaces) {
-          os << sfc.pathLength() << "  ";
-        }
-        logger().log(Logging::VERBOSE, os.str());
+    std::sort(state.navigation.navSurfaces.begin(),
+              state.navigation.navSurfaces.end(),
+              SurfaceIntersection::pathLengthOrder);
+
+    // Print surface information
+    if (logger().doPrint(Logging::VERBOSE)) {
+      std::ostringstream os;
+      os << state.navigation.navSurfaces.size();
+      os << " surface candidates found at path(s): ";
+      for (auto& sfc : state.navigation.navSurfaces) {
+        os << sfc.pathLength() << "  ";
       }
+      logger().log(Logging::VERBOSE, os.str());
+    }
 
+    // Surface candidates have been found
+    if (!state.navigation.navSurfaces.empty()) {
       // set the index
       state.navigation.navSurfaceIndex = 0;
       // The stepper updates the step size ( single / multi component)
@@ -1070,6 +1078,7 @@ class Navigator {
                                   << stepper.outputStepSize(state.stepping));
       return true;
     }
+
     state.navigation.navSurfaceIndex = state.navigation.navSurfaces.size();
     ACTS_VERBOSE(volInfo(state) << "No surface candidates found.");
     return false;
@@ -1104,25 +1113,32 @@ class Navigator {
     navOpts.nearLimit = state.options.surfaceTolerance;
     navOpts.farLimit =
         stepper.getStepSize(state.stepping, ConstrainedStep::aborter);
+
     // Request the compatible layers
     state.navigation.navLayers =
         state.navigation.currentVolume->compatibleLayers(
             state.geoContext, stepper.position(state.stepping),
             state.options.direction * stepper.direction(state.stepping),
             navOpts);
+    std::sort(state.navigation.navLayers.begin(),
+              state.navigation.navLayers.end(),
+              [](const auto& a, const auto& b) {
+                return SurfaceIntersection::pathLengthOrder(a.first, b.first);
+              });
+
+    // Print layer information
+    if (logger().doPrint(Logging::VERBOSE)) {
+      std::ostringstream os;
+      os << state.navigation.navLayers.size();
+      os << " layer candidates found at path(s): ";
+      for (auto& lc : state.navigation.navLayers) {
+        os << lc.first.pathLength() << "  ";
+      }
+      logger().log(Logging::VERBOSE, os.str());
+    }
 
     // Layer candidates have been found
     if (!state.navigation.navLayers.empty()) {
-      // Screen output where they are
-      if (logger().doPrint(Logging::VERBOSE)) {
-        std::ostringstream os;
-        os << state.navigation.navLayers.size();
-        os << " layer candidates found at path(s): ";
-        for (auto& lc : state.navigation.navLayers) {
-          os << lc.first.pathLength() << "  ";
-        }
-        logger().log(Logging::VERBOSE, os.str());
-      }
       // Set the index to the first
       state.navigation.navLayerIndex = 0;
       // Setting the step size towards first
diff --git a/Core/include/Acts/Propagator/StandardAborters.hpp b/Core/include/Acts/Propagator/StandardAborters.hpp
index b9af52d0b9d9..a8c07f073df0 100644
--- a/Core/include/Acts/Propagator/StandardAborters.hpp
+++ b/Core/include/Acts/Propagator/StandardAborters.hpp
@@ -72,10 +72,10 @@ struct SurfaceReached {
   /// Distance limit to discard intersections "behind us"
   /// @note this is only necessary because some surfaces have more than one
   ///       intersection
-  std::optional<double> overrideNearLimit;
+  double nearLimit = -100 * UnitConstants::um;
 
   SurfaceReached() = default;
-  SurfaceReached(double nearLimit) : overrideNearLimit(nearLimit) {}
+  SurfaceReached(double nLimit) : nearLimit(nLimit) {}
 
   /// boolean operator for abort condition without using the result
   ///
@@ -101,14 +101,11 @@ struct SurfaceReached {
       return true;
     }
 
-    // not blindly using the stepper overstep limit here because it does not
-    // always work for perigee surfaces.
+    // not using the stepper overstep limit here because it does not always work
+    // for perigee surfaces
     // note: the near limit is necessary for surfaces with more than one
-    // intersections in order to discard the ones which are behind us.
-    const double nearLimit =
-        overrideNearLimit.value_or(stepper.overstepLimit(state.stepping));
-    const double farLimit =
-        state.stepping.stepSize.value(ConstrainedStep::aborter);
+    // intersection in order to discard the ones which are behind us
+    const double farLimit = std::numeric_limits<double>::max();
     const double tolerance = state.options.surfaceTolerance;
 
     const auto sIntersection = surface->intersect(
diff --git a/Core/include/Acts/Propagator/StepperConcept.hpp b/Core/include/Acts/Propagator/StepperConcept.hpp
index c62f31274060..fdb43394dff8 100644
--- a/Core/include/Acts/Propagator/StepperConcept.hpp
+++ b/Core/include/Acts/Propagator/StepperConcept.hpp
@@ -44,7 +44,6 @@ METHOD_TRAIT(absolute_momentum_t, absoluteMomentum);
 METHOD_TRAIT(momentum_t, momentum);
 METHOD_TRAIT(charge_t, charge);
 METHOD_TRAIT(time_t, time);
-METHOD_TRAIT(overstep_t, overstepLimit);
 METHOD_TRAIT(bound_state_method_t, boundState);
 METHOD_TRAIT(curvilinear_state_method_t, curvilinearState);
 METHOD_TRAIT(update_t, update);
@@ -114,8 +113,6 @@ constexpr bool MultiStepperStateConcept= require<
         static_assert(charge_exists, "charge method not found");
         constexpr static bool time_exists = has_method<const S, double, time_t, const state&>;
         static_assert(time_exists, "time method not found");
-        constexpr static bool overstep_exists = has_method<const S, double, overstep_t, const state&>;
-        static_assert(overstep_exists, "overstepLimit method not found");
         constexpr static bool bound_state_method_exists= has_method<const S, Result<typename S::BoundState>, bound_state_method_t, state&, const Surface&, bool, const FreeToBoundCorrection&>;
         static_assert(bound_state_method_exists, "boundState method not found");
         constexpr static bool curvilinear_state_method_exists = has_method<const S, typename S::CurvilinearState, curvilinear_state_method_t, state&, bool>;
diff --git a/Core/include/Acts/Propagator/StraightLineStepper.hpp b/Core/include/Acts/Propagator/StraightLineStepper.hpp
index 4a8bec6396e7..231d7ba77885 100644
--- a/Core/include/Acts/Propagator/StraightLineStepper.hpp
+++ b/Core/include/Acts/Propagator/StraightLineStepper.hpp
@@ -218,11 +218,6 @@ class StraightLineStepper {
   /// @param state [in] The stepping state (thread-local cache)
   double time(const State& state) const { return state.pars[eFreeTime]; }
 
-  /// Overstep limit
-  double overstepLimit(const State& /*state*/) const {
-    return -m_overstepLimit;
-  }
-
   /// Update surface status
   ///
   /// This method intersects the provided surface and update the navigation
@@ -451,9 +446,6 @@ class StraightLineStepper {
     // return h
     return h;
   }
-
- private:
-  double m_overstepLimit = s_onSurfaceTolerance;
 };
 
 template <typename navigator_t>
diff --git a/Core/include/Acts/Propagator/detail/SteppingHelper.hpp b/Core/include/Acts/Propagator/detail/SteppingHelper.hpp
index 23a24c62191d..cd8f330f1efe 100644
--- a/Core/include/Acts/Propagator/detail/SteppingHelper.hpp
+++ b/Core/include/Acts/Propagator/detail/SteppingHelper.hpp
@@ -17,6 +17,8 @@
 #include "Acts/Utilities/Intersection.hpp"
 #include "Acts/Utilities/Logger.hpp"
 
+#include <limits>
+
 namespace Acts::detail {
 
 /// Update surface status - Single component
@@ -50,8 +52,7 @@ Acts::Intersection3D::Status updateSingleSurfaceStatus(
     return Intersection3D::Status::onSurface;
   }
 
-  // Path and overstep limit checking
-  const double nearLimit = stepper.overstepLimit(state);
+  const double nearLimit = std::numeric_limits<double>::lowest();
   const double farLimit = state.stepSize.value(ConstrainedStep::aborter);
 
   if (sIntersection && detail::checkIntersection(sIntersection.intersection(),
diff --git a/Core/include/Acts/Utilities/Intersection.hpp b/Core/include/Acts/Utilities/Intersection.hpp
index d76c50a78b84..f30da706c071 100644
--- a/Core/include/Acts/Utilities/Intersection.hpp
+++ b/Core/include/Acts/Utilities/Intersection.hpp
@@ -217,6 +217,10 @@ class ObjectMultiIntersection {
     return {m_intersections[index], m_object, index};
   }
 
+  constexpr const MultiIntersection3D& intersections() const {
+    return m_intersections;
+  }
+
   constexpr std::size_t size() const { return m_intersections.size(); }
 
   constexpr const object_t* object() const { return m_object; }
@@ -274,7 +278,7 @@ bool checkIntersection(const intersection_t& intersection, double nearLimit,
                << nearLimit << ", " << farLimit << ", " << distance);
 
   const bool coCriterion = distance > nearLimit;
-  const bool cpCriterion = std::abs(distance) < std::abs(farLimit) + tolerance;
+  const bool cpCriterion = distance < farLimit + tolerance;
 
   const bool accept = coCriterion && cpCriterion;
 
@@ -288,9 +292,9 @@ bool checkIntersection(const intersection_t& intersection, double nearLimit,
     }
     if (!cpCriterion) {
       ACTS_VERBOSE("- intersection path length "
-                   << std::abs(distance) << " is over the far limit "
-                   << (std::abs(farLimit) + tolerance)
-                   << " (including tolerance of " << tolerance << ")");
+                   << distance << " is over the far limit "
+                   << (farLimit + tolerance) << " (including tolerance of "
+                   << tolerance << ")");
     }
   }
 
diff --git a/Core/src/Geometry/Layer.cpp b/Core/src/Geometry/Layer.cpp
index deeeaed56775..d420f2919a68 100644
--- a/Core/src/Geometry/Layer.cpp
+++ b/Core/src/Geometry/Layer.cpp
@@ -119,11 +119,13 @@ Acts::Layer::compatibleSurfaces(
     return sIntersections;
   }
 
+  double nearLimit = options.nearLimit;
+  double farLimit = options.farLimit;
+
+  // TODO this looks like a major hack; is this really needed?
   // (0) End surface check
   // @todo: - we might be able to skip this by use of options.pathLimit
   // check if you have to stop at the endSurface
-  double nearLimit = options.nearLimit;
-  double farLimit = options.farLimit;
   if (options.endObject != nullptr) {
     // intersect the end surface
     // - it is the final one don't use the boundary check at all
@@ -150,6 +152,14 @@ Acts::Layer::compatibleSurfaces(
     farLimit = 1.5 * thickness() * pCorrection;
   }
 
+  auto isUnique = [&](const SurfaceIntersection& b) {
+    auto find_it = std::find_if(
+        sIntersections.begin(), sIntersections.end(), [&b](const auto& a) {
+          return a.object() == b.object() && a.index() == b.index();
+        });
+    return find_it == sIntersections.end();
+  };
+
   // lemma 0 : accept the surface
   auto acceptSurface = [&options](const Surface& sf,
                                   bool sensitive = false) -> bool {
@@ -168,8 +178,8 @@ Acts::Layer::compatibleSurfaces(
   // lemma 1 : check and fill the surface
   // [&sIntersections, &options, &parameters
   auto processSurface = [&](const Surface& sf, bool sensitive = false) {
-    // veto if it's start or end surface
-    if (options.startObject == &sf || options.endObject == &sf) {
+    // veto if it's start surface
+    if (options.startObject == &sf) {
       return;
     }
     // veto if it doesn't fit the prescription
@@ -183,14 +193,13 @@ Acts::Layer::compatibleSurfaces(
       boundaryCheck = false;
     }
     // the surface intersection
-    SurfaceMultiIntersection sfmi =
-        sf.intersect(gctx, position, direction, BoundaryCheck(boundaryCheck));
-    for (const auto& sfi : sfmi.split()) {
-      // check if intersection is valid and pathLimit has not been exceeded
-      if (sfi &&
-          detail::checkIntersection(sfi.intersection(), nearLimit, farLimit)) {
-        sIntersections.push_back(sfi);
-      }
+    SurfaceIntersection sfi =
+        sf.intersect(gctx, position, direction, BoundaryCheck(boundaryCheck))
+            .closest();
+    if (sfi &&
+        detail::checkIntersection(sfi.intersection(), nearLimit, farLimit) &&
+        isUnique(sfi)) {
+      sIntersections.push_back(sfi);
     }
   };
 
@@ -234,26 +243,6 @@ Acts::Layer::compatibleSurfaces(
   const Surface* layerSurface = &surfaceRepresentation();
   processSurface(*layerSurface);
 
-  // Sort by object address
-  std::sort(
-      sIntersections.begin(), sIntersections.end(),
-      [](const auto& a, const auto& b) { return a.object() < b.object(); });
-  // Now look for duplicates. As we just sorted by path length, duplicates
-  // should be subsequent
-  auto it = std::unique(
-      sIntersections.begin(), sIntersections.end(),
-      [](const SurfaceIntersection& a, const SurfaceIntersection& b) -> bool {
-        return a.object() == b.object();
-      });
-
-  // resize to remove all items that are past the unique range
-  sIntersections.resize(std::distance(sIntersections.begin(), it),
-                        SurfaceIntersection::invalid());
-
-  // sort according to the path length
-  std::sort(sIntersections.begin(), sIntersections.end(),
-            SurfaceIntersection::pathLengthOrder);
-
   return sIntersections;
 }
 
@@ -270,7 +259,7 @@ Acts::SurfaceIntersection Acts::Layer::surfaceOnApproach(
                    (m_ssSensitiveSurfaces > 1 || m_ssApproachSurfaces > 1 ||
                     (surfaceRepresentation().surfaceMaterial() != nullptr));
 
-  // The Limits: current path & overstepping
+  // The Limits
   double nearLimit = options.nearLimit;
   double farLimit = options.farLimit;
 
diff --git a/Core/src/Geometry/TrackingVolume.cpp b/Core/src/Geometry/TrackingVolume.cpp
index 81dd4445a4e4..bef62b1e59c8 100644
--- a/Core/src/Geometry/TrackingVolume.cpp
+++ b/Core/src/Geometry/TrackingVolume.cpp
@@ -24,11 +24,9 @@
 #include <algorithm>
 #include <ostream>
 #include <string>
-#include <tuple>
 #include <utility>
 
 namespace Acts {
-class ISurfaceMaterial;
 
 // constructor for arguments
 TrackingVolume::TrackingVolume(
@@ -431,114 +429,98 @@ TrackingVolume::compatibleBoundaries(const GeometryContext& gctx,
                                      const NavigationOptions<Surface>& options,
                                      const Logger& logger) const {
   ACTS_VERBOSE("Finding compatibleBoundaries");
-  // Loop over boundarySurfaces and calculate the intersection
-  auto excludeObject = options.startObject;
-  boost::container::small_vector<BoundaryIntersection, 4> bIntersections;
 
-  // The Limits: current, path & overstepping
-  double nearLimit = 0;
+  boost::container::small_vector<BoundaryIntersection, 4> intersections;
+
+  // The limits for this navigation step
+  double nearLimit = options.nearLimit;
   double farLimit = options.farLimit;
 
   // Helper function to test intersection
   auto checkIntersection =
-      [&](SurfaceMultiIntersection& smIntersection,
-          const BoundarySurface* bSurface) -> BoundaryIntersection {
-    for (const auto& sIntersection : smIntersection.split()) {
-      if (!sIntersection) {
+      [&](SurfaceMultiIntersection& candidates,
+          const BoundarySurface* boundary) -> BoundaryIntersection {
+    for (const auto& intersection : candidates.split()) {
+      if (!intersection) {
         continue;
       }
 
       if (options.forceIntersectBoundaries) {
         const bool coCriterion =
-            std::abs(sIntersection.pathLength()) < std::abs(nearLimit);
+            std::abs(intersection.pathLength()) < std::abs(nearLimit);
         ACTS_VERBOSE("Forcing intersection with surface "
-                     << bSurface->surfaceRepresentation().geometryId());
+                     << boundary->surfaceRepresentation().geometryId());
         if (coCriterion) {
           ACTS_VERBOSE("Intersection forced successfully ");
           ACTS_VERBOSE("- intersection path length "
-                       << std::abs(sIntersection.pathLength())
+                       << std::abs(intersection.pathLength())
                        << " < overstep limit " << std::abs(nearLimit));
-          return BoundaryIntersection(sIntersection, bSurface);
+          return BoundaryIntersection(intersection, boundary);
         }
         ACTS_VERBOSE("Can't force intersection: ");
         ACTS_VERBOSE("- intersection path length "
-                     << std::abs(sIntersection.pathLength())
+                     << std::abs(intersection.pathLength())
                      << " > overstep limit " << std::abs(nearLimit));
       }
 
       ACTS_VERBOSE("Check intersection with surface "
-                   << bSurface->surfaceRepresentation().geometryId());
-      if (detail::checkIntersection(sIntersection.intersection(), nearLimit,
+                   << boundary->surfaceRepresentation().geometryId());
+      if (detail::checkIntersection(intersection.intersection(), nearLimit,
                                     farLimit, logger)) {
-        return BoundaryIntersection(sIntersection, bSurface);
+        return BoundaryIntersection(intersection, boundary);
       }
     }
 
     ACTS_VERBOSE("No intersection accepted");
-    return BoundaryIntersection(SurfaceIntersection::invalid(), bSurface);
+    return BoundaryIntersection(SurfaceIntersection::invalid(), nullptr);
   };
 
   /// Helper function to process boundary surfaces
   auto processBoundaries =
-      [&](const TrackingVolumeBoundaries& bSurfaces) -> void {
-    ACTS_VERBOSE("Processing boundaries");
+      [&](const TrackingVolumeBoundaries& boundaries) -> void {
     // Loop over the boundary surfaces
-    for (auto& bsIter : bSurfaces) {
+    for (auto& boundary : boundaries) {
       // Get the boundary surface pointer
-      const auto& bSurfaceRep = bsIter->surfaceRepresentation();
-      ACTS_VERBOSE("Consider boundary surface " << bSurfaceRep.geometryId()
-                                                << " :\n"
-                                                << std::tie(bSurfaceRep, gctx));
+      const auto& surface = boundary->surfaceRepresentation();
+      ACTS_VERBOSE("Consider boundary surface " << surface.geometryId());
 
       // Exclude the boundary where you are on
-      if (excludeObject != &bSurfaceRep) {
-        auto bCandidate = bSurfaceRep.intersect(gctx, position, direction,
-                                                options.boundaryCheck);
-        // Intersect and continue
-        auto bIntersection = checkIntersection(bCandidate, bsIter.get());
-        if (bIntersection.first) {
-          ACTS_VERBOSE(" - Proceed with surface");
-          bIntersections.push_back(bIntersection);
-        } else {
-          ACTS_VERBOSE(" - Surface intersecion invalid");
-        }
-      } else {
+      // TODO this is not optimal as we might exit via the same boundary (e.g.
+      // cylinder)
+      if (&surface == options.startObject) {
         ACTS_VERBOSE(" - Surface is excluded surface");
+        continue;
+      }
+
+      auto candidates =
+          surface.intersect(gctx, position, direction, options.boundaryCheck);
+      // Intersect and continue
+      auto intersection = checkIntersection(candidates, boundary.get());
+      if (intersection.first) {
+        ACTS_VERBOSE(" - Proceed with surface");
+        intersections.push_back(intersection);
+      } else {
+        ACTS_VERBOSE(" - Surface intersecion invalid");
       }
     }
   };
 
   // Process the boundaries of the current volume
-  auto& bSurfaces = boundarySurfaces();
-  ACTS_VERBOSE("Volume reports " << bSurfaces.size() << " boundary surfaces");
-  processBoundaries(bSurfaces);
+  const auto& surfaces = boundarySurfaces();
+  ACTS_VERBOSE("Volume reports " << surfaces.size() << " boundary surfaces");
+  processBoundaries(surfaces);
 
   // Process potential boundaries of contained volumes
   auto confinedDenseVolumes = denseVolumes();
   ACTS_VERBOSE("Volume reports " << confinedDenseVolumes.size()
                                  << " confined dense volumes");
-  for (const auto& dv : confinedDenseVolumes) {
-    auto& bSurfacesConfined = dv->boundarySurfaces();
-    ACTS_VERBOSE(" -> " << bSurfacesConfined.size() << " boundary surfaces");
-    processBoundaries(bSurfacesConfined);
+  for (const auto& volume : confinedDenseVolumes) {
+    const auto& surfacesConfined = volume->boundarySurfaces();
+    ACTS_VERBOSE(" -> " << surfacesConfined.size() << " boundary surfaces");
+    processBoundaries(surfacesConfined);
   }
 
-  auto comparator = [](double a, double b) {
-    // sign function would be nice but ...
-    if ((a > 0 && b > 0) || (a < 0 && b < 0)) {
-      return a < b;
-    }
-    if (a > 0) {  // b < 0
-      return true;
-    }
-    return false;
-  };
-
-  std::sort(bIntersections.begin(), bIntersections.end(),
-            [&](const BoundaryIntersection& a, const BoundaryIntersection& b) {
-              return comparator(a.first.pathLength(), b.first.pathLength());
-            });
-  return bIntersections;
+  return intersections;
 }
 
 boost::container::small_vector<LayerIntersection, 10>
@@ -549,41 +531,49 @@ TrackingVolume::compatibleLayers(
   boost::container::small_vector<LayerIntersection, 10> lIntersections;
 
   // the confinedLayers
-  if (m_confinedLayers != nullptr) {
-    // start layer given or not - test layer
-    const Layer* tLayer = options.startObject != nullptr
-                              ? options.startObject
-                              : associatedLayer(gctx, position);
-    while (tLayer != nullptr) {
-      // check if the layer needs resolving
-      // - resolveSensitive -> always take layer if it has a surface array
-      // - resolveMaterial -> always take layer if it has material
-      // - resolvePassive -> always take, unless it's a navigation layer
-      // skip the start object
-      if (tLayer != options.startObject && tLayer->resolve(options)) {
-        // if it's a resolveable start layer, you are by definition on it
-        // layer on approach intersection
-        auto atIntersection =
-            tLayer->surfaceOnApproach(gctx, position, direction, options);
-        auto path = atIntersection.pathLength();
-        bool withinLimit = std::abs(path) <= std::abs(options.farLimit);
-        // Intersection is ok - take it (move to surface on approach)
-        if (atIntersection && withinLimit) {
-          // create a layer intersection
-          lIntersections.push_back(LayerIntersection(atIntersection, tLayer));
-        }
+  if (m_confinedLayers == nullptr) {
+    return {};
+  }
+
+  // start layer given or not - test layer
+  const Layer* tLayer = options.startObject != nullptr
+                            ? static_cast<const Layer*>(options.startObject)
+                            : associatedLayer(gctx, position);
+  while (tLayer != nullptr) {
+    // check if the layer needs resolving
+    // - resolveSensitive -> always take layer if it has a surface array
+    // - resolveMaterial -> always take layer if it has material
+    // - resolvePassive -> always take, unless it's a navigation layer
+    // skip the start object
+    if (tLayer != options.startObject && tLayer->resolve(options)) {
+      // if it's a resolveable start layer, you are by definition on it
+      // layer on approach intersection
+      auto atIntersection =
+          tLayer->surfaceOnApproach(gctx, position, direction, options);
+      // Intersection is ok - take it (move to surface on approach)
+      if (atIntersection) {
+        // create a layer intersection
+        lIntersections.push_back(LayerIntersection(atIntersection, tLayer));
       }
-      // move to next one or break because you reached the end layer
-      tLayer = (tLayer == options.endObject)
-                   ? nullptr
-                   : tLayer->nextLayer(gctx, position, direction);
     }
-    std::sort(lIntersections.begin(), lIntersections.end(),
-              [](const LayerIntersection& a, const LayerIntersection& b) {
-                return SurfaceIntersection::pathLengthOrder(a.first, b.first);
-              });
+    // move to next one or break because you reached the end layer
+    tLayer = (tLayer == options.endObject)
+                 ? nullptr
+                 : tLayer->nextLayer(gctx, position, direction);
   }
-  // and return
+
+  // In case of cylindrical layers we might resolve far intersection solutions
+  // which are not valid for navigation. These are discarded here by checking
+  // against the minimum path length.
+  auto min = std::min_element(
+      lIntersections.begin(), lIntersections.end(),
+      [](const LayerIntersection& a, const LayerIntersection& b) {
+        return a.first.pathLength() < b.first.pathLength();
+      });
+  std::rotate(lIntersections.begin(), min, lIntersections.end());
+  lIntersections.resize(std::distance(min, lIntersections.end()),
+                        {SurfaceIntersection::invalid(), nullptr});
+
   return lIntersections;
 }
 
diff --git a/Examples/Python/src/Propagation.cpp b/Examples/Python/src/Propagation.cpp
index 9876481d5531..a044ac82b0f4 100644
--- a/Examples/Python/src/Propagation.cpp
+++ b/Examples/Python/src/Propagation.cpp
@@ -126,19 +126,7 @@ void addPropagation(Context& ctx) {
 
   {
     auto stepper = py::class_<Acts::EigenStepper<>>(m, "EigenStepper");
-    stepper.def(
-        // Add custom constructor lambda so that not specifying the overstep
-        // limit takes the default from C++ EigenStepper
-        py::init(
-            [](std::shared_ptr<const Acts::MagneticFieldProvider> bField,
-               std::optional<double> overstepLimit) -> Acts::EigenStepper<> {
-              if (overstepLimit) {
-                return {std::move(bField), overstepLimit.value()};
-              } else {
-                return {std::move(bField)};
-              }
-            }),
-        py::arg("bField"), py::arg("overstepLimit") = std::nullopt);
+    stepper.def(py::init<std::shared_ptr<const Acts::MagneticFieldProvider>>());
 
     addPropagator<Acts::EigenStepper<>, Acts::Navigator>(prop, "Eigen");
   }
diff --git a/Examples/Python/tests/root_file_hashes.txt b/Examples/Python/tests/root_file_hashes.txt
index f2f038b2a5e4..fb8472ef6a5e 100644
--- a/Examples/Python/tests/root_file_hashes.txt
+++ b/Examples/Python/tests/root_file_hashes.txt
@@ -29,10 +29,10 @@ test_truth_tracking_kalman[odd-0.0]__performance_track_finder.root: 39aec6316cce
 test_truth_tracking_kalman[odd-1000.0]__trackstates_fitter.root: 72c79be1458c4f9c9a1661778c900f0875d257f2d391c4183a698825448919a1
 test_truth_tracking_kalman[odd-1000.0]__tracksummary_fitter.root: 3d424dec9b172f253c8c4ffbda470f678fd1081a3d36dcfea517ab0f94995ae4
 test_truth_tracking_kalman[odd-1000.0]__performance_track_finder.root: 39aec6316cceb90e314e16b02947faa691c18f57c3a851a25e547a8fc05a4593
-test_truth_tracking_gsf[generic]__trackstates_gsf.root: d4160a87b2f0f21cb19b8e3b0d07f97c343c830a535de4d33a186d5179ab5f06
-test_truth_tracking_gsf[generic]__tracksummary_gsf.root: c110627dd7015c20b302273b5cba49f64a51e3450a722b064be1605f85a6bd6b
-test_truth_tracking_gsf[odd]__trackstates_gsf.root: a3b08c5c0497c85ce2280b4487bbbbb1c4a665d9917c0b1ca040f51524c9d1d4
-test_truth_tracking_gsf[odd]__tracksummary_gsf.root: 84310dd790f3f43b5cf9cee0b08c4deed389171f4caaaa00412f62f7e9944e3b
+test_truth_tracking_gsf[generic]__trackstates_gsf.root: b843bcac04e7647aef716bd86f865ae72fa01448705ae933fd75ea5e340df95c
+test_truth_tracking_gsf[generic]__tracksummary_gsf.root: 7de085a3b03136fc7ccf80410457a4ea30dcd702463e7c4e806b4a85a51c05ac
+test_truth_tracking_gsf[odd]__trackstates_gsf.root: ef12f504effc8dd00c378fde348291ce6d7e264c0f186b7ec82c8262f35c9f71
+test_truth_tracking_gsf[odd]__tracksummary_gsf.root: 07c72a4fe5383aeeada41ded42fc3337db67e1b75eccb9e85dfe783379869566
 test_particle_gun__particles.root: 5fe7dda2933ee6b9615b064d192322fe07831133cd998e5ed99a3b992b713a10
 test_material_mapping__material-map_tracks.root: b1138566d8d51579dce07f80166f05986942cf78c76b36875aea9b4b8b9bb957
 test_material_mapping__propagation-material.root: 2ec28ec8c10860ce0ca03b0b8abbfc9b52cef779d0113fddfe05949a8c362ec0
diff --git a/Tests/UnitTests/Core/Propagator/AtlasStepperTests.cpp b/Tests/UnitTests/Core/Propagator/AtlasStepperTests.cpp
index a2895ebb50a5..f0cb47ff0fa0 100644
--- a/Tests/UnitTests/Core/Propagator/AtlasStepperTests.cpp
+++ b/Tests/UnitTests/Core/Propagator/AtlasStepperTests.cpp
@@ -557,9 +557,6 @@ BOOST_AUTO_TEST_CASE(StepSize) {
                                  particleHypothesis),
       stepSize, tolerance);
 
-  // TODO figure out why this fails and what it should be
-  // BOOST_CHECK_EQUAL(stepper.overstepLimit(state), tolerance);
-
   stepper.updateStepSize(state, -5_cm, ConstrainedStep::actor);
   BOOST_CHECK_EQUAL(state.previousStepSize, stepSize);
   BOOST_CHECK_EQUAL(state.stepSize.value(), -5_cm);
diff --git a/Tests/UnitTests/Core/Propagator/EigenStepperTests.cpp b/Tests/UnitTests/Core/Propagator/EigenStepperTests.cpp
index 34b4b4a33481..7ca61a3587b7 100644
--- a/Tests/UnitTests/Core/Propagator/EigenStepperTests.cpp
+++ b/Tests/UnitTests/Core/Propagator/EigenStepperTests.cpp
@@ -681,9 +681,7 @@ BOOST_AUTO_TEST_CASE(step_extension_material_test) {
       tgb.trackingGeometry(tgContext);
 
   // Build navigator
-  Navigator naviMat(
-      {material, true, true, true},
-      Acts::getDefaultLogger("Navigator", Acts::Logging::VERBOSE));
+  Navigator naviMat({material, true, true, true});
 
   // Set initial parameters for the particle track
   Covariance cov = Covariance::Identity();
@@ -919,7 +917,7 @@ BOOST_AUTO_TEST_CASE(step_extension_vacmatvac_test) {
 
   PropagatorOptions<ActionList<StepCollector>, AbortList<EndOfWorld>>
       propOptsDef(tgContext, mfContext);
-  abortList.get<EndOfWorld>().maxX = 1_m;
+  abortList.get<EndOfWorld>().maxX = 3_m;
   propOptsDef.abortList = abortList;
   propOptsDef.maxSteps = 1000;
   propOptsDef.maxStepSize = 1.5_m;
@@ -930,8 +928,7 @@ BOOST_AUTO_TEST_CASE(step_extension_vacmatvac_test) {
       propDef(esDef, naviDet);
 
   // Launch and collect results
-  const auto& resultDef =
-      propDef.propagate(sbtp, *(surs[0]), propOptsDef).value();
+  const auto& resultDef = propDef.propagate(sbtp, propOptsDef).value();
   const StepCollector::this_result& stepResultDef =
       resultDef.get<typename StepCollector::result_type>();
 
@@ -965,41 +962,40 @@ BOOST_AUTO_TEST_CASE(step_extension_vacmatvac_test) {
   // Build launcher through material
   // Set initial parameters for the particle track by using the result of the
   // first volume
-  CurvilinearTrackParameters sbtpPiecewise(
-      makeVector4(endParams.first, 0), endParams.second,
-      1_e / endParams.second.norm(), std::nullopt, ParticleHypothesis::pion());
 
   // Set options for propagator
   DenseStepperPropagatorOptions<ActionList<StepCollector>,
                                 AbortList<EndOfWorld>>
       propOptsDense(tgContext, mfContext);
-  abortList.get<EndOfWorld>().maxX = 2_m;
+  abortList.get<EndOfWorld>().maxX = 3_m;
   propOptsDense.abortList = abortList;
   propOptsDense.maxSteps = 1000;
   propOptsDense.maxStepSize = 1.5_m;
 
   // Build stepper and propagator
-  EigenStepper<StepperExtensionList<DenseEnvironmentExtension>> esDense(bField);
-  Propagator<EigenStepper<StepperExtensionList<DenseEnvironmentExtension>>,
+  EigenStepper<
+      StepperExtensionList<DefaultExtension, DenseEnvironmentExtension>>
+      esDense(bField);
+  Propagator<EigenStepper<StepperExtensionList<DefaultExtension,
+                                               DenseEnvironmentExtension>>,
              Navigator>
       propDense(esDense, naviDet);
 
   // Launch and collect results
-  const auto& resultDense =
-      propDense.propagate(sbtpPiecewise, *(surs[1]), propOptsDense).value();
+  const auto& resultDense = propDense.propagate(sbtp, propOptsDense).value();
   const StepCollector::this_result& stepResultDense =
       resultDense.get<typename StepCollector::result_type>();
 
   // Check the exit situation of the second volume
   for (unsigned int i = 0; i < stepResultDense.position.size(); i++) {
-    if (2_m - stepResultDense.position[i].x() < 1e-4) {
+    if (1_m - stepResultDense.position[i].x() < 1e-4) {
       endParams = std::make_pair(stepResultDense.position[i],
                                  stepResultDense.momentum[i]);
       break;
     }
   }
   for (unsigned int i = 0; i < stepResult.position.size(); i++) {
-    if (2_m - stepResult.position[i].x() < 1e-4) {
+    if (1_m - stepResult.position[i].x() < 1e-4) {
       endParamsControl =
           std::make_pair(stepResult.position[i], stepResult.momentum[i]);
       break;
diff --git a/Tests/UnitTests/Core/Propagator/MaterialCollectionTests.cpp b/Tests/UnitTests/Core/Propagator/MaterialCollectionTests.cpp
index f14679f77db7..ea9c06d72423 100644
--- a/Tests/UnitTests/Core/Propagator/MaterialCollectionTests.cpp
+++ b/Tests/UnitTests/Core/Propagator/MaterialCollectionTests.cpp
@@ -56,10 +56,6 @@ MagneticFieldContext mfContext = MagneticFieldContext();
 CylindricalTrackingGeometry cGeometry(tgContext);
 auto tGeometry = cGeometry();
 
-// create a navigator for this tracking geometry
-Navigator navigatorES({tGeometry});
-Navigator navigatorSL({tGeometry});
-
 using BField = ConstantBField;
 using EigenStepper = Acts::EigenStepper<>;
 using EigenPropagator = Propagator<EigenStepper, Navigator>;
@@ -67,60 +63,37 @@ using StraightLinePropagator = Propagator<StraightLineStepper, Navigator>;
 
 const double Bz = 2_T;
 auto bField = std::make_shared<BField>(Vector3{0, 0, Bz});
+
 EigenStepper estepper(bField);
-EigenPropagator epropagator(std::move(estepper), std::move(navigatorES));
+Navigator esnavigator({tGeometry});
+EigenPropagator epropagator(std::move(estepper), std::move(esnavigator));
 
 StraightLineStepper slstepper;
-StraightLinePropagator slpropagator(slstepper, std::move(navigatorSL));
-const int ntests = 500;
-const int skip = 0;
-bool debugModeFwd = false;
-bool debugModeBwd = false;
-bool debugModeFwdStep = false;
-bool debugModeBwdStep = false;
-
-/// the actual test nethod that runs the test
-/// can be used with several propagator types
+Navigator slnavigator({tGeometry});
+StraightLinePropagator slpropagator(slstepper, std::move(slnavigator));
+
+int ntests = 500;
+int skip = 0;
+bool debugMode = false;
+
+/// the actual test nethod that runs the test can be used with several
+/// propagator types
+///
 /// @tparam propagator_t is the actual propagator type
 ///
 /// @param prop is the propagator instance
-/// @param pT the transverse momentum
-/// @param phi the azimuthal angle of the track at creation
-/// @param theta the polar angle of the track at creation
-/// @param charge is the charge of the particle
-/// @param index is the run index from the test
+/// @param start the start parameters
 template <typename propagator_t>
-void runTest(const propagator_t& prop, double pT, double phi, double theta,
-             int charge, int index) {
-  double p = pT / sin(theta);
-  double q = -1 + 2 * charge;
-
-  if (index < skip) {
-    return;
-  }
-
-  // define start parameters
-  BoundSquareMatrix cov;
-  // take some major correlations (off-diagonals)
-  // clang-format off
-    cov <<
-     10_mm, 0, 0.123, 0, 0.5, 0,
-     0, 10_mm, 0, 0.162, 0, 0,
-     0.123, 0, 0.1, 0, 0, 0,
-     0, 0.162, 0, 0.1, 0, 0,
-     0.5, 0, 0, 0, 1_e / 10_GeV, 0,
-     0, 0, 0, 0, 0, 1_us;
-  // clang-format on
-  CurvilinearTrackParameters start(Vector4(0, 0, 0, 0), phi, theta, q / p, cov,
-                                   ParticleHypothesis::pion());
-
+void runTest(const propagator_t& prop,
+             const CurvilinearTrackParameters& start) {
   // Action list and abort list
   using ActionListType = ActionList<MaterialInteractor>;
   using AbortListType = AbortList<>;
 
   using Options = PropagatorOptions<ActionListType, AbortListType>;
-  Options fwdOptions(tgContext, mfContext);
 
+  // forward material test
+  Options fwdOptions(tgContext, mfContext);
   fwdOptions.maxStepSize = 25_cm;
   fwdOptions.pathLimit = 25_cm;
 
@@ -131,12 +104,13 @@ void runTest(const propagator_t& prop, double pT, double phi, double theta,
   fwdMaterialInteractor.energyLoss = false;
   fwdMaterialInteractor.multipleScattering = false;
 
-  if (debugModeFwd) {
+  if (debugMode) {
     std::cout << ">>> Forward Propagation : start." << std::endl;
   }
   // forward material test
   const auto& fwdResult = prop.propagate(start, fwdOptions).value();
-  auto& fwdMaterial = fwdResult.template get<MaterialInteractor::result_type>();
+  const auto& fwdMaterial =
+      fwdResult.template get<MaterialInteractor::result_type>();
   // check that the collected material is not zero
   BOOST_CHECK_NE(fwdMaterial.materialInX0, 0.);
   BOOST_CHECK_NE(fwdMaterial.materialInL0, 0.);
@@ -152,7 +126,7 @@ void runTest(const propagator_t& prop, double pT, double phi, double theta,
   CHECK_CLOSE_REL(fwdMaterial.materialInL0, fwdStepMaterialInL0, 1e-3);
 
   // get the forward output to the screen
-  if (debugModeFwd) {
+  if (debugMode) {
     // check if the surfaces are free
     std::cout << ">>> Material steps found on ..." << std::endl;
     for (auto& fwdStepsC : fwdMaterial.materialInteractions) {
@@ -176,18 +150,16 @@ void runTest(const propagator_t& prop, double pT, double phi, double theta,
 
   const auto& startSurface = start.referenceSurface();
 
-  if (debugModeBwd) {
+  if (debugMode) {
     std::cout << ">>> Backward Propagation : start." << std::endl;
   }
   const auto& bwdResult =
       prop.propagate(*fwdResult.endParameters, startSurface, bwdOptions)
           .value();
-
-  if (debugModeBwd) {
+  if (debugMode) {
     std::cout << ">>> Backward Propagation : end." << std::endl;
   }
-
-  auto& bwdMaterial =
+  const auto& bwdMaterial =
       bwdResult.template get<typename MaterialInteractor::result_type>();
   // check that the collected material is not zero
   BOOST_CHECK_NE(bwdMaterial.materialInX0, 0.);
@@ -204,7 +176,7 @@ void runTest(const propagator_t& prop, double pT, double phi, double theta,
   CHECK_CLOSE_REL(bwdMaterial.materialInL0, bwdStepMaterialInL0, 1e-3);
 
   // get the backward output to the screen
-  if (debugModeBwd) {
+  if (debugMode) {
     // check if the surfaces are free
     std::cout << ">>> Material steps found on ..." << std::endl;
     for (auto& bwdStepsC : bwdMaterial.materialInteractions) {
@@ -218,7 +190,7 @@ void runTest(const propagator_t& prop, double pT, double phi, double theta,
                     fwdMaterial.materialInteractions.size());
 
   CHECK_CLOSE_REL(bwdMaterial.materialInX0, fwdMaterial.materialInX0, 1e-3);
-  CHECK_CLOSE_REL(bwdMaterial.materialInL0, bwdMaterial.materialInL0, 1e-3);
+  CHECK_CLOSE_REL(bwdMaterial.materialInL0, fwdMaterial.materialInL0, 1e-3);
 
   // stepping from one surface to the next
   // now go from surface to surface and check
@@ -236,7 +208,7 @@ void runTest(const propagator_t& prop, double pT, double phi, double theta,
   double fwdStepStepMaterialInX0 = 0.;
   double fwdStepStepMaterialInL0 = 0.;
 
-  if (debugModeFwdStep) {
+  if (debugMode) {
     // check if the surfaces are free
     std::cout << ">>> Forward steps to be processed sequentially ..."
               << std::endl;
@@ -250,7 +222,7 @@ void runTest(const propagator_t& prop, double pT, double phi, double theta,
   BoundTrackParameters sParameters = start;
   std::vector<BoundTrackParameters> stepParameters;
   for (auto& fwdSteps : fwdMaterial.materialInteractions) {
-    if (debugModeFwdStep) {
+    if (debugMode) {
       std::cout << ">>> Forward step : "
                 << sParameters.referenceSurface().geometryId() << " --> "
                 << fwdSteps.surface->geometryId() << std::endl;
@@ -275,7 +247,7 @@ void runTest(const propagator_t& prop, double pT, double phi, double theta,
   // final destination surface
   const Surface& dSurface = fwdResult.endParameters->referenceSurface();
 
-  if (debugModeFwdStep) {
+  if (debugMode) {
     std::cout << ">>> Forward step : "
               << sParameters.referenceSurface().geometryId() << " --> "
               << dSurface.geometryId() << std::endl;
@@ -296,7 +268,6 @@ void runTest(const propagator_t& prop, double pT, double phi, double theta,
   // stepping from one surface to the next : backwards
   // now go from surface to surface and check
   Options bwdStepOptions(tgContext, mfContext);
-
   bwdStepOptions.maxStepSize = 25_cm;
   bwdStepOptions.pathLimit = -25_cm;
   bwdStepOptions.direction = Direction::Backward;
@@ -311,7 +282,7 @@ void runTest(const propagator_t& prop, double pT, double phi, double theta,
   double bwdStepStepMaterialInX0 = 0.;
   double bwdStepStepMaterialInL0 = 0.;
 
-  if (debugModeBwdStep) {
+  if (debugMode) {
     // check if the surfaces are free
     std::cout << ">>> Backward steps to be processed sequentially ..."
               << std::endl;
@@ -324,7 +295,7 @@ void runTest(const propagator_t& prop, double pT, double phi, double theta,
   // move forward step by step through the surfaces
   sParameters = *fwdResult.endParameters;
   for (auto& bwdSteps : bwdMaterial.materialInteractions) {
-    if (debugModeBwdStep) {
+    if (debugMode) {
       std::cout << ">>> Backward step : "
                 << sParameters.referenceSurface().geometryId() << " --> "
                 << bwdSteps.surface->geometryId() << std::endl;
@@ -348,7 +319,7 @@ void runTest(const propagator_t& prop, double pT, double phi, double theta,
   // final destination surface
   const Surface& dbSurface = start.referenceSurface();
 
-  if (debugModeBwdStep) {
+  if (debugMode) {
     std::cout << ">>> Backward step : "
               << sParameters.referenceSurface().geometryId() << " --> "
               << dSurface.geometryId() << std::endl;
@@ -378,7 +349,7 @@ void runTest(const propagator_t& prop, double pT, double phi, double theta,
   const auto& covfwdResult = prop.propagate(start, fwdOptions).value();
 
   BOOST_CHECK_LE(
-      cov.determinant(),
+      start.covariance()->determinant(),
       covfwdResult.endParameters->covariance().value().determinant());
 }
 
@@ -402,8 +373,30 @@ BOOST_DATA_TEST_CASE(
                            std::uniform_int_distribution<std::uint8_t>(0, 1))) ^
         bdata::xrange(ntests),
     pT, phi, theta, charge, index) {
-  runTest(epropagator, pT, phi, theta, charge, index);
-  runTest(slpropagator, pT, phi, theta, charge, index);
+  if (index < skip) {
+    return;
+  }
+
+  double p = pT / sin(theta);
+  double q = -1 + 2 * charge;
+
+  // define start parameters
+  BoundSquareMatrix cov;
+  // take some major correlations (off-diagonals)
+  // clang-format off
+    cov <<
+     10_mm, 0, 0.123, 0, 0.5, 0,
+     0, 10_mm, 0, 0.162, 0, 0,
+     0.123, 0, 0.1, 0, 0, 0,
+     0, 0.162, 0, 0.1, 0, 0,
+     0.5, 0, 0, 0, 1_e / 10_GeV, 0,
+     0, 0, 0, 0, 0, 1_us;
+  // clang-format on
+  CurvilinearTrackParameters start(Vector4(0, 0, 0, 0), phi, theta, q / p, cov,
+                                   ParticleHypothesis::pion());
+
+  runTest(epropagator, start);
+  runTest(slpropagator, start);
 }
 
 }  // namespace Acts::Test
diff --git a/Tests/UnitTests/Core/Propagator/MultiStepperTests.cpp b/Tests/UnitTests/Core/Propagator/MultiStepperTests.cpp
index 9d6f13e0e332..4ecdb3ca59ab 100644
--- a/Tests/UnitTests/Core/Propagator/MultiStepperTests.cpp
+++ b/Tests/UnitTests/Core/Propagator/MultiStepperTests.cpp
@@ -399,11 +399,14 @@ void test_multi_stepper_surface_status_update() {
     BOOST_CHECK_EQUAL(status, Intersection3D::Status::reachable);
 
     auto cmp_iterable = multi_stepper.constComponentIterable(multi_state);
+    auto cmp_1 = *cmp_iterable.begin();
+    auto cmp_2 = *(++cmp_iterable.begin());
 
-    BOOST_CHECK_EQUAL((*cmp_iterable.begin()).status(),
-                      Intersection3D::Status::reachable);
-    BOOST_CHECK_EQUAL((*(++cmp_iterable.begin())).status(),
-                      Intersection3D::Status::missed);
+    BOOST_CHECK_EQUAL(cmp_1.status(), Intersection3D::Status::reachable);
+    BOOST_CHECK_EQUAL(cmp_2.status(), Intersection3D::Status::reachable);
+
+    BOOST_CHECK_EQUAL(cmp_1.cmp.state.stepSize.value(), 1.0);
+    BOOST_CHECK_EQUAL(cmp_2.cmp.state.stepSize.value(), -1.0);
   }
 
   // Step forward now
@@ -425,27 +428,30 @@ void test_multi_stepper_surface_status_update() {
     BOOST_CHECK_EQUAL(status, Intersection3D::Status::onSurface);
 
     auto cmp_iterable = multi_stepper.constComponentIterable(multi_state);
+    auto cmp_1 = *cmp_iterable.begin();
+    auto cmp_2 = *(++cmp_iterable.begin());
 
-    BOOST_CHECK_EQUAL((*cmp_iterable.begin()).status(),
-                      Intersection3D::Status::onSurface);
-    BOOST_CHECK_EQUAL((*(++cmp_iterable.begin())).status(),
-                      Intersection3D::Status::missed);
+    BOOST_CHECK_EQUAL(cmp_1.status(), Intersection3D::Status::onSurface);
+    BOOST_CHECK_EQUAL(cmp_2.status(), Intersection3D::Status::onSurface);
   }
 
-  // Start surface should be unreachable
+  // Start surface should be reachable
   {
     auto status = multi_stepper.updateSurfaceStatus(multi_state, *start_surface,
                                                     0, Direction::Forward,
                                                     BoundaryCheck(false));
 
-    BOOST_CHECK_EQUAL(status, Intersection3D::Status::unreachable);
+    BOOST_CHECK_EQUAL(status, Intersection3D::Status::reachable);
 
     auto cmp_iterable = multi_stepper.constComponentIterable(multi_state);
+    auto cmp_1 = *cmp_iterable.begin();
+    auto cmp_2 = *(++cmp_iterable.begin());
+
+    BOOST_CHECK_EQUAL(cmp_1.status(), Intersection3D::Status::reachable);
+    BOOST_CHECK_EQUAL(cmp_2.status(), Intersection3D::Status::reachable);
 
-    BOOST_CHECK_EQUAL((*cmp_iterable.begin()).status(),
-                      Intersection3D::Status::unreachable);
-    BOOST_CHECK_EQUAL((*(++cmp_iterable.begin())).status(),
-                      Intersection3D::Status::unreachable);
+    BOOST_CHECK_EQUAL(cmp_1.cmp.state.stepSize.value(), -1.0);
+    BOOST_CHECK_EQUAL(cmp_2.cmp.state.stepSize.value(), 1.0);
   }
 }
 
@@ -515,17 +521,17 @@ void test_component_bound_state() {
     BOOST_REQUIRE(single_bound_state.ok());
 
     auto cmp_iterable = multi_stepper.componentIterable(multi_state);
+    auto cmp_1 = *cmp_iterable.begin();
+    auto cmp_2 = *(++cmp_iterable.begin());
 
-    auto ok_bound_state =
-        (*cmp_iterable.begin())
-            .boundState(*right_surface, true, FreeToBoundCorrection(false));
-    BOOST_REQUIRE(ok_bound_state.ok());
-    BOOST_CHECK(*single_bound_state == *ok_bound_state);
+    auto bound_state_1 =
+        cmp_1.boundState(*right_surface, true, FreeToBoundCorrection(false));
+    BOOST_REQUIRE(bound_state_1.ok());
+    BOOST_CHECK(*single_bound_state == *bound_state_1);
 
-    auto failed_bound_state =
-        (*(++cmp_iterable.begin()))
-            .boundState(*right_surface, true, FreeToBoundCorrection(false));
-    BOOST_CHECK(!failed_bound_state.ok());
+    auto bound_state_2 =
+        cmp_2.boundState(*right_surface, true, FreeToBoundCorrection(false));
+    BOOST_CHECK(bound_state_2.ok());
   }
 }
 
diff --git a/Tests/UnitTests/Core/Propagator/NavigatorTests.cpp b/Tests/UnitTests/Core/Propagator/NavigatorTests.cpp
index 3f0aec61513c..9c8f75f25083 100644
--- a/Tests/UnitTests/Core/Propagator/NavigatorTests.cpp
+++ b/Tests/UnitTests/Core/Propagator/NavigatorTests.cpp
@@ -21,15 +21,18 @@
 #include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/Geometry/TrackingGeometry.hpp"
 #include "Acts/Geometry/TrackingVolume.hpp"
+#include "Acts/MagneticField/ConstantBField.hpp"
 #include "Acts/Propagator/ConstrainedStep.hpp"
+#include "Acts/Propagator/EigenStepper.hpp"
 #include "Acts/Propagator/Navigator.hpp"
 #include "Acts/Propagator/StepperConcept.hpp"
+#include "Acts/Propagator/StraightLineStepper.hpp"
+#include "Acts/Propagator/SurfaceCollector.hpp"
 #include "Acts/Propagator/detail/SteppingHelper.hpp"
 #include "Acts/Surfaces/BoundaryCheck.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Tests/CommonHelpers/CylindricalTrackingGeometry.hpp"
 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
-#include "Acts/Utilities/Helpers.hpp"
 #include "Acts/Utilities/Intersection.hpp"
 #include "Acts/Utilities/Logger.hpp"
 #include "Acts/Utilities/Result.hpp"
@@ -37,7 +40,6 @@
 #include <cstddef>
 #include <cstdint>
 #include <iostream>
-#include <map>
 #include <memory>
 #include <string>
 #include <system_error>
@@ -49,6 +51,7 @@ class Layer;
 struct FreeToBoundCorrection;
 }  // namespace Acts
 
+namespace bdata = boost::unit_test::data;
 using namespace Acts::UnitLiterals;
 using Acts::VectorHelpers::perp;
 
@@ -56,6 +59,7 @@ namespace Acts::Test {
 
 // Create a test context
 GeometryContext tgContext = GeometryContext();
+MagneticFieldContext mfContext = MagneticFieldContext();
 
 /// This is a simple cache struct to mimic the
 /// Propagator cache
@@ -319,6 +323,9 @@ bool testNavigatorStatePointers(
 CylindricalTrackingGeometry cGeometry(tgContext);
 auto tGeometry = cGeometry();
 
+const double Bz = 2_T;
+auto bField = std::make_shared<ConstantBField>(Vector3{0, 0, Bz});
+
 // the debug boolean
 bool debug = true;
 
@@ -692,4 +699,361 @@ BOOST_AUTO_TEST_CASE(Navigator_target_methods) {
   }
 }
 
+using SurfaceCollector = SurfaceCollector<SurfaceSelector>;
+
+std::vector<GeometryIdentifier> collectRelevantGeoIds(
+    const SurfaceCollector::result_type& surfaceHits) {
+  std::vector<GeometryIdentifier> geoIds;
+  for (const auto& surfaceHit : surfaceHits.collected) {
+    auto geoId = surfaceHit.surface->geometryId();
+    auto material = surfaceHit.surface->surfaceMaterial();
+    if (geoId.sensitive() == 0 && material == nullptr) {
+      continue;
+    }
+    geoIds.push_back(geoId);
+  }
+  return geoIds;
+}
+
+/// the actual test method that runs the test can be used with several
+/// propagator types
+///
+/// @tparam propagator_t is the actual propagator type
+///
+/// @param prop is the propagator instance
+/// @param start start parameters for propagation
+/// @param debugMode toggle debug mode
+template <typename propagator_t>
+void runSelfConsistencyTest(const propagator_t& prop,
+                            const CurvilinearTrackParameters& start,
+                            bool debugMode) {
+  // Action list and abort list
+  using ActionListType = ActionList<SurfaceCollector>;
+  using AbortListType = AbortList<>;
+  using Options = PropagatorOptions<ActionListType, AbortListType>;
+
+  // forward surface test
+  Options fwdOptions(tgContext, mfContext);
+  fwdOptions.pathLimit = 25_cm;
+  fwdOptions.maxStepSize = 1_cm;
+
+  // get the surface collector and configure it
+  auto& fwdSurfaceCollector =
+      fwdOptions.actionList.template get<SurfaceCollector>();
+  fwdSurfaceCollector.selector.selectSensitive = true;
+  fwdSurfaceCollector.selector.selectMaterial = true;
+  fwdSurfaceCollector.selector.selectPassive = true;
+
+  if (debugMode) {
+    std::cout << ">>> Forward Propagation : start." << std::endl;
+  }
+  auto fwdResult = prop.propagate(start, fwdOptions).value();
+  auto fwdSurfaceHits =
+      fwdResult.template get<SurfaceCollector::result_type>().collected;
+  auto fwdSurfaces = collectRelevantGeoIds(
+      fwdResult.template get<SurfaceCollector::result_type>());
+
+  // get the forward output to the screen
+  if (debugMode) {
+    // check if the surfaces are free
+    std::cout << ">>> Surface hits found on ..." << std::endl;
+    for (const auto& fwdSteps : fwdSurfaces) {
+      std::cout << "--> Surface with " << fwdSteps << std::endl;
+    }
+    std::cout << ">>> Forward Propagation : end." << std::endl;
+  }
+
+  // backward surface test
+  Options bwdOptions(tgContext, mfContext);
+  bwdOptions.pathLimit = 25_cm;
+  bwdOptions.maxStepSize = 1_cm;
+  bwdOptions.direction = Direction::Backward;
+
+  // get the surface collector and configure it
+  auto& bwdMSurfaceCollector =
+      bwdOptions.actionList.template get<SurfaceCollector>();
+  bwdMSurfaceCollector.selector.selectSensitive = true;
+  bwdMSurfaceCollector.selector.selectMaterial = true;
+  bwdMSurfaceCollector.selector.selectPassive = true;
+
+  const auto& startSurface = start.referenceSurface();
+
+  if (debugMode) {
+    std::cout << ">>> Backward Propagation : start." << std::endl;
+  }
+  auto bwdResult =
+      prop.propagate(*fwdResult.endParameters, startSurface, bwdOptions)
+          .value();
+  auto bwdSurfaceHits =
+      bwdResult.template get<SurfaceCollector::result_type>().collected;
+  auto bwdSurfaces = collectRelevantGeoIds(
+      bwdResult.template get<SurfaceCollector::result_type>());
+
+  // get the backward output to the screen
+  if (debugMode) {
+    // check if the surfaces are free
+    std::cout << ">>> Surface hits found on ..." << std::endl;
+    for (auto& bwdSteps : bwdSurfaces) {
+      std::cout << "--> Surface with " << bwdSteps << std::endl;
+    }
+    std::cout << ">>> Backward Propagation : end." << std::endl;
+  }
+
+  // forward-backward compatibility test
+  {
+    std::reverse(bwdSurfaces.begin(), bwdSurfaces.end());
+    BOOST_CHECK_EQUAL_COLLECTIONS(bwdSurfaces.begin(), bwdSurfaces.end(),
+                                  fwdSurfaces.begin(), fwdSurfaces.end());
+  }
+
+  // stepping from one surface to the next
+  // now go from surface to surface and check
+  Options fwdStepOptions(tgContext, mfContext);
+  fwdStepOptions.maxStepSize = 1_cm;
+
+  // get the surface collector and configure it
+  auto& fwdStepSurfaceCollector =
+      fwdOptions.actionList.template get<SurfaceCollector>();
+  fwdStepSurfaceCollector.selector.selectSensitive = true;
+  fwdStepSurfaceCollector.selector.selectMaterial = true;
+  fwdStepSurfaceCollector.selector.selectPassive = true;
+
+  std::vector<GeometryIdentifier> fwdStepSurfaces;
+
+  // move forward step by step through the surfaces
+  BoundTrackParameters sParameters = start;
+  std::vector<BoundTrackParameters> stepParameters;
+  for (auto& fwdSteps : fwdSurfaceHits) {
+    if (debugMode) {
+      std::cout << ">>> Forward step : "
+                << sParameters.referenceSurface().geometryId() << " --> "
+                << fwdSteps.surface->geometryId() << std::endl;
+    }
+
+    // make a forward step
+    auto fwdStep =
+        prop.propagate(sParameters, *fwdSteps.surface, fwdStepOptions).value();
+
+    auto fwdStepSurfacesTmp = collectRelevantGeoIds(
+        fwdStep.template get<SurfaceCollector::result_type>());
+    fwdStepSurfaces.insert(fwdStepSurfaces.end(), fwdStepSurfacesTmp.begin(),
+                           fwdStepSurfacesTmp.end());
+
+    if (fwdStep.endParameters.has_value()) {
+      // make sure the parameters do not run out of scope
+      stepParameters.push_back(*fwdStep.endParameters);
+      sParameters = stepParameters.back();
+    }
+  }
+  // final destination surface
+  const Surface& dSurface = fwdResult.endParameters->referenceSurface();
+  if (debugMode) {
+    std::cout << ">>> Forward step : "
+              << sParameters.referenceSurface().geometryId() << " --> "
+              << dSurface.geometryId() << std::endl;
+  }
+  auto fwdStepFinal =
+      prop.propagate(sParameters, dSurface, fwdStepOptions).value();
+  auto fwdStepSurfacesTmp = collectRelevantGeoIds(
+      fwdStepFinal.template get<SurfaceCollector::result_type>());
+  fwdStepSurfaces.insert(fwdStepSurfaces.end(), fwdStepSurfacesTmp.begin(),
+                         fwdStepSurfacesTmp.end());
+
+  // TODO forward-forward step compatibility test
+
+  // stepping from one surface to the next : backwards
+  // now go from surface to surface and check
+  Options bwdStepOptions(tgContext, mfContext);
+  bwdStepOptions.maxStepSize = 1_cm;
+  bwdStepOptions.direction = Direction::Backward;
+
+  // get the surface collector and configure it
+  auto& bwdStepSurfaceCollector =
+      bwdOptions.actionList.template get<SurfaceCollector>();
+  bwdStepSurfaceCollector.selector.selectSensitive = true;
+  bwdStepSurfaceCollector.selector.selectMaterial = true;
+  bwdStepSurfaceCollector.selector.selectPassive = true;
+
+  std::vector<GeometryIdentifier> bwdStepSurfaces;
+
+  // move forward step by step through the surfaces
+  sParameters = *fwdResult.endParameters;
+  for (auto& bwdSteps : bwdSurfaceHits) {
+    if (debugMode) {
+      std::cout << ">>> Backward step : "
+                << sParameters.referenceSurface().geometryId() << " --> "
+                << bwdSteps.surface->geometryId() << std::endl;
+    }
+
+    // make a forward step
+    auto bwdStep =
+        prop.propagate(sParameters, *bwdSteps.surface, bwdStepOptions).value();
+
+    auto bwdStepSurfacesTmp = collectRelevantGeoIds(
+        bwdStep.template get<SurfaceCollector::result_type>());
+    bwdStepSurfaces.insert(bwdStepSurfaces.end(), bwdStepSurfacesTmp.begin(),
+                           bwdStepSurfacesTmp.end());
+
+    if (bwdStep.endParameters.has_value()) {
+      // make sure the parameters do not run out of scope
+      stepParameters.push_back(*bwdStep.endParameters);
+      sParameters = stepParameters.back();
+    }
+  }
+  // final destination surface
+  const Surface& dbSurface = start.referenceSurface();
+  if (debugMode) {
+    std::cout << ">>> Backward step : "
+              << sParameters.referenceSurface().geometryId() << " --> "
+              << dSurface.geometryId() << std::endl;
+  }
+  auto bwdStepFinal =
+      prop.propagate(sParameters, dbSurface, bwdStepOptions).value();
+  auto bwdStepSurfacesTmp = collectRelevantGeoIds(
+      bwdStepFinal.template get<SurfaceCollector::result_type>());
+  bwdStepSurfaces.insert(bwdStepSurfaces.end(), bwdStepSurfacesTmp.begin(),
+                         bwdStepSurfacesTmp.end());
+
+  // TODO backward-backward step compatibility test
+
+  std::reverse(bwdStepSurfaces.begin(), bwdStepSurfaces.end());
+  BOOST_CHECK_EQUAL_COLLECTIONS(bwdStepSurfaces.begin(), bwdStepSurfaces.end(),
+                                fwdStepSurfaces.begin(), fwdStepSurfaces.end());
+}
+
+/// the actual test method that runs the test can be used with several
+/// propagator types
+///
+/// @tparam propagator_probe_t is the probe propagator type
+/// @tparam propagator_ref_t is the reference propagator type
+///
+/// @param propProbe is the probe propagator instance
+/// @param propRef is the reference propagator instance
+/// @param start start parameters for propagation
+/// @param debugMode toggle debug mode
+template <typename propagator_probe_t, typename propagator_ref_t>
+void runConsistencyTest(const propagator_probe_t& propProbe,
+                        const propagator_ref_t& propRef,
+                        const CurvilinearTrackParameters& start,
+                        bool debugMode) {
+  // Action list and abort list
+  using ActionListType = ActionList<SurfaceCollector>;
+  using AbortListType = AbortList<>;
+  using Options = PropagatorOptions<ActionListType, AbortListType>;
+
+  auto run = [&](const auto& prop) {
+    // forward surface test
+    Options fwdOptions(tgContext, mfContext);
+    fwdOptions.pathLimit = 25_cm;
+    fwdOptions.maxStepSize = 1_cm;
+
+    // get the surface collector and configure it
+    auto& fwdSurfaceCollector =
+        fwdOptions.actionList.template get<SurfaceCollector>();
+    fwdSurfaceCollector.selector.selectSensitive = true;
+    fwdSurfaceCollector.selector.selectMaterial = true;
+    fwdSurfaceCollector.selector.selectPassive = true;
+
+    auto fwdResult = prop.propagate(start, fwdOptions).value();
+    auto fwdSurfaces = collectRelevantGeoIds(
+        fwdResult.template get<SurfaceCollector::result_type>());
+
+    // get the forward output to the screen
+    if (debugMode) {
+      // check if the surfaces are free
+      std::cout << ">>> Surface hits found on ..." << std::endl;
+      for (const auto& fwdSteps : fwdSurfaces) {
+        std::cout << "--> Surface with " << fwdSteps << std::endl;
+      }
+    }
+
+    return fwdSurfaces;
+  };
+
+  if (debugMode) {
+    std::cout << ">>> Probe Propagation : start." << std::endl;
+  }
+  const auto& probeSurfaces = run(propProbe);
+  if (debugMode) {
+    std::cout << ">>> Probe Propagation : end." << std::endl;
+  }
+
+  if (debugMode) {
+    std::cout << ">>> Reference Propagation : start." << std::endl;
+  }
+  const auto& refSurfaces = run(propRef);
+  if (debugMode) {
+    std::cout << ">>> Reference Propagation : end." << std::endl;
+  }
+
+  // probe-ref compatibility test
+  BOOST_CHECK_EQUAL_COLLECTIONS(probeSurfaces.begin(), probeSurfaces.end(),
+                                refSurfaces.begin(), refSurfaces.end());
+}
+
+int ntests = 500;
+int skip = 0;
+bool debugMode = false;
+
+using EigenStepper = Acts::EigenStepper<>;
+using EigenPropagator = Propagator<EigenStepper, Navigator>;
+using StraightLinePropagator = Propagator<StraightLineStepper, Navigator>;
+
+EigenStepper estepper(bField);
+StraightLineStepper slstepper;
+
+EigenPropagator epropagator(estepper,
+                            Navigator({tGeometry, true, true, false},
+                                      getDefaultLogger("e_nav", Logging::INFO)),
+                            getDefaultLogger("e_prop", Logging::INFO));
+StraightLinePropagator slpropagator(slstepper,
+                                    Navigator({tGeometry, true, true, false},
+                                              getDefaultLogger("sl_nav",
+                                                               Logging::INFO)),
+                                    getDefaultLogger("sl_prop", Logging::INFO));
+
+BOOST_DATA_TEST_CASE(
+    Navigator_random,
+    bdata::random((bdata::engine = std::mt19937(), bdata::seed = 20,
+                   bdata::distribution = std::uniform_real_distribution<double>(
+                       0.5_GeV, 10_GeV))) ^
+        bdata::random((bdata::engine = std::mt19937(), bdata::seed = 21,
+                       bdata::distribution =
+                           std::uniform_real_distribution<double>(-M_PI,
+                                                                  M_PI))) ^
+        bdata::random(
+            (bdata::engine = std::mt19937(), bdata::seed = 22,
+             bdata::distribution =
+                 std::uniform_real_distribution<double>(1.0, M_PI - 1.0))) ^
+        bdata::random(
+            (bdata::engine = std::mt19937(), bdata::seed = 23,
+             bdata::distribution = std::uniform_int_distribution<int>(0, 1))) ^
+        bdata::xrange(ntests),
+    pT, phi, theta, charge, index) {
+  if (index < skip) {
+    return;
+  }
+
+  double p = pT / std::sin(theta);
+  double q = -1 + 2 * charge;
+  CurvilinearTrackParameters start(Vector4(0, 0, 0, 0), phi, theta, q / p,
+                                   std::nullopt, ParticleHypothesis::pion());
+
+  if (debugMode) {
+    std::cout << ">>> Run navigation tests with pT = " << pT
+              << "; phi = " << phi << "; theta = " << theta
+              << "; charge = " << charge << "; index = " << index << ";"
+              << std::endl;
+  }
+
+  if (debugMode) {
+    std::cout << ">>> Test self consistency epropagator" << std::endl;
+  }
+  runSelfConsistencyTest(epropagator, start, debugMode);
+  if (debugMode) {
+    std::cout << ">>> Test self consistency slpropagator" << std::endl;
+  }
+  runSelfConsistencyTest(slpropagator, start, debugMode);
+}
+
 }  // namespace Acts::Test
diff --git a/Tests/UnitTests/Core/Propagator/StraightLineStepperTests.cpp b/Tests/UnitTests/Core/Propagator/StraightLineStepperTests.cpp
index 814806086abe..65a3eb2f28b5 100644
--- a/Tests/UnitTests/Core/Propagator/StraightLineStepperTests.cpp
+++ b/Tests/UnitTests/Core/Propagator/StraightLineStepperTests.cpp
@@ -150,8 +150,6 @@ BOOST_AUTO_TEST_CASE(straight_line_stepper_test) {
   BOOST_CHECK_EQUAL(sls.charge(slsState), charge);
   BOOST_CHECK_EQUAL(sls.time(slsState), time);
 
-  //~ BOOST_CHECK_EQUAL(sls.overstepLimit(slsState), tolerance);
-
   // Step size modifies
   const std::string originalStepSize = slsState.stepSize.toString();