diff --git a/cpp/dolfinx/geometry/BoundingBoxTree.cpp b/cpp/dolfinx/geometry/BoundingBoxTree.cpp
index f692c58b1e5..76493cf4e77 100644
--- a/cpp/dolfinx/geometry/BoundingBoxTree.cpp
+++ b/cpp/dolfinx/geometry/BoundingBoxTree.cpp
@@ -1,5 +1,5 @@
-// Copyright (C) 2013-2021 Chris N. Richardson, Anders Logg, Garth N. Wells and
-// Jørgen S. Dokken
+// Copyright (C) 2013-2022 Chris N. Richardson, Anders Logg, Garth N. Wells,
+// Jørgen S. Dokken, Sarah Roggendorf
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -33,9 +33,10 @@ std::vector<std::int32_t> range(const mesh::Mesh& mesh, int tdim)
   return r;
 }
 //-----------------------------------------------------------------------------
-// Compute bounding box of mesh entity
-std::array<std::array<double, 3>, 2>
-compute_bbox_of_entity(const mesh::Mesh& mesh, int dim, std::int32_t index)
+// Compute bounding box of mesh entity. The bounding box is defined by (lower
+// left corner, top right corner). Storage flattened row-major
+std::array<double, 6> compute_bbox_of_entity(const mesh::Mesh& mesh, int dim,
+                                             std::int32_t index)
 {
   // Get the geometrical indices for the mesh entity
   std::span<const double> xg = mesh.geometry().x();
@@ -45,40 +46,41 @@ compute_bbox_of_entity(const mesh::Mesh& mesh, int dim, std::int32_t index)
   const std::vector<std::int32_t> vertex_indices
       = mesh::entities_to_geometry(mesh, dim, entity, false);
 
-  std::array<std::array<double, 3>, 2> b;
-  b[0] = {xg[3 * vertex_indices.front()], xg[3 * vertex_indices.front() + 1],
-          xg[3 * vertex_indices.front() + 2]};
-  b[1] = b[0];
+  std::array<double, 6> b;
+  auto b0 = std::span(b).subspan<0, 3>();
+  auto b1 = std::span(b).subspan<3, 3>();
+  std::copy_n(std::next(xg.begin(), 3 * vertex_indices.front()), 3, b0.begin());
+  std::copy_n(std::next(xg.begin(), 3 * vertex_indices.front()), 3, b1.begin());
 
   // Compute min and max over vertices
   for (int local_vertex : vertex_indices)
   {
     for (int j = 0; j < 3; ++j)
     {
-      b[0][j] = std::min(b[0][j], xg[3 * local_vertex + j]);
-      b[1][j] = std::max(b[1][j], xg[3 * local_vertex + j]);
+      b0[j] = std::min(b0[j], xg[3 * local_vertex + j]);
+      b1[j] = std::max(b1[j], xg[3 * local_vertex + j]);
     }
   }
 
   return b;
 }
 //-----------------------------------------------------------------------------
-// Compute bounding box of bounding boxes
-std::array<std::array<double, 3>, 2> compute_bbox_of_bboxes(
-    const std::span<const std::pair<std::array<std::array<double, 3>, 2>,
-                                    std::int32_t>>& leaf_bboxes)
+// Compute bounding box of bounding boxes. Each bounding box is defined as a
+// tuple (corners, entity_index). The corners of the bounding box is flattened
+// row-major as (lower left corner, top right corner).
+std::array<double, 6> compute_bbox_of_bboxes(
+    std::span<const std::pair<std::array<double, 6>, std::int32_t>> leaf_bboxes)
 {
   // Compute min and max over remaining boxes
-  std::array<std::array<double, 3>, 2> b;
-  b[0] = leaf_bboxes[0].first[0];
-  b[1] = leaf_bboxes[0].first[1];
-  for (auto& box : leaf_bboxes)
+  std::array<double, 6> b = leaf_bboxes.front().first;
+
+  for (auto [box, _] : leaf_bboxes)
   {
-    std::transform(box.first[0].begin(), box.first[0].end(), b[0].begin(),
-                   b[0].begin(),
+    std::transform(box.cbegin(), std::next(box.cbegin(), 3), b.cbegin(),
+                   b.begin(),
                    [](double a, double b) { return std::min(a, b); });
-    std::transform(box.first[1].begin(), box.first[1].end(), b[1].begin(),
-                   b[1].begin(),
+    std::transform(std::next(box.cbegin(), 3), box.cend(),
+                   std::next(b.cbegin(), 3), std::next(b.begin(), 3),
                    [](double a, double b) { return std::max(a, b); });
   }
 
@@ -86,38 +88,31 @@ std::array<std::array<double, 3>, 2> compute_bbox_of_bboxes(
 }
 //------------------------------------------------------------------------------
 int _build_from_leaf(
-    std::span<std::pair<std::array<std::array<double, 3>, 2>, std::int32_t>>
-        leaf_bboxes,
-    std::vector<std::array<int, 2>>& bboxes,
-    std::vector<double>& bbox_coordinates)
+    std::span<std::pair<std::array<double, 6>, std::int32_t>> leaf_bboxes,
+    std::vector<int>& bboxes, std::vector<double>& bbox_coordinates)
 {
   if (leaf_bboxes.size() == 1)
   {
     // Reached leaf
 
     // Get bounding box coordinates for leaf
-    const std::int32_t entity_index = leaf_bboxes[0].second;
-    const std::array<double, 3> b0 = leaf_bboxes[0].first[0];
-    const std::array<double, 3> b1 = leaf_bboxes[0].first[1];
+    const auto [b, entity_index] = leaf_bboxes.front();
 
     // Store bounding box data
-    bboxes.push_back({entity_index, entity_index});
-    bbox_coordinates.insert(bbox_coordinates.end(), b0.begin(), b0.end());
-    bbox_coordinates.insert(bbox_coordinates.end(), b1.begin(), b1.end());
-    return bboxes.size() - 1;
+    bboxes.push_back(entity_index);
+    bboxes.push_back(entity_index);
+    std::copy_n(b.begin(), 6, std::back_inserter(bbox_coordinates));
+    return bboxes.size() / 2 - 1;
   }
   else
   {
     // Compute bounding box of all bounding boxes
-    std::array<std::array<double, 3>, 2> b
-        = compute_bbox_of_bboxes(leaf_bboxes);
-    std::array<double, 3> b0 = b[0];
-    std::array<double, 3> b1 = b[1];
+    std::array<double, 6> b = compute_bbox_of_bboxes(leaf_bboxes);
 
     // Sort bounding boxes along longest axis
     std::array<double, 3> b_diff;
-    std::transform(b1.begin(), b1.end(), b0.begin(), b_diff.begin(),
-                   std::minus<double>());
+    std::transform(std::next(b.cbegin(), 3), b.cend(), b.cbegin(),
+                   b_diff.begin(), std::minus<double>());
     const std::size_t axis = std::distance(
         b_diff.begin(), std::max_element(b_diff.begin(), b_diff.end()));
 
@@ -125,49 +120,39 @@ int _build_from_leaf(
     std::nth_element(leaf_bboxes.begin(), middle, leaf_bboxes.end(),
                      [axis](auto& p0, auto& p1) -> bool
                      {
-                       double x0 = p0.first[0][axis] + p0.first[1][axis];
-                       double x1 = p1.first[0][axis] + p1.first[1][axis];
+                       double x0 = p0.first[axis] + p0.first[3 + axis];
+                       double x1 = p1.first[axis] + p1.first[3 + axis];
                        return x0 < x1;
                      });
 
     // Split bounding boxes into two groups and call recursively
     assert(!leaf_bboxes.empty());
     std::size_t part = leaf_bboxes.size() / 2;
-    std::array bbox{
-        _build_from_leaf(leaf_bboxes.first(part), bboxes, bbox_coordinates),
-        _build_from_leaf(leaf_bboxes.last(leaf_bboxes.size() - part), bboxes,
-                         bbox_coordinates)};
+    int bbox0
+        = _build_from_leaf(leaf_bboxes.first(part), bboxes, bbox_coordinates);
+    int bbox1 = _build_from_leaf(leaf_bboxes.last(leaf_bboxes.size() - part),
+                                 bboxes, bbox_coordinates);
 
     // Store bounding box data. Note that root box will be added last.
-    bboxes.push_back(bbox);
-    bbox_coordinates.insert(bbox_coordinates.end(), b0.begin(), b0.end());
-    bbox_coordinates.insert(bbox_coordinates.end(), b1.begin(), b1.end());
-    return bboxes.size() - 1;
+    bboxes.push_back(bbox0);
+    bboxes.push_back(bbox1);
+    std::copy_n(b.begin(), 6, std::back_inserter(bbox_coordinates));
+    return bboxes.size() / 2 - 1;
   }
 }
 //-----------------------------------------------------------------------------
 std::pair<std::vector<std::int32_t>, std::vector<double>> build_from_leaf(
-    std::vector<std::pair<std::array<std::array<double, 3>, 2>, std::int32_t>>
-        leaf_bboxes)
+    std::vector<std::pair<std::array<double, 6>, std::int32_t>> leaf_bboxes)
 {
-  std::vector<std::array<std::int32_t, 2>> bboxes;
+  std::vector<std::int32_t> bboxes;
   std::vector<double> bbox_coordinates;
   _build_from_leaf(leaf_bboxes, bboxes, bbox_coordinates);
-
-  std::vector<std::int32_t> bbox_array(2 * bboxes.size());
-  for (std::size_t i = 0; i < bboxes.size(); ++i)
-  {
-    bbox_array[2 * i] = bboxes[i][0];
-    bbox_array[2 * i + 1] = bboxes[i][1];
-  }
-
-  return {std::move(bbox_array), std::move(bbox_coordinates)};
+  return {std::move(bboxes), std::move(bbox_coordinates)};
 }
 //-----------------------------------------------------------------------------
 int _build_from_point(
     std::span<std::pair<std::array<double, 3>, std::int32_t>> points,
-    std::vector<std::array<std::int32_t, 2>>& bboxes,
-    std::vector<double>& bbox_coordinates)
+    std::vector<std::int32_t>& bboxes, std::vector<double>& bbox_coordinates)
 {
   // Reached leaf
   if (points.size() == 1)
@@ -176,12 +161,13 @@ int _build_from_point(
 
     // Index of entity contained in leaf
     const std::int32_t c1 = points[0].second;
-    bboxes.push_back({c1, c1});
+    bboxes.push_back(c1);
+    bboxes.push_back(c1);
     bbox_coordinates.insert(bbox_coordinates.end(), points[0].first.begin(),
                             points[0].first.end());
     bbox_coordinates.insert(bbox_coordinates.end(), points[0].first.begin(),
                             points[0].first.end());
-    return bboxes.size() - 1;
+    return bboxes.size() / 2 - 1;
   }
 
   // Compute bounding box of all points
@@ -204,16 +190,17 @@ int _build_from_point(
   // Split bounding boxes into two groups and call recursively
   assert(!points.empty());
   std::size_t part = points.size() / 2;
-  std::array bbox{
-      _build_from_point(points.first(part), bboxes, bbox_coordinates),
-      _build_from_point(points.last(points.size() - part), bboxes,
-                        bbox_coordinates)};
+  std::int32_t bbox0
+      = _build_from_point(points.first(part), bboxes, bbox_coordinates);
+  std::int32_t bbox1 = _build_from_point(points.last(points.size() - part),
+                                         bboxes, bbox_coordinates);
 
   // Store bounding box data. Note that root box will be added last.
-  bboxes.push_back(bbox);
+  bboxes.push_back(bbox0);
+  bboxes.push_back(bbox1);
   bbox_coordinates.insert(bbox_coordinates.end(), b0.begin(), b0.end());
   bbox_coordinates.insert(bbox_coordinates.end(), b1.begin(), b1.end());
-  return bboxes.size() - 1;
+  return bboxes.size() / 2 - 1;
 }
 //-----------------------------------------------------------------------------
 } // namespace
@@ -242,16 +229,14 @@ BoundingBoxTree::BoundingBoxTree(const mesh::Mesh& mesh, int tdim,
   mesh.topology_mutable().create_connectivity(tdim, mesh.topology().dim());
 
   // Create bounding boxes for all mesh entities (leaves)
-  std::vector<std::pair<std::array<std::array<double, 3>, 2>, std::int32_t>>
-      leaf_bboxes;
+  std::vector<std::pair<std::array<double, 6>, std::int32_t>> leaf_bboxes;
   leaf_bboxes.reserve(entities.size());
   for (std::int32_t e : entities)
   {
-    std::array<std::array<double, 3>, 2> b
-        = compute_bbox_of_entity(mesh, tdim, e);
-    std::transform(b[0].begin(), b[0].end(), b[0].begin(),
+    std::array<double, 6> b = compute_bbox_of_entity(mesh, tdim, e);
+    std::transform(b.cbegin(), std::next(b.cbegin(), 3), b.begin(),
                    [padding](double x) { return x - padding; });
-    std::transform(b[1].begin(), b[1].end(), b[1].begin(),
+    std::transform(std::next(b.begin(), 3), b.end(), std::next(b.begin(), 3),
                    [padding](double& x) { return x + padding; });
     leaf_bboxes.emplace_back(b, e);
   }
@@ -269,16 +254,10 @@ BoundingBoxTree::BoundingBoxTree(
     : _tdim(0)
 {
   // Recursively build the bounding box tree from the leaves
-  std::vector<std::array<int, 2>> bboxes;
   if (!points.empty())
   {
-    _build_from_point(std::span(points), bboxes, _bbox_coordinates);
-    _bboxes.resize(2 * bboxes.size());
-    for (std::size_t i = 0; i < bboxes.size(); ++i)
-    {
-      _bboxes[2 * i] = bboxes[i][0];
-      _bboxes[2 * i + 1] = bboxes[i][1];
-    }
+    _bboxes.clear();
+    _build_from_point(std::span(points), _bboxes, _bbox_coordinates);
   }
 
   LOG(INFO) << "Computed bounding box tree with " << num_bboxes()
@@ -305,14 +284,13 @@ BoundingBoxTree BoundingBoxTree::create_global_tree(MPI_Comm comm) const
   MPI_Allgather(send_bbox.data(), 6, MPI_DOUBLE, recv_bbox.data(), 6,
                 MPI_DOUBLE, comm);
 
-  std::vector<std::pair<std::array<std::array<double, 3>, 2>, std::int32_t>>
-      _recv_bbox(mpi_size);
+  std::vector<std::pair<std::array<double, 6>, std::int32_t>> _recv_bbox(
+      mpi_size);
   for (std::size_t i = 0; i < _recv_bbox.size(); ++i)
   {
-    std::copy_n(std::next(recv_bbox.begin(), 6 * i), 3,
-                _recv_bbox[i].first[0].begin());
-    std::copy_n(std::next(recv_bbox.begin(), 6 * i + 3), 3,
-                _recv_bbox[i].first[1].begin());
+    std::copy_n(std::next(recv_bbox.begin(), 6 * i), 6,
+                _recv_bbox[i].first.begin());
+
     _recv_bbox[i].second = i;
   }
 
@@ -362,13 +340,11 @@ void BoundingBoxTree::tree_print(std::stringstream& s, int i) const
   }
 }
 //-----------------------------------------------------------------------------
-std::array<std::array<double, 3>, 2>
-BoundingBoxTree::get_bbox(std::size_t node) const
+std::array<double, 6> BoundingBoxTree::get_bbox(std::size_t node) const
 {
-  std::array<std::array<double, 3>, 2> x;
-  std::copy_n(std::next(_bbox_coordinates.begin(), 6 * node), 3, x[0].begin());
-  std::copy_n(std::next(_bbox_coordinates.begin(), 6 * node + 3), 3,
-              x[1].begin());
+  std::array<double, 6> x;
+  std::copy_n(std::next(_bbox_coordinates.begin(), 6 * node), 6, x.begin());
+
   return x;
 }
 //-----------------------------------------------------------------------------
diff --git a/cpp/dolfinx/geometry/BoundingBoxTree.h b/cpp/dolfinx/geometry/BoundingBoxTree.h
index c694d7af590..26a3d3543d9 100644
--- a/cpp/dolfinx/geometry/BoundingBoxTree.h
+++ b/cpp/dolfinx/geometry/BoundingBoxTree.h
@@ -70,9 +70,8 @@ class BoundingBoxTree
 
   /// Return bounding box coordinates for a given node in the tree
   /// @param[in] node The bounding box node index
-  /// @return The bounding box where [0] is the lower corner and [1] is
-  /// the upper corner
-  std::array<std::array<double, 3>, 2> get_bbox(std::size_t node) const;
+  /// @return A copy of bounding box (lower_corner, upper_corner). Shape (2,3).
+  std::array<double, 6> get_bbox(std::size_t node) const;
 
   /// Compute a global bounding tree (collective on comm)
   /// This can be used to find which process a point might have a
diff --git a/cpp/dolfinx/geometry/gjk.cpp b/cpp/dolfinx/geometry/gjk.cpp
index 0a64e17f225..612984278ad 100644
--- a/cpp/dolfinx/geometry/gjk.cpp
+++ b/cpp/dolfinx/geometry/gjk.cpp
@@ -27,8 +27,8 @@ nearest_simplex(const std::span<const double>& s)
   case 2:
   {
     // Compute lm = dot(s0, ds / |ds|)
-    auto s0 = s.subspan(0, 3);
-    auto s1 = s.subspan(3, 3);
+    auto s0 = s.subspan<0, 3>();
+    auto s1 = s.subspan<3, 3>();
     std::array ds = {s1[0] - s0[0], s1[1] - s0[1], s1[2] - s0[2]};
     const double lm = -(s0[0] * ds[0] + s0[1] * ds[1] + s0[2] * ds[2])
                       / (ds[0] * ds[0] + ds[1] * ds[1] + ds[2] * ds[2]);
@@ -48,9 +48,9 @@ nearest_simplex(const std::span<const double>& s)
   }
   case 3:
   {
-    auto a = s.subspan(0, 3);
-    auto b = s.subspan(3, 3);
-    auto c = s.subspan(6, 3);
+    auto a = s.subspan<0, 3>();
+    auto b = s.subspan<3, 3>();
+    auto c = s.subspan<6, 3>();
     auto length = [](auto& x, auto& y)
     {
       return std::transform_reduce(
@@ -155,10 +155,10 @@ nearest_simplex(const std::span<const double>& s)
   }
   case 4:
   {
-    auto s0 = s.subspan(0, 3);
-    auto s1 = s.subspan(3, 3);
-    auto s2 = s.subspan(6, 3);
-    auto s3 = s.subspan(9, 3);
+    auto s0 = s.subspan<0, 3>();
+    auto s1 = s.subspan<3, 3>();
+    auto s2 = s.subspan<6, 3>();
+    auto s3 = s.subspan<9, 3>();
     auto W1 = math::cross(s0, s1);
     auto W2 = math::cross(s2, s3);
 
@@ -178,7 +178,7 @@ nearest_simplex(const std::span<const double>& s)
       if (f_inside[0]) // The origin is inside the tetrahedron
         return {std::vector<double>(s.begin(), s.end()), {0, 0, 0}};
       else // The origin projection P faces BCD
-        return nearest_simplex(s.subspan(0, 3 * 3));
+        return nearest_simplex(s.subspan<0, 3 * 3>());
     }
 
     // Test ACD, ABD and/or ABC
diff --git a/cpp/dolfinx/geometry/utils.cpp b/cpp/dolfinx/geometry/utils.cpp
index f26a6361696..fe482656337 100644
--- a/cpp/dolfinx/geometry/utils.cpp
+++ b/cpp/dolfinx/geometry/utils.cpp
@@ -1,5 +1,5 @@
-// Copyright (C) 2006-2021 Chris N. Richardson, Anders Logg, Garth N. Wells and
-// Jørgen S. Dokken
+// Copyright (C) 2006-2022 Chris N. Richardson, Anders Logg, Garth N. Wells,
+// Jørgen S. Dokken, Sarah Roggendorf
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -20,7 +20,7 @@ namespace
 {
 //-----------------------------------------------------------------------------
 // Check whether bounding box is a leaf node
-constexpr bool is_leaf(const std::array<int, 2>& bbox)
+constexpr bool is_leaf(std::span<const int, 2> bbox)
 {
   // Leaf nodes are marked by setting child_0 equal to child_1
   return bbox[0] == bbox[1];
@@ -29,16 +29,19 @@ constexpr bool is_leaf(const std::array<int, 2>& bbox)
 /// A point `x` is inside a bounding box `b` if each component of its
 /// coordinates lies within the range `[b(0,i), b(1,i)]` that defines the bounds
 /// of the bounding box, b(0,i) <= x[i] <= b(1,i) for i = 0, 1, 2
-constexpr bool point_in_bbox(const std::array<std::array<double, 3>, 2>& b,
-                             const std::array<double, 3>& x)
+constexpr bool point_in_bbox(std::span<const double, 6> b,
+                             std::span<const double, 3> x)
 {
+  assert(b.size() == 6);
   constexpr double rtol = 1e-14;
   bool in = true;
+  auto b0 = b.subspan<0, 3>();
+  auto b1 = b.subspan<3, 3>();
   for (int i = 0; i < 3; i++)
   {
-    double eps = rtol * (b[1][i] - b[0][i]);
-    in &= x[i] >= (b[0][i] - eps);
-    in &= x[i] <= (b[1][i] + eps);
+    double eps = rtol * (b1[i] - b0[i]);
+    in &= x[i] >= (b0[i] - eps);
+    in &= x[i] <= (b1[i] + eps);
   }
 
   return in;
@@ -47,16 +50,22 @@ constexpr bool point_in_bbox(const std::array<std::array<double, 3>, 2>& b,
 /// A bounding box "a" is contained inside another bounding box "b", if each
 /// of its intervals [a(0,i), a(1,i)] is contained in [b(0,i), b(1,i)],
 /// a(0,i) <= b(1, i) and a(1,i) >= b(0, i)
-constexpr bool bbox_in_bbox(const std::array<std::array<double, 3>, 2>& a,
-                            const std::array<std::array<double, 3>, 2>& b)
+constexpr bool bbox_in_bbox(std::span<const double> a,
+                            std::span<const double> b)
 {
+  assert(b.size() == 6);
+  assert(a.size() == 6);
   constexpr double rtol = 1e-14;
   bool in = true;
+  auto a0 = a.subspan<0, 3>();
+  auto a1 = a.subspan<3, 3>();
+  auto b0 = b.subspan<0, 3>();
+  auto b1 = b.subspan<3, 3>();
   for (int i = 0; i < 3; i++)
   {
-    double eps = rtol * (b[1][i] - b[0][i]);
-    in &= a[1][i] >= (b[0][i] - eps);
-    in &= a[0][i] <= (b[1][i] + eps);
+    double eps = rtol * (b1[i] - b0[i]);
+    in &= a1[i] >= (b0[i] - eps);
+    in &= a0[i] <= (b1[i] + eps);
   }
 
   return in;
@@ -64,19 +73,21 @@ constexpr bool bbox_in_bbox(const std::array<std::array<double, 3>, 2>& a,
 //-----------------------------------------------------------------------------
 // Compute closest entity {closest_entity, R2} (recursive)
 std::pair<std::int32_t, double> _compute_closest_entity(
-    const geometry::BoundingBoxTree& tree, const std::array<double, 3>& point,
+    const geometry::BoundingBoxTree& tree, std::span<const double, 3> point,
     int node, const mesh::Mesh& mesh, std::int32_t closest_entity, double R2)
 {
   // Get children of current bounding box node (child_1 denotes entity
   // index for leaves)
-  const std::array bbox = tree.bbox(node);
+  const std::array<int, 2> bbox = tree.bbox(node);
   double r2;
   if (is_leaf(bbox))
   {
     // If point cloud tree the exact distance is easy to compute
     if (tree.tdim() == 0)
     {
-      std::array<double, 3> diff = tree.get_bbox(node)[0];
+      std::array<double, 6> bbox_coord = tree.get_bbox(node);
+      auto diff = std::span(bbox_coord).subspan<0, 3>();
+
       for (std::size_t k = 0; k < 3; ++k)
         diff[k] -= point[k];
       r2 = diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2];
@@ -88,9 +99,9 @@ std::pair<std::int32_t, double> _compute_closest_entity(
       // obtain exact distance to the convex hull of the entity
       if (r2 <= R2)
       {
-        r2 = geometry::squared_distance(mesh, tree.tdim(),
-                                        std::span(&bbox[1], 1),
-                                        {{point[0], point[1], point[2]}})
+        r2 = geometry::squared_distance(
+                 mesh, tree.tdim(), std::span(std::next(bbox.begin(), 1), 1),
+                 point)
                  .front();
       }
     }
@@ -127,7 +138,7 @@ std::pair<std::int32_t, double> _compute_closest_entity(
 /// @param[in] points The points (shape=(num_points, 3))
 /// @param[in, out] entities The list of colliding entities (local to process)
 void _compute_collisions_point(const geometry::BoundingBoxTree& tree,
-                               const std::array<double, 3>& p,
+                               std::span<const double, 3> p,
                                std::vector<int>& entities)
 {
   std::deque<std::int32_t> stack;
@@ -135,7 +146,7 @@ void _compute_collisions_point(const geometry::BoundingBoxTree& tree,
 
   while (next != -1)
   {
-    std::array bbox = tree.bbox(next);
+    const std::array<int, 2> bbox = tree.bbox(next);
     next = -1;
 
     if (is_leaf(bbox))
@@ -180,16 +191,15 @@ void _compute_collisions_point(const geometry::BoundingBoxTree& tree,
 // Compute collisions with tree (recursive)
 void _compute_collisions_tree(const geometry::BoundingBoxTree& A,
                               const geometry::BoundingBoxTree& B, int node_A,
-                              int node_B,
-                              std::vector<std::array<int, 2>>& entities)
+                              int node_B, std::vector<int>& entities)
 {
   // If bounding boxes don't collide, then don't search further
   if (!bbox_in_bbox(A.get_bbox(node_A), B.get_bbox(node_B)))
     return;
 
   // Get bounding boxes for current nodes
-  const std::array bbox_A = A.bbox(node_A);
-  const std::array bbox_B = B.bbox(node_B);
+  const std::array<int, 2> bbox_A = A.bbox(node_A);
+  const std::array<int, 2> bbox_B = B.bbox(node_B);
 
   // Check whether we've reached a leaf in A or B
   const bool is_leaf_A = is_leaf(bbox_A);
@@ -198,7 +208,8 @@ void _compute_collisions_tree(const geometry::BoundingBoxTree& A,
   {
     // If both boxes are leaves (which we know collide), then add them
     // child_1 denotes entity for leaves
-    entities.push_back({bbox_A[1], bbox_B[1]});
+    entities.push_back(bbox_A[1]);
+    entities.push_back(bbox_B[1]);
   }
   else if (is_leaf_A)
   {
@@ -257,12 +268,11 @@ geometry::create_midpoint_tree(const mesh::Mesh& mesh, int tdim,
   return geometry::BoundingBoxTree(points);
 }
 //-----------------------------------------------------------------------------
-std::vector<std::array<int, 2>>
-geometry::compute_collisions(const BoundingBoxTree& tree0,
-                             const BoundingBoxTree& tree1)
+std::vector<int> geometry::compute_collisions(const BoundingBoxTree& tree0,
+                                              const BoundingBoxTree& tree1)
 {
   // Call recursive find function
-  std::vector<std::array<int, 2>> entities;
+  std::vector<int> entities;
   if (tree0.num_bboxes() > 0 and tree1.num_bboxes() > 0)
   {
     _compute_collisions_tree(tree0, tree1, tree0.num_bboxes() - 1,
@@ -283,8 +293,7 @@ geometry::compute_collisions(const BoundingBoxTree& tree,
     for (std::size_t p = 0; p < points.size() / 3; ++p)
     {
       _compute_collisions_point(
-          tree, {points[3 * p + 0], points[3 * p + 1], points[3 * p + 2]},
-          entities);
+          tree, std::span<const double, 3>(points.data() + 3 * p, 3), entities);
       offsets[p + 1] = entities.size();
     }
 
@@ -320,7 +329,8 @@ std::vector<std::int32_t> geometry::compute_closest_entity(
       leaves = midpoint_tree.bbox(0);
       assert(is_leaf(leaves));
       initial_entity = leaves[0];
-      std::array<double, 3> diff = midpoint_tree.get_bbox(0)[0];
+      std::array<double, 6> bbox_coord = midpoint_tree.get_bbox(0);
+      auto diff = std::span(bbox_coord).subspan<0, 3>();
       for (std::size_t k = 0; k < 3; ++k)
         diff[k] -= points[3 * i + k];
       R2 = diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2];
@@ -330,8 +340,7 @@ std::vector<std::int32_t> geometry::compute_closest_entity(
       // As the midpoint tree only consist of points, the distance
       // queries are lightweight.
       const auto [m_index, m_distance2] = _compute_closest_entity(
-          midpoint_tree,
-          {points[3 * i + 0], points[3 * i + 1], points[3 * i + 2]},
+          midpoint_tree, std::span<const double, 3>(points.data() + 3 * i, 3),
           midpoint_tree.num_bboxes() - 1, mesh, initial_entity, R2);
 
       // Use a recursive search through the bounding box tree to
@@ -340,7 +349,7 @@ std::vector<std::int32_t> geometry::compute_closest_entity(
       // the distance from the midpoint to the point of interest as the
       // initial search radius.
       const auto [index, distance2] = _compute_closest_entity(
-          tree, {points[3 * i + 0], points[3 * i + 1], points[3 * i + 2]},
+          tree, std::span<const double, 3>(points.data() + 3 * i, 3),
           tree.num_bboxes() - 1, mesh, m_index, m_distance2);
 
       entities.push_back(index);
@@ -351,12 +360,13 @@ std::vector<std::int32_t> geometry::compute_closest_entity(
 }
 
 //-----------------------------------------------------------------------------
-double geometry::compute_squared_distance_bbox(
-    const std::array<std::array<double, 3>, 2>& b,
-    const std::array<double, 3>& x)
+double geometry::compute_squared_distance_bbox(std::span<const double, 6> b,
+                                               std::span<const double, 3> x)
 {
-  auto& b0 = b[0];
-  auto& b1 = b[1];
+  assert(b.size() == 6);
+  auto b0 = b.subspan<0, 3>();
+  auto b1 = b.subspan<3, 3>();
+
   return std::transform_reduce(x.begin(), x.end(), b0.begin(), 0.0,
                                std::plus<>{},
                                [](auto x, auto b)
diff --git a/cpp/dolfinx/geometry/utils.h b/cpp/dolfinx/geometry/utils.h
index 0ddaa16f5a2..ae982700089 100644
--- a/cpp/dolfinx/geometry/utils.h
+++ b/cpp/dolfinx/geometry/utils.h
@@ -35,9 +35,10 @@ create_midpoint_tree(const mesh::Mesh& mesh, int tdim,
 /// process)
 /// @param[in] tree0 First BoundingBoxTree
 /// @param[in] tree1 Second BoundingBoxTree
-/// @return List of pairs of intersecting box indices from each tree
-std::vector<std::array<int, 2>>
-compute_collisions(const BoundingBoxTree& tree0, const BoundingBoxTree& tree1);
+/// @return List of pairs of intersecting box indices from each tree, flattened
+/// as a vector of size num_intersections*2
+std::vector<int> compute_collisions(const BoundingBoxTree& tree0,
+                                    const BoundingBoxTree& tree1);
 
 /// Compute all collisions between bounding boxes and for a set of
 /// points
@@ -69,9 +70,8 @@ std::vector<std::int32_t> compute_closest_entity(
 /// @param[in] x A point
 /// @return The shortest distance between the bounding box `b` and the
 /// point `x`. Returns zero if `x` is inside box.
-double
-compute_squared_distance_bbox(const std::array<std::array<double, 3>, 2>& b,
-                              const std::array<double, 3>& x);
+double compute_squared_distance_bbox(std::span<const double, 6> b,
+                                     std::span<const double, 3> x);
 
 /// Compute the shortest vector from a mesh entity to a point
 /// @param[in] mesh The mesh
diff --git a/python/dolfinx/wrappers/geometry.cpp b/python/dolfinx/wrappers/geometry.cpp
index ca0b0e749e7..1cbe3f63d2b 100644
--- a/python/dolfinx/wrappers/geometry.cpp
+++ b/python/dolfinx/wrappers/geometry.cpp
@@ -6,6 +6,7 @@
 
 #include "array.h"
 #include "caster_mpi.h"
+#include <dolfinx/common/utils.h>
 #include <dolfinx/geometry/BoundingBoxTree.h>
 #include <dolfinx/geometry/gjk.h>
 #include <dolfinx/geometry/utils.h>
@@ -89,11 +90,18 @@ void geometry(py::module& m)
         return dolfinx::geometry::compute_collisions(tree, _p);
       },
       py::arg("tree"), py::arg("points"));
-  m.def("compute_collisions",
-        py::overload_cast<const dolfinx::geometry::BoundingBoxTree&,
-                          const dolfinx::geometry::BoundingBoxTree&>(
-            &dolfinx::geometry::compute_collisions),
-        py::arg("tree0"), py::arg("tree1"));
+  m.def(
+      "compute_collisions",
+      [](const dolfinx::geometry::BoundingBoxTree& treeA,
+         const dolfinx::geometry::BoundingBoxTree& treeB)
+      {
+        std::vector<std::int32_t> coll
+            = dolfinx::geometry::compute_collisions(treeA, treeB);
+
+        std::array<py::ssize_t, 2> shape = {py::ssize_t(coll.size() / 2), 2};
+        return as_pyarray(std::move(coll), shape);
+      },
+      py::arg("tree0"), py::arg("tree1"));
 
   m.def(
       "compute_distance_gjk",
@@ -231,7 +239,14 @@ void geometry(py::module& m)
       .def(
           "get_bbox",
           [](const dolfinx::geometry::BoundingBoxTree& self,
-             const std::size_t i) { return self.get_bbox(i); },
+             const std::size_t i)
+          {
+            std::array<double, 6> bbox = self.get_bbox(i);
+            std::array<std::size_t, 2> shape = {2, 3};
+            std::vector<double> bbox_out(6);
+            std::copy_n(bbox.begin(), 6, bbox_out.begin());
+            return dolfinx_wrappers::as_pyarray(std::move(bbox_out), shape);
+          },
           py::arg("i"))
       .def("__repr__", &dolfinx::geometry::BoundingBoxTree::str)
       .def(