diff --git a/Core/CMakeLists.txt b/Core/CMakeLists.txt
index b760ba576ed..eeb6bc8b92a 100644
--- a/Core/CMakeLists.txt
+++ b/Core/CMakeLists.txt
@@ -108,6 +108,7 @@ add_subdirectory(src/Detector)
 add_subdirectory(src/Geometry)
 add_subdirectory(src/MagneticField)
 add_subdirectory(src/Material)
+add_subdirectory(src/Navigation)
 add_subdirectory(src/Propagator)
 add_subdirectory(src/Surfaces)
 add_subdirectory(src/TrackFinding)
diff --git a/Core/include/Acts/Navigation/NavigationStream.hpp b/Core/include/Acts/Navigation/NavigationStream.hpp
new file mode 100644
index 00000000000..f23fed9087c
--- /dev/null
+++ b/Core/include/Acts/Navigation/NavigationStream.hpp
@@ -0,0 +1,178 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2024 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Surfaces/BoundaryTolerance.hpp"
+#include "Acts/Utilities/Intersection.hpp"
+
+#include <tuple>
+#include <vector>
+
+namespace Acts {
+
+// To be removed when the namespace Experimental is omitted
+namespace Experimental {
+class Portal;
+}
+using namespace Experimental;
+
+class Surface;
+
+/// The NavigationStream is a container for the navigation candidates that
+/// are currentlu processed in a given context. The context could be local to a
+/// volume, or global to an entire track following.
+///
+/// The current candidates are stored in a vector of candidates, where an index
+/// is used to indicate the current active candidate.
+class NavigationStream {
+ public:
+  /// The query point for the navigation stream
+  ///
+  /// This holds the position and direction from which the navigation stream
+  /// should either be initialized or updated.
+  struct QueryPoint {
+    /// The position of the query point
+    Vector3 position = Vector3::Zero();
+    /// The direction of the query point
+    Vector3 direction = Vector3::Zero();
+  };
+
+  /// This is a candidate object of the navigation stream, it holds:
+  ///
+  /// - a Surface intersection
+  /// - a Portal : set if the surface represents a portal
+  /// - a BoundaryTolerance : the boundary tolerance used for the intersection
+  struct Candidate {
+    /// The intersection
+    ObjectIntersection<Surface> intersection =
+        ObjectIntersection<Surface>::invalid();
+    /// The portal
+    const Portal* portal = nullptr;
+    /// The boundary tolerance
+    BoundaryTolerance bTolerance = BoundaryTolerance::None();
+    /// Convenience access to surface
+    const Surface& surface() const { return *intersection.object(); }
+    /// Cinvencience access to the path length
+    ActsScalar pathLength() const { return intersection.pathLength(); }
+
+    /// Order along the path length
+    ///
+    /// @param aCandidate is the first candidate
+    /// @param bCandidate is the second candidate
+    ///
+    /// @return true if aCandidate is closer to the origin
+    constexpr static bool pathLengthOrder(const Candidate& aCandidate,
+                                          const Candidate& bCandidate) {
+      return ObjectIntersection<Surface>::pathLengthOrder(
+          aCandidate.intersection, bCandidate.intersection);
+    }
+  };
+
+  /// Switch to next next candidate
+  ///
+  /// @return true if a next candidate is available
+  bool switchToNextCandidate() {
+    if (m_currentIndex < m_candidates.size()) {
+      ++m_currentIndex;
+      return true;
+    }
+    return false;
+  }
+
+  /// Const access the current candidate
+  const Candidate& currentCandidate() const {
+    return m_candidates.at(m_currentIndex);
+  }
+
+  /// Current Index
+  std::size_t currentIndex() const { return m_currentIndex; }
+
+  /// Non-cost access the candidate vector
+  std::vector<Candidate>& candidates() { return m_candidates; }
+
+  /// Const access the candidate vector
+  const std::vector<Candidate>& candidates() const { return m_candidates; }
+
+  /// Non-cost access the current candidate
+  ///
+  /// This will throw and out of bounds exception if the stream is not
+  /// valid anymore.
+  Candidate& currentCandidate() { return m_candidates.at(m_currentIndex); }
+
+  /// The number of active candidates
+  std::size_t remainingCandidates() const {
+    return (m_candidates.size() - m_currentIndex);
+  }
+
+  /// Fill one surface into the candidate vector
+  ///
+  /// @param surface the surface to be filled
+  /// @param bTolerance the boundary tolerance used for the intersection
+  void addSurfaceCandidate(const Surface* surface,
+                           const BoundaryTolerance& bTolerance);
+
+  /// Fill n surfaces into the candidate vector
+  ///
+  /// @param surfaces the surfaces that are filled in
+  /// @param bTolerance the boundary tolerance used for the intersection
+  void addSurfaceCandidates(const std::vector<const Surface*>& surfaces,
+                            const BoundaryTolerance& bTolerance);
+
+  /// Fill one portal into the candidate vector
+  ///
+  /// @param portal the portals that are filled in
+  void addPortalCandidate(const Portal* portal);
+
+  /// Fill n portals into the candidate vector
+  ///
+  /// @param portals the portals that are filled in
+  void addPortalCandidates(const std::vector<const Portal*>& portals);
+
+  /// Initialize the stream from a query point
+  ///
+  /// @param gctx is the geometry context
+  /// @param queryPoint holds current position, direction, etc.
+  /// @param cTolerance is the candidate search tolerance
+  /// @param onSurfaceTolerance is the tolerance for on-surface intersections
+  ///
+  /// This method will first de-duplicate the candidates on basis of the surface
+  /// pointer to make sure that the multi-intersections are handled correctly.
+  /// This will allow intializeStream() to be called even as a re-initialization
+  /// and still work correctly with at one time valid candidates.
+  ///
+  /// @return true if the stream is active, false indicates that there are no valid candidates
+  bool initialize(const GeometryContext& gctx,
+                  const NavigationStream::QueryPoint& queryPoint,
+                  const BoundaryTolerance& cTolerance,
+                  ActsScalar onSurfaceTolerance = s_onSurfaceTolerance);
+
+  /// Convenience method to update a stream from a new query point,
+  /// this could be called from navigation delegates that do not require
+  /// a local state or from the navigator on the target stream
+  ///
+  /// @param gctx is the geometry context
+  /// @param queryPoint holds current position, direction, etc.
+  /// @param onSurfaceTolerance is the tolerance for on-surface intersections
+  ///
+  /// @return true if the stream is active, false indicate no valid candidates left
+  bool update(const GeometryContext& gctx,
+              const NavigationStream::QueryPoint& queryPoint,
+              ActsScalar onSurfaceTolerance = s_onSurfaceTolerance);
+
+ private:
+  /// The candidates of this navigation stream
+  std::vector<Candidate> m_candidates;
+
+  /// The currently active candidate
+  std::size_t m_currentIndex = 0u;
+};
+
+}  // namespace Acts
diff --git a/Core/include/Acts/Utilities/Intersection.hpp b/Core/include/Acts/Utilities/Intersection.hpp
index 49dc23517b7..30b986f6d24 100644
--- a/Core/include/Acts/Utilities/Intersection.hpp
+++ b/Core/include/Acts/Utilities/Intersection.hpp
@@ -69,12 +69,16 @@ class Intersection {
   /// Returns whether the intersection was successful or not
   constexpr bool isValid() const { return m_status != Status::missed; }
 
+  /// Returns the position of the interseciton
   constexpr const Position& position() const { return m_position; }
 
+  /// Returns the path length to the interseciton
   constexpr ActsScalar pathLength() const { return m_pathLength; }
 
+  /// Returns the intersection status enum
   constexpr Status status() const { return m_status; }
 
+  /// Static factory to creae an invalid instesection
   constexpr static Intersection invalid() { return Intersection(); }
 
   /// Comparison function for path length order i.e. intersection closest to
@@ -155,27 +159,35 @@ class ObjectIntersection {
   /// Returns whether the intersection was successful or not
   constexpr bool isValid() const { return m_intersection.isValid(); }
 
+  /// Returns the intersection
   constexpr const Intersection3D& intersection() const {
     return m_intersection;
   }
 
+  /// Returns the position of the interseciton
   constexpr const Intersection3D::Position& position() const {
     return m_intersection.position();
   }
 
+  /// Returns the path length to the interseciton
   constexpr ActsScalar pathLength() const {
     return m_intersection.pathLength();
   }
 
+  /// Returns the status of the interseciton
   constexpr Intersection3D::Status status() const {
     return m_intersection.status();
   }
 
+  /// Returns the object that has been intersected
   constexpr const object_t* object() const { return m_object; }
 
   constexpr std::uint8_t index() const { return m_index; }
 
-  constexpr static ObjectIntersection invalid() { return ObjectIntersection(); }
+  constexpr static ObjectIntersection invalid(
+      const object_t* object = nullptr) {
+    return ObjectIntersection(Intersection3D::invalid(), object);
+  }
 
   constexpr static bool pathLengthOrder(
       const ObjectIntersection& aIntersection,
diff --git a/Core/src/Navigation/CMakeLists.txt b/Core/src/Navigation/CMakeLists.txt
new file mode 100644
index 00000000000..81cbf616d44
--- /dev/null
+++ b/Core/src/Navigation/CMakeLists.txt
@@ -0,0 +1 @@
+target_sources(ActsCore PRIVATE NavigationStream.cpp)
diff --git a/Core/src/Navigation/NavigationStream.cpp b/Core/src/Navigation/NavigationStream.cpp
new file mode 100644
index 00000000000..3d90c5bb0d4
--- /dev/null
+++ b/Core/src/Navigation/NavigationStream.cpp
@@ -0,0 +1,153 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2024 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "Acts/Navigation/NavigationStream.hpp"
+
+#include "Acts/Detector/Portal.hpp"
+#include "Acts/Surfaces/Surface.hpp"
+
+#include <algorithm>
+
+bool Acts::NavigationStream::initialize(const GeometryContext& gctx,
+                                        const QueryPoint& queryPoint,
+                                        const BoundaryTolerance& cTolerance,
+                                        ActsScalar onSurfaceTolerance) {
+  // Position and direction from the query point
+  const Vector3& position = queryPoint.position;
+  const Vector3& direction = queryPoint.direction;
+
+  // De-duplicate first (necessary to deal correctly with multiple
+  // intersections) - sort them by surface pointer
+  std::ranges::sort(m_candidates, [](const Candidate& a, const Candidate& b) {
+    return (&a.surface()) < (&b.surface());
+  });
+  // Remove duplicates on basis of the surface pointer
+  m_candidates.erase(std::unique(m_candidates.begin(), m_candidates.end(),
+                                 [](const Candidate& a, const Candidate& b) {
+                                   return (&a.surface()) == (&b.surface());
+                                 }),
+                     m_candidates.end());
+
+  // A container collecting additional candidates from multiple
+  // valid interseciton
+  std::vector<Candidate> additionalCandidates = {};
+  for (auto& [sIntersection, portal, bTolerance] : m_candidates) {
+    // Get the surface from the object intersection
+    const Surface* surface = sIntersection.object();
+    // Intersect the surface
+    auto multiIntersection = surface->intersect(gctx, position, direction,
+                                                cTolerance, onSurfaceTolerance);
+
+    // Split them into valid intersections, keep track of potentially
+    // additional candidates
+    bool originalCandidateUpdated = false;
+    for (const auto& rsIntersection : multiIntersection.split()) {
+      // Skip negative solutions, respecting the on surface tolerance
+      if (rsIntersection.pathLength() < -onSurfaceTolerance) {
+        continue;
+      }
+      // Valid solution is either on surface or updates the distance
+      if (rsIntersection.isValid()) {
+        if (!originalCandidateUpdated) {
+          sIntersection = rsIntersection;
+          originalCandidateUpdated = true;
+        } else {
+          additionalCandidates.emplace_back(
+              Candidate{rsIntersection, portal, bTolerance});
+        }
+      }
+    }
+  }
+
+  // Append the multi intersection candidates
+  m_candidates.insert(m_candidates.end(), additionalCandidates.begin(),
+                      additionalCandidates.end());
+
+  // Sort the candidates by path length
+  std::ranges::sort(m_candidates, Candidate::pathLengthOrder);
+
+  // The we find the first invalid candidate
+  auto firstInvalid =
+      std::ranges::find_if(m_candidates, [](const Candidate& a) {
+        const auto& [aIntersection, aPortal, aTolerance] = a;
+        return !aIntersection.isValid();
+      });
+
+  // Set the range and initialize
+  m_candidates.resize(std::distance(m_candidates.begin(), firstInvalid));
+
+  m_currentIndex = 0;
+  if (m_candidates.empty()) {
+    return false;
+  }
+  return true;
+}
+
+bool Acts::NavigationStream::update(const GeometryContext& gctx,
+                                    const QueryPoint& queryPoint,
+                                    ActsScalar onSurfaceTolerance) {
+  // Position and direction from the query point
+  const Vector3& position = queryPoint.position;
+  const Vector3& direction = queryPoint.direction;
+
+  // Loop over the (currently valid) candidates and update
+  for (; m_currentIndex < m_candidates.size(); ++m_currentIndex) {
+    // Get the candidate, and resolve the tuple
+    Candidate& candidate = currentCandidate();
+    auto& [sIntersection, portal, bTolerance] = candidate;
+    // Get the surface from the object intersection
+    const Surface* surface = sIntersection.object();
+    // (re-)Intersect the surface
+    auto multiIntersection = surface->intersect(gctx, position, direction,
+                                                bTolerance, onSurfaceTolerance);
+    // Split them into valid intersections
+    for (const auto& rsIntersection : multiIntersection.split()) {
+      // Skip wrong index solution
+      if (rsIntersection.index() != sIntersection.index()) {
+        continue;
+      }
+      // Valid solution is either on surface or updates the distance
+      if (rsIntersection.isValid()) {
+        sIntersection = rsIntersection;
+        return true;
+      }
+    }
+  }
+  // No candidate was reachable
+  return false;
+}
+
+void Acts::NavigationStream::addSurfaceCandidate(
+    const Surface* surface, const BoundaryTolerance& bTolerance) {
+  m_candidates.emplace_back(Candidate{
+      ObjectIntersection<Surface>::invalid(surface), nullptr, bTolerance});
+}
+
+void Acts::NavigationStream::addSurfaceCandidates(
+    const std::vector<const Surface*>& surfaces,
+    const BoundaryTolerance& bTolerance) {
+  std::ranges::for_each(surfaces, [&](const auto* surface) {
+    m_candidates.emplace_back(Candidate{
+        ObjectIntersection<Surface>::invalid(surface), nullptr, bTolerance});
+  });
+}
+
+void Acts::NavigationStream::addPortalCandidate(const Portal* portal) {
+  m_candidates.emplace_back(
+      Candidate{ObjectIntersection<Surface>::invalid(&(portal->surface())),
+                portal, BoundaryTolerance::None()});
+}
+
+void Acts::NavigationStream::addPortalCandidates(
+    const std::vector<const Portal*>& portals) {
+  std::ranges::for_each(portals, [&](const auto& portal) {
+    m_candidates.emplace_back(
+        Candidate{ObjectIntersection<Surface>::invalid(&(portal->surface())),
+                  portal, BoundaryTolerance::None()});
+  });
+}
diff --git a/Tests/UnitTests/Core/Navigation/CMakeLists.txt b/Tests/UnitTests/Core/Navigation/CMakeLists.txt
index 69ebc582396..edca953321d 100644
--- a/Tests/UnitTests/Core/Navigation/CMakeLists.txt
+++ b/Tests/UnitTests/Core/Navigation/CMakeLists.txt
@@ -1,6 +1,7 @@
 add_unittest(PortalNavigation PortalNavigationTests.cpp)
 add_unittest(DetectorVolumeFinders DetectorVolumeFindersTests.cpp)
 add_unittest(NavigationState NavigationStateTests.cpp)
+add_unittest(NavigationStream NavigationStreamTests.cpp)
 add_unittest(NavigationStateUpdaters NavigationStateUpdatersTests.cpp)
 add_unittest(DetectorNavigator DetectorNavigatorTests.cpp)
 add_unittest(MultiWireNavigation MultiWireNavigationTests.cpp)
diff --git a/Tests/UnitTests/Core/Navigation/NavigationStreamTests.cpp b/Tests/UnitTests/Core/Navigation/NavigationStreamTests.cpp
new file mode 100644
index 00000000000..5ac91ca72bb
--- /dev/null
+++ b/Tests/UnitTests/Core/Navigation/NavigationStreamTests.cpp
@@ -0,0 +1,285 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2024 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include <boost/test/unit_test.hpp>
+
+#include "Acts/Navigation/NavigationStream.hpp"
+#include "Acts/Surfaces/CylinderSurface.hpp"
+#include "Acts/Surfaces/PlaneSurface.hpp"
+#include "Acts/Surfaces/RectangleBounds.hpp"
+#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
+
+namespace {
+
+// This creates a set of plane surfaces along the z axis
+std::vector<std::shared_ptr<Acts::Surface>> createPlaneSurfaces() {
+  auto rectangle = std::make_shared<Acts::RectangleBounds>(10., 10.);
+  // Surface A:
+  // This surface should not be reachable from (0.,0.,0.) position along z
+  Acts::Transform3 aTransform = Acts::Transform3::Identity();
+  aTransform.pretranslate(Acts::Vector3(0., 0., -20.));
+  auto surfaceA =
+      Acts::Surface::makeShared<Acts::PlaneSurface>(aTransform, rectangle);
+  // Surface B:
+  // This surface should not be reachable from (0.,0.,0.) position along z with
+  // boundary check
+  Acts::Transform3 bTransform = Acts::Transform3::Identity();
+  bTransform.pretranslate(Acts::Vector3(50., 50., 100.));
+  auto surfaceB =
+      Acts::Surface::makeShared<Acts::PlaneSurface>(bTransform, rectangle);
+  // Surface C:
+  Acts::Transform3 cTransform = Acts::Transform3::Identity();
+  cTransform.pretranslate(Acts::Vector3(0., 0., 200.));
+  auto surfaceC =
+      Acts::Surface::makeShared<Acts::PlaneSurface>(cTransform, rectangle);
+  // Surface D:
+  Acts::Transform3 dTransform = Acts::Transform3::Identity();
+  dTransform.pretranslate(Acts::Vector3(0., 0., 400.));
+  auto surfaceD =
+      Acts::Surface::makeShared<Acts::PlaneSurface>(dTransform, rectangle);
+
+  // Let's fill them shuffled
+  return {surfaceC, surfaceA, surfaceD, surfaceB};
+}
+
+// This creates a set of cylinder surfaces
+std::vector<std::shared_ptr<Acts::Surface>> createCylinders() {
+  // Surface A:
+  // A concentric cylinder with a radius of 10 and a half length of 20
+  Acts::Transform3 aTransform = Acts::Transform3::Identity();
+  auto surfaceA =
+      Acts::Surface::makeShared<Acts::CylinderSurface>(aTransform, 10., 20);
+
+  // Surface B:
+  // A  small cylinder sitting at 20, 20
+  Acts::Transform3 bTransform = Acts::Transform3::Identity();
+  bTransform.pretranslate(Acts::Vector3(20., 20., 0.));
+  auto surfaceB =
+      Acts::Surface::makeShared<Acts::CylinderSurface>(bTransform, 2., 10);
+
+  // Surface C:
+  // A concentric cylinder with a radius of 40 and a half length of 20
+  Acts::Transform3 cTransform = Acts::Transform3::Identity();
+  auto surfaceC =
+      Acts::Surface::makeShared<Acts::CylinderSurface>(cTransform, 40., 20);
+
+  // Surface C:
+  // A concentric, but shifted cylinder with a radius of 50 and a half length of
+  // 5
+  Acts::Transform3 dTransform = Acts::Transform3::Identity();
+  dTransform.pretranslate(Acts::Vector3(0., 0., 10.));
+  auto surfaceD =
+      Acts::Surface::makeShared<Acts::CylinderSurface>(dTransform, 50., 5.);
+
+  // Return in a shuffled order
+  return {surfaceC, surfaceB, surfaceA, surfaceD};
+}
+
+}  // namespace
+
+using namespace Acts;
+
+auto gContext = GeometryContext();
+
+BOOST_AUTO_TEST_SUITE(Navigation)
+
+BOOST_AUTO_TEST_CASE(NavigationStream_InitializePlanes) {
+  // Create the punch of surfaces
+  auto surfaces = createPlaneSurfaces();
+
+  NavigationStream nStreamTemplate;
+  for (const auto& surface : surfaces) {
+    nStreamTemplate.addSurfaceCandidate(surface.get(),
+                                        Acts::BoundaryTolerance::None());
+  }
+  BOOST_CHECK_EQUAL(nStreamTemplate.remainingCandidates(), 4u);
+
+  // (1) Run an initial update
+  // - from a position where all are reachable and valid
+  // - with infinite boundary tolerance
+  NavigationStream nStream = nStreamTemplate;
+  BOOST_CHECK(nStream.initialize(gContext,
+                                 {Vector3(0., 0., -30.), Vector3(0., 0., 1.)},
+                                 BoundaryTolerance::Infinite()));
+
+  BOOST_CHECK_EQUAL(nStream.remainingCandidates(), 4u);
+  BOOST_CHECK_EQUAL(&nStream.currentCandidate().surface(), surfaces[1u].get());
+
+  // (2) Run an initial update
+  // - from a position where all but one are reachable
+  // - with infinite boundary tolerance
+  nStream = nStreamTemplate;
+  BOOST_CHECK(nStream.initialize(gContext,
+                                 {Vector3(0., 0., 0.), Vector3(0., 0., 1.)},
+                                 BoundaryTolerance::Infinite()));
+  BOOST_CHECK_EQUAL(nStream.remainingCandidates(), 3u);
+  BOOST_CHECK_EQUAL(&nStream.currentCandidate().surface(), surfaces[3u].get());
+
+  // (3) Run an initial update
+  // - from a position where all would be reachable, but
+  // - with no boundary tolerance
+  nStream = nStreamTemplate;
+  BOOST_CHECK(nStream.initialize(gContext,
+                                 {Vector3(0., 0., -100.), Vector3(0., 0., 1.)},
+                                 BoundaryTolerance::None()));
+  BOOST_CHECK_EQUAL(nStream.remainingCandidates(), 3u);
+
+  // (4) Run an initial update
+  // - none of the surfaces should be reachable
+  nStream = nStreamTemplate;
+  BOOST_CHECK(!nStream.initialize(gContext,
+                                  {Vector3(0., 0., 0.), Vector3(1., 0., 0.)},
+                                  BoundaryTolerance::Infinite()));
+  BOOST_CHECK_EQUAL(nStream.remainingCandidates(), 0u);
+  BOOST_CHECK_THROW(nStream.currentCandidate(), std::out_of_range);
+
+  // (5) Test de-duplication
+  nStream = nStreamTemplate;
+  nStreamTemplate.addSurfaceCandidate(surfaces[0].get(),
+                                      Acts::BoundaryTolerance::None());
+  // One surface is duplicated in the stream
+  BOOST_CHECK_EQUAL(nStreamTemplate.remainingCandidates(), 5u);
+  // Initialize stream reaches all surfaces, but also de-duplicates
+  BOOST_CHECK(nStream.initialize(gContext,
+                                 {Vector3(0., 0., -100.), Vector3(0., 0., 1.)},
+                                 BoundaryTolerance::Infinite()));
+  BOOST_CHECK_EQUAL(nStream.remainingCandidates(), 4u);
+}
+
+BOOST_AUTO_TEST_CASE(NavigationStream_UpdatePlanes) {
+  // Create the punch of surfaces
+  auto surfaces = createPlaneSurfaces();
+
+  // Surfaces are filled with no boundary tolerance, we require them to be
+  // reachable and intersections inside bounds
+  NavigationStream nStreamTemplate;
+  for (const auto& surface : surfaces) {
+    nStreamTemplate.addSurfaceCandidate(surface.get(),
+                                        Acts::BoundaryTolerance::None());
+  }
+  BOOST_CHECK_EQUAL(nStreamTemplate.remainingCandidates(), 4u);
+
+  // Run an initial update
+  // - from a position where all are reachable and valid
+  // - with infinite boundary tolerance
+  NavigationStream::QueryPoint qPoint = {Vector3(0., 0., -30.),
+                                         Vector3(0., 0., 1.)};
+
+  NavigationStream nStream = nStreamTemplate;
+  BOOST_CHECK(
+      nStream.initialize(gContext, qPoint, BoundaryTolerance::Infinite()));
+  BOOST_CHECK_EQUAL(nStream.remainingCandidates(), 4u);
+  BOOST_CHECK_EQUAL(&nStream.currentCandidate().surface(), surfaces[1u].get());
+  CHECK_CLOSE_ABS(nStream.currentCandidate().pathLength(), 10.,
+                  std::numeric_limits<ActsScalar>::epsilon());
+
+  // Let's push a bit closer to the surface
+  qPoint.position = Vector3(0., 0., -22.);
+  BOOST_CHECK(nStream.update(gContext, qPoint));
+  // Surface unchanged, but the intersection should be closer
+  BOOST_CHECK_EQUAL(&nStream.currentCandidate().surface(), surfaces[1u].get());
+  CHECK_CLOSE_ABS(nStream.currentCandidate().pathLength(), 2.,
+                  std::numeric_limits<ActsScalar>::epsilon());
+
+  // Uuuups, an overstep
+  qPoint.position = Vector3(0., 0., -19.5);
+  BOOST_CHECK(nStream.update(gContext, qPoint));
+  // Surface still unchanged, but pathLength is now negative
+  BOOST_CHECK_EQUAL(&nStream.currentCandidate().surface(), surfaces[1u].get());
+  CHECK_CLOSE_ABS(nStream.currentCandidate().pathLength(), -0.5,
+                  std::numeric_limits<ActsScalar>::epsilon());
+
+  // Finally hit it
+  qPoint.position = Vector3(0., 0., -20.);
+  BOOST_CHECK(nStream.update(gContext, qPoint));
+  // Surface still unchanged, however, now withL
+  // - pathlength smaller on surface tolerance, intersection status onSurface
+  BOOST_CHECK_EQUAL(&nStream.currentCandidate().surface(), surfaces[1u].get());
+  CHECK_CLOSE_ABS(
+      nStream.currentCandidate().pathLength(), s_onSurfaceTolerance,
+      std::numeric_limits<ActsScalar>::epsilon() + s_onSurfaceTolerance);
+  BOOST_CHECK_EQUAL(nStream.currentCandidate().intersection.status(),
+                    IntersectionStatus::onSurface);
+  // Let's say the stepper confirms this
+  BOOST_CHECK(nStream.switchToNextCandidate());
+  // Surface is now surfaceB
+  BOOST_CHECK_EQUAL(&nStream.currentCandidate().surface(), surfaces[3u].get());
+  // Distance should be the initial estimate from the intialializeStream() call
+  CHECK_CLOSE_ABS(nStream.currentCandidate().pathLength(), 130.,
+                  std::numeric_limits<ActsScalar>::epsilon());
+  // Query update will re-evaluate this one: however, we will miss the surface
+  // due to outside bounds - and will switch to the next candidate: which sits
+  // at 200 and then will yield 220
+  BOOST_CHECK(nStream.update(gContext, qPoint));
+  CHECK_CLOSE_ABS(nStream.currentCandidate().pathLength(), 220.,
+                  std::numeric_limits<ActsScalar>::epsilon());
+  // Oh noooo, an actor just kicked in and changed the direction
+  qPoint.direction = Vector3(0., 1., 1.).normalized();
+  // All is lost, no surface is reachable anymore
+  BOOST_CHECK(!nStream.update(gContext, qPoint));
+}
+
+BOOST_AUTO_TEST_CASE(NavigationStream_InitializeCylinders) {
+  // Create the cylinder setup
+  auto surfaces = createCylinders();
+
+  // Let us fill the surfaces into the navigation stream
+  NavigationStream nStreamTemplate;
+  for (const auto& surface : surfaces) {
+    nStreamTemplate.addSurfaceCandidates({surface.get()},
+                                         Acts::BoundaryTolerance::None());
+  }
+  BOOST_CHECK_EQUAL(nStreamTemplate.remainingCandidates(), 4u);
+
+  // (1) Run an initial update - from a position/direction where all are
+  // reachable
+  // - with infinite boundary tolerance
+  NavigationStream nStream = nStreamTemplate;
+  BOOST_CHECK(nStream.initialize(
+      gContext, {Vector3(0., 0., 0.), Vector3(1., 1., 0.).normalized()},
+      BoundaryTolerance::Infinite()));
+  // We should have 5 candidates, as one cylinder is reachable twice
+  BOOST_CHECK_EQUAL(nStream.remainingCandidates(), 5u);
+  // First one is inner candidate
+  BOOST_CHECK_EQUAL(&nStream.currentCandidate().surface(), surfaces[2].get());
+  // Surface of 2nd and 3rd candidate should be the same
+  BOOST_CHECK_EQUAL(&nStream.candidates()[1u].surface(), surfaces[1].get());
+  BOOST_CHECK_EQUAL(&nStream.candidates()[2u].surface(), surfaces[1].get());
+
+  // (2) Run an initial update - from a position/direction where only
+  // the concentric ones are reachable
+  // - with infinite boundary tolerance
+  nStream = nStreamTemplate;
+  BOOST_CHECK(nStream.initialize(gContext,
+                                 {Vector3(0., 0., 0.), Vector3(1., 0., 0.)},
+                                 BoundaryTolerance::Infinite()));
+  // We should have 3 candidates
+  BOOST_CHECK_EQUAL(nStream.remainingCandidates(), 3u);
+
+  // (3) Run an initial update - from a position/direction where only the
+  // concentric ones within bounds are reachable
+  nStream = nStreamTemplate;
+  BOOST_CHECK(nStream.initialize(gContext,
+                                 {Vector3(0., 0., 0.), Vector3(1., 0., 0.)},
+                                 BoundaryTolerance::None()));
+  // We should have 2 candidates
+  BOOST_CHECK_EQUAL(nStream.remainingCandidates(), 2u);
+
+  // (4) Run an initial update - from a position/direction where none are
+  // reachable
+  // - (even) with infinite boundary tolerance
+  nStream = nStreamTemplate;
+  BOOST_CHECK(!nStream.initialize(gContext,
+                                  {Vector3(0., 0., 0.), Vector3(0., 0., 1.)},
+                                  BoundaryTolerance::None()));
+  // We should have 0 candidates
+  BOOST_CHECK_EQUAL(nStream.remainingCandidates(), 0u);
+  BOOST_CHECK_THROW(nStream.currentCandidate(), std::out_of_range);
+}
+
+BOOST_AUTO_TEST_SUITE_END()