diff --git a/.codecov.yml b/.codecov.yml
index fc495410..db247200 100644
--- a/.codecov.yml
+++ b/.codecov.yml
@@ -1 +1 @@
-comments: off
\ No newline at end of file
+comment: off
diff --git a/README.md b/README.md
index 840fd94b..77fd9f5a 100644
--- a/README.md
+++ b/README.md
@@ -37,11 +37,11 @@ int main() {
     auto h = make_histogram( axis::regular<>(4, 0.0, 2.0) );
 
     // push some values into the histogram
-    for (auto value : { 0.4, 1.1, 0.3, 1.7, 10. })
+    for (auto&& value : { 0.4, 1.1, 0.3, 1.7, 10. })
       h(value);
 
     // iterate over bins
-    for (auto x : indexed(h)) {
+    for (auto&& x : indexed(h)) {
       std::cout << boost::format("bin %i [ %.1f, %.1f ): %i\n")
         % x.index() % x.bin().lower() % x.bin().upper() % *x;
     }
diff --git a/doc/benchmarks.qbk b/doc/benchmarks.qbk
index 17c06030..f1541283 100644
--- a/doc/benchmarks.qbk
+++ b/doc/benchmarks.qbk
@@ -46,7 +46,7 @@ for (int i = 0; i < h.axis(0).size(); ++i) {
 }
 
 // same, with indexed range generator
-for (auto x : boost::histogram::indexed(h)) {
+for (auto&& x : boost::histogram::indexed(h)) {
   std::cout << x.index(0) << " " << x.index(1) << " " << *x << std::endl;
 }
 ```
diff --git a/doc/doxygen_postprocessing.py b/doc/doxygen_postprocessing.py
index edfbc2f7..03c25906 100644
--- a/doc/doxygen_postprocessing.py
+++ b/doc/doxygen_postprocessing.py
@@ -30,7 +30,14 @@ def is_detail(x):
 
     p = x.find("purpose")
     if p is not None:
-        return p.text.lower().strip() == "implementation detail"
+        return p.text.lower().lstrip().startswith("implementation detail")
+    return False
+
+
+def is_deprecated(x):
+    p = x.find("purpose")
+    if p is not None:
+        return p.text.lower().lstrip().startswith("deprecated")
     return False
 
 
@@ -60,6 +67,12 @@ for item in select(is_detail, "type"):
     item.clear()
     item.append(unspecified)
 
+# hide everything that's deprecated
+for item in select(is_deprecated, "typedef"):
+    parent = parent_map[item]
+    log("removing deprecated", item.tag, item.get("name"), "from", parent.tag, parent.get("name"))
+    parent.remove(item)
+
 # hide private member functions
 for item in select(lambda x: x.get("name") == "private member functions", "method-group"):
     parent = parent_map[item]
diff --git a/doc/guide.qbk b/doc/guide.qbk
index 4c3013f8..f0fcfb8f 100644
--- a/doc/guide.qbk
+++ b/doc/guide.qbk
@@ -33,7 +33,7 @@ The factory function named [headerref boost/histogram/make_histogram.hpp make_hi
 The library provides a number of useful axis types. The builtin axis types can be configured to fit many needs. If you still need something more exotic, no problem, it is easy to write your own axis types, see [link histogram.guide.expert Advanced usage]. In the following, we give some advice when to use which of the builtin axis types.
 
 [section Overview]
- 
+
 [variablelist
   [
     [
@@ -203,7 +203,7 @@ When the weights come from a stochastic process, it is useful to keep track of t
 [import ../examples/guide_fill_weighted_histogram.cpp]
 [guide_fill_weighted_histogram]
 
-When iterating over all cells, using [memberref boost::histogram::histogram::at histogram::at] can be inconvenient. The [funcref boost::histogram::indexed indexed] range generator is provided for this case, which is very convenient and /faster/ than naive for-loops.
+To iterate over all cells, the [funcref boost::histogram::indexed indexed] range generator is very convenient and also efficient. For almost all configurations, the range generator iterates /faster/ than a naive for-loop. Under- and overflow are skipped by default.
 
 [import ../examples/guide_indexed_access.cpp]
 [guide_indexed_access]
@@ -212,7 +212,7 @@ When iterating over all cells, using [memberref boost::histogram::histogram::at
 
 [section Using profiles]
 
-Histograms from this library can do more than counting, they can hold arbitrary accumulators which accept samples. A histogram is called a /profile/, if it computes the means of the samples in each cell. 
+Histograms from this library can do more than counting, they can hold arbitrary accumulators which accept samples. A histogram is called a /profile/, if it computes the means of the samples in each cell.
 
 Profiles can be generated with the factory function [funcref boost::histogram::make_profile make_profile].
 
diff --git a/examples/getting_started_listing_01.cpp b/examples/getting_started_listing_01.cpp
index eaf8adb3..54e9bf0c 100644
--- a/examples/getting_started_listing_01.cpp
+++ b/examples/getting_started_listing_01.cpp
@@ -70,7 +70,7 @@ int main() {
   */
 
   std::ostringstream os;
-  for (auto x : indexed(h, coverage::all)) {
+  for (auto&& x : indexed(h, coverage::all)) {
     os << boost::format("bin %2i [%4.1f, %4.1f): %i\n") % x.index() % x.bin().lower() %
               x.bin().upper() % *x;
   }
diff --git a/examples/getting_started_listing_02.cpp b/examples/getting_started_listing_02.cpp
index 80bee986..11418bb8 100644
--- a/examples/getting_started_listing_02.cpp
+++ b/examples/getting_started_listing_02.cpp
@@ -48,7 +48,7 @@ int main() {
   // get reference to category axis, performs a run-time checked static cast
   const auto& cat_axis = axis::get<cat>(h.axis(0));
   std::ostringstream os;
-  for (auto x : indexed(h)) {
+  for (auto&& x : indexed(h)) {
     os << boost::format("(%i, %i, %i) %4s [%3.1f, %3.1f) [%3.1f, %3.1f) %3.0f\n") %
               x.index(0) % x.index(1) % x.index(2) % cat_axis.bin(x.index(0)) %
               x.bin(1).lower() % x.bin(1).upper() % x.bin(2).lower() % x.bin(2).upper() %
diff --git a/examples/getting_started_listing_03.cpp b/examples/getting_started_listing_03.cpp
index 6f7977f7..3b5e2eda 100644
--- a/examples/getting_started_listing_03.cpp
+++ b/examples/getting_started_listing_03.cpp
@@ -35,7 +35,7 @@ int main() {
     Iterate over bins and print profile.
   */
   std::ostringstream os;
-  for (auto x : indexed(p)) {
+  for (auto&& x : indexed(p)) {
     os << boost::format("bin %i [%3.1f, %3.1f) count %i mean %g\n") % x.index() %
               x.bin().lower() % x.bin().upper() % x->count() % x->value();
   }
diff --git a/examples/guide_custom_accumulators.cpp b/examples/guide_custom_accumulators.cpp
index a2c41a8e..eea68dfe 100644
--- a/examples/guide_custom_accumulators.cpp
+++ b/examples/guide_custom_accumulators.cpp
@@ -27,7 +27,7 @@ int main() {
   h1(0.5, sample(4)); // sample 4 goes to second bin
 
   std::ostringstream os1;
-  for (auto x : indexed(h1)) {
+  for (auto&& x : indexed(h1)) {
     // Accumulators usually have methods to access their state. Use the arrow
     // operator to access them. Here, `count()` gives the number of samples,
     // `value()` the mean, and `variance()` the variance estimate of the mean.
@@ -57,7 +57,7 @@ int main() {
   h2(0.5, sample(4));   // sample 4 goes to second bin
 
   std::ostringstream os2;
-  for (auto x : indexed(h2)) {
+  for (auto&& x : indexed(h2)) {
     os2 << boost::format("%i value %.1f\n") % x.index() % x->value;
   }
   std::cout << os2.str() << std::flush;
diff --git a/examples/guide_fill_profile.cpp b/examples/guide_fill_profile.cpp
index cda9d11a..687582a3 100644
--- a/examples/guide_fill_profile.cpp
+++ b/examples/guide_fill_profile.cpp
@@ -33,7 +33,7 @@ int main() {
 
   // builtin accumulators have methods to access their state
   std::ostringstream os;
-  for (auto x : indexed(h)) {
+  for (auto&& x : indexed(h)) {
     // use `.` to access methods of accessor, like `index()`
     // use `->` to access methods of accumulator
     const auto i = x.index();
diff --git a/examples/guide_histogram_projection.cpp b/examples/guide_histogram_projection.cpp
index 4f598583..54f45212 100644
--- a/examples/guide_histogram_projection.cpp
+++ b/examples/guide_histogram_projection.cpp
@@ -31,7 +31,7 @@ int main() {
   assert(algorithm::sum(h) == 3 && algorithm::sum(hr0) == 3 && algorithm::sum(hr1) == 3);
 
   std::ostringstream os1;
-  for (auto x : indexed(h))
+  for (auto&& x : indexed(h))
     os1 << "(" << x.index(0) << ", " << x.index(1) << "): " << *x << "\n";
   std::cout << os1.str() << std::flush;
   assert(os1.str() == "(0, 0): 1\n"
@@ -42,14 +42,14 @@ int main() {
                       "(2, 1): 1\n");
 
   std::ostringstream os2;
-  for (auto x : indexed(hr0)) os2 << "(" << x.index(0) << ", -): " << *x << "\n";
+  for (auto&& x : indexed(hr0)) os2 << "(" << x.index(0) << ", -): " << *x << "\n";
   std::cout << os2.str() << std::flush;
   assert(os2.str() == "(0, -): 1\n"
                       "(1, -): 1\n"
                       "(2, -): 1\n");
 
   std::ostringstream os3;
-  for (auto x : indexed(hr1)) os3 << "(- ," << x.index(0) << "): " << *x << "\n";
+  for (auto&& x : indexed(hr1)) os3 << "(- ," << x.index(0) << "): " << *x << "\n";
   std::cout << os3.str() << std::flush;
   assert(os3.str() == "(- ,0): 2\n"
                       "(- ,1): 1\n");
diff --git a/examples/guide_indexed_access.cpp b/examples/guide_indexed_access.cpp
index 65297a7a..5850c1d3 100644
--- a/examples/guide_indexed_access.cpp
+++ b/examples/guide_indexed_access.cpp
@@ -27,7 +27,7 @@ int main() {
   // use the `indexed` range adaptor to iterate over all bins;
   // it is not only more convenient but also faster than a hand-crafted loop!
   std::ostringstream os;
-  for (auto x : indexed(h)) {
+  for (auto&& x : indexed(h)) {
     // x is a special accessor object
     const auto i = x.index(0); // current index along first axis
     const auto j = x.index(1); // current index along second axis
@@ -48,7 +48,7 @@ int main() {
   // `indexed` skips underflow and overflow bins by default, but can be called
   // with the second argument `coverage::all` to walk over all bins
   std::ostringstream os2;
-  for (auto x : indexed(h, coverage::all)) {
+  for (auto&& x : indexed(h, coverage::all)) {
     os2 << boost::format("%2i %2i: %i\n") % x.index(0) % x.index(1) % *x;
   }
 
diff --git a/examples/guide_stdlib_algorithms.cpp b/examples/guide_stdlib_algorithms.cpp
index 73691bb0..c010a2f5 100644
--- a/examples/guide_stdlib_algorithms.cpp
+++ b/examples/guide_stdlib_algorithms.cpp
@@ -19,35 +19,33 @@ int main() {
   // make histogram that represents a probability density function (PDF)
   auto h1 = make_histogram(axis::regular<>(4, 1.0, 3.0));
 
-  // use std::fill to set all counters to 0.25, including *flow cells
-  std::fill(h1.begin(), h1.end(), 0.25);
-  // reset *flow cells to zero
-  h1.at(-1) = h1.at(4) = 0;
+  // make indexed range to skip underflow and overflow cells
+  auto ind = indexed(h1);
+
+  // use std::fill to set all counters to 0.25 (except under- and overflow counters)
+  std::fill(ind.begin(), ind.end(), 0.25);
 
   // compute the cumulative density function (CDF), overriding cell values
-  std::partial_sum(h1.begin(), h1.end(), h1.begin());
+  std::partial_sum(ind.begin(), ind.end(), ind.begin());
 
-  assert(h1.at(-1) == 0.0);
+  assert(h1.at(-1) == 0.00);
   assert(h1.at(0) == 0.25);
   assert(h1.at(1) == 0.50);
   assert(h1.at(2) == 0.75);
   assert(h1.at(3) == 1.00);
-  assert(h1.at(4) == 1.00);
+  assert(h1.at(4) == 0.00);
 
   // use any_of to check if any cell values are smaller than 0.1,
-  // and use indexed() to skip underflow and overflow cells
-  auto h1_ind = indexed(h1);
-  const auto any_small =
-      std::any_of(h1_ind.begin(), h1_ind.end(), [](const auto& x) { return *x < 0.1; });
-  assert(any_small == false); // underflow and overflow are zero, but skipped
-
-  // find maximum element
-  const auto max_it = std::max_element(h1.begin(), h1.end());
-  assert(max_it == h1.end() - 2);
+  auto b = std::any_of(ind.begin(), ind.end(), [](const auto& x) { return *x < 0.1; });
+  assert(b == false); // under- and overflow cells are zero, but skipped
 
   // find minimum element
-  const auto min_it = std::min_element(h1.begin(), h1.end());
-  assert(min_it == h1.begin());
+  auto min_it = std::min_element(ind.begin(), ind.end());
+  assert(*min_it == 0.25); // under- and overflow cells are skipped
+
+  // find maximum element
+  auto max_it = std::max_element(ind.begin(), ind.end());
+  assert(*max_it == 1.0);
 
   // make second PDF
   auto h2 = make_histogram(axis::regular<>(4, 1.0, 4.0));
@@ -57,12 +55,13 @@ int main() {
   h2.at(3) = 0.4;
 
   // computing cosine similiarity: cos(theta) = A dot B / sqrt((A dot A) * (B dot B))
-  const auto aa = std::inner_product(h1.begin(), h1.end(), h1.begin(), 0.0);
-  const auto bb = std::inner_product(h2.begin(), h2.end(), h2.begin(), 0.0);
-  const auto ab = std::inner_product(h1.begin(), h1.end(), h2.begin(), 0.0);
+  auto ind2 = indexed(h2);
+  const auto aa = std::inner_product(ind.begin(), ind.end(), ind.begin(), 0.0);
+  const auto bb = std::inner_product(ind2.begin(), ind2.end(), ind2.begin(), 0.0);
+  const auto ab = std::inner_product(ind.begin(), ind.end(), ind2.begin(), 0.0);
   const auto cos_sim = ab / std::sqrt(aa * bb);
 
-  assert(std::abs(cos_sim - 0.78) < 1e-2);
+  assert(std::abs(cos_sim - 0.967) < 1e-2);
 }
 
 //]
diff --git a/include/boost/histogram/algorithm/project.hpp b/include/boost/histogram/algorithm/project.hpp
index d0ae25b5..6c0495c8 100644
--- a/include/boost/histogram/algorithm/project.hpp
+++ b/include/boost/histogram/algorithm/project.hpp
@@ -53,7 +53,7 @@ auto project(const histogram<A, S>& h, std::integral_constant<unsigned, N>, Ns..
   using A2 = decltype(axes);
   auto result = histogram<A2, S>(std::move(axes), detail::make_default(old_storage));
   auto idx = detail::make_stack_buffer<int>(unsafe_access::axes(result));
-  for (auto x : indexed(h, coverage::all)) {
+  for (auto&& x : indexed(h, coverage::all)) {
     auto i = idx.begin();
     mp11::mp_for_each<LN>([&i, &x](auto J) { *i++ = x.index(J); });
     result.at(idx) += *x;
@@ -86,7 +86,7 @@ auto project(const histogram<A, S>& h, const Iterable& c) {
   auto result =
       histogram<decltype(axes), S>(std::move(axes), detail::make_default(old_storage));
   auto idx = detail::make_stack_buffer<int>(unsafe_access::axes(result));
-  for (auto x : indexed(h, coverage::all)) {
+  for (auto&& x : indexed(h, coverage::all)) {
     auto i = idx.begin();
     for (auto d : c) *i++ = x.index(d);
     result.at(idx) += *x;
diff --git a/include/boost/histogram/algorithm/reduce.hpp b/include/boost/histogram/algorithm/reduce.hpp
index 8236d5cf..ba71cb72 100644
--- a/include/boost/histogram/algorithm/reduce.hpp
+++ b/include/boost/histogram/algorithm/reduce.hpp
@@ -269,7 +269,7 @@ decltype(auto) reduce(const Histogram& hist, const Iterable& options) {
   auto result = Histogram(std::move(axes), std::move(storage));
 
   auto idx = detail::make_stack_buffer<int>(unsafe_access::axes(result));
-  for (auto x : indexed(hist, coverage::all)) {
+  for (auto&& x : indexed(hist, coverage::all)) {
     auto i = idx.begin();
     auto o = opts.begin();
     for (auto j : x.indices()) {
diff --git a/include/boost/histogram/algorithm/sum.hpp b/include/boost/histogram/algorithm/sum.hpp
index 63316ade..5ff437dd 100644
--- a/include/boost/histogram/algorithm/sum.hpp
+++ b/include/boost/histogram/algorithm/sum.hpp
@@ -30,7 +30,7 @@ auto sum(const histogram<A, S>& h) {
   using T = typename histogram<A, S>::value_type;
   using Sum = mp11::mp_if<std::is_arithmetic<T>, accumulators::sum<double>, T>;
   Sum sum;
-  for (auto x : h) sum += x;
+  for (auto&& x : h) sum += x;
   using R = mp11::mp_if<std::is_arithmetic<T>, double, T>;
   return static_cast<R>(sum);
 }
diff --git a/include/boost/histogram/detail/axes.hpp b/include/boost/histogram/detail/axes.hpp
index 7e484665..1372936d 100644
--- a/include/boost/histogram/detail/axes.hpp
+++ b/include/boost/histogram/detail/axes.hpp
@@ -223,19 +223,25 @@ using buffer_size = mp11::mp_eval_or<
 
 template <class T, std::size_t N>
 class sub_array : public std::array<T, N> {
+  using base_type = std::array<T, N>;
+
 public:
-  explicit sub_array(std::size_t s) : size_(s) {
+  explicit sub_array(std::size_t s) noexcept(
+      std::is_nothrow_default_constructible<T>::value)
+      : size_(s) {
     BOOST_ASSERT_MSG(size_ <= N, "requested size exceeds size of static buffer");
   }
 
-  sub_array(std::size_t s, const T& value) : size_(s) {
+  sub_array(std::size_t s,
+            const T& value) noexcept(std::is_nothrow_copy_constructible<T>::value)
+      : size_(s) {
     BOOST_ASSERT_MSG(size_ <= N, "requested size exceeds size of static buffer");
     std::array<T, N>::fill(value);
   }
 
   // need to override both versions of std::array
-  auto end() noexcept { return std::array<T, N>::begin() + size_; }
-  auto end() const noexcept { return std::array<T, N>::begin() + size_; }
+  auto end() noexcept { return base_type::begin() + size_; }
+  auto end() const noexcept { return base_type::begin() + size_; }
 
   auto size() const noexcept { return size_; }
 
diff --git a/include/boost/histogram/indexed.hpp b/include/boost/histogram/indexed.hpp
index c5d38fef..8cba2ca6 100644
--- a/include/boost/histogram/indexed.hpp
+++ b/include/boost/histogram/indexed.hpp
@@ -7,6 +7,7 @@
 #ifndef BOOST_HISTOGRAM_INDEXED_HPP
 #define BOOST_HISTOGRAM_INDEXED_HPP
 
+#include <array>
 #include <boost/histogram/axis/traits.hpp>
 #include <boost/histogram/detail/attribute.hpp>
 #include <boost/histogram/detail/axes.hpp>
@@ -38,54 +39,47 @@ class BOOST_HISTOGRAM_NODISCARD indexed_range {
 private:
   using histogram_type = Histogram;
   static constexpr std::size_t buffer_size =
-      detail::buffer_size<typename std::decay_t<histogram_type>::axes_type>::value;
+      detail::buffer_size<typename std::remove_const_t<histogram_type>::axes_type>::value;
 
 public:
-  using value_iterator = decltype(std::declval<Histogram>().begin());
+  using value_iterator = decltype(std::declval<histogram_type>().begin());
   using value_reference = typename value_iterator::reference;
-  class range_iterator;
+  using value_type = typename value_iterator::value_type;
 
-private:
-  struct state_type {
-    struct index_data {
-      axis::index_type idx, begin, end, extent;
-    };
-
-    state_type(histogram_type& h) : hist_(h) {}
-
-    histogram_type& hist_;
-    index_data indices_[buffer_size];
-  };
-
-public:
-  class detached_accessor; // forward declaration
+  class iterator;
+  using range_iterator = iterator; ///< deprecated
 
   /** Pointer-like class to access value and index of current cell.
 
-    Its methods allow one to query the current indices and bins. Furthermore, it acts
-    like a pointer to the cell value. The accessor is coupled to the current
-    range_iterator. Moving the range_iterator forward invalidates the accessor. Use the
-    detached_accessor class if you must store accessors for later use, but be aware
-    that a detached_accessor has a state many times larger than a pointer.
+    Its methods provide access to the current indices and bins and it acts like a pointer
+    to the cell value. To interoperate with the algorithms of the standard library, the
+    accessor is implicitly convertible to a cell value. Assignments and comparisons
+    are passed through to the cell. The accessor is coupled to its parent
+    iterator. Moving the parent iterator forward also updates the linked
+    accessor. Accessors are not copyable. They cannot be stored in containers, but
+    range_iterators can be stored.
   */
   class accessor {
   public:
     /// Array-like view into the current multi-dimensional index.
     class index_view {
-      using index_pointer = const typename state_type::index_data*;
+      using index_pointer = const typename iterator::index_data*;
 
     public:
-      using reference = const axis::index_type&;
+      using const_reference = const axis::index_type&;
+      using reference = const_reference; ///< deprecated
 
       /// implementation detail
       class const_iterator
-          : public detail::iterator_adaptor<const_iterator, index_pointer, reference> {
+          : public detail::iterator_adaptor<const_iterator, index_pointer,
+                                            const_reference> {
       public:
-        reference operator*() const noexcept { return const_iterator::base()->idx; }
+        const_reference operator*() const noexcept { return const_iterator::base()->idx; }
 
       private:
         explicit const_iterator(index_pointer i) noexcept
             : const_iterator::iterator_adaptor_(i) {}
+
         friend class index_view;
       };
 
@@ -94,8 +88,8 @@ public:
       std::size_t size() const noexcept {
         return static_cast<std::size_t>(end_ - begin_);
       }
-      reference operator[](unsigned d) const noexcept { return begin_[d].idx; }
-      reference at(unsigned d) const { return begin_[d].idx; }
+      const_reference operator[](unsigned d) const noexcept { return begin_[d].idx; }
+      const_reference at(unsigned d) const { return begin_[d].idx; }
 
     private:
       /// implementation detail
@@ -105,34 +99,54 @@ public:
       friend class accessor;
     };
 
+    // assignment is pass-through
+    accessor& operator=(const accessor& o) {
+      get() = o.get();
+      return *this;
+    }
+
+    // assignment is pass-through
+    template <class T>
+    accessor& operator=(const T& x) {
+      get() = x;
+      return *this;
+    }
+
     /// Returns the cell reference.
-    value_reference get() const noexcept { return *iter_; }
+    value_reference get() const noexcept { return *(iter_.iter_); }
     /// @copydoc get()
     value_reference operator*() const noexcept { return get(); }
     /// Access fields and methods of the cell object.
-    value_iterator operator->() const noexcept { return iter_; }
+    value_iterator operator->() const noexcept { return iter_.iter_; }
 
     /// Access current index.
     /// @param d axis dimension.
     axis::index_type index(unsigned d = 0) const noexcept {
-      return state_.indices_[d].idx;
+      return iter_.indices_[d].idx;
     }
 
     /// Access indices as an iterable range.
     index_view indices() const noexcept {
-      return {state_.indices_, state_.indices_ + state_.hist_.rank()};
+      BOOST_ASSERT(iter_.indices_.hist_);
+      return {iter_.indices_.data(),
+              iter_.indices_.data() + iter_.indices_.hist_->rank()};
     }
 
     /// Access current bin.
     /// @tparam N axis dimension.
     template <unsigned N = 0>
     decltype(auto) bin(std::integral_constant<unsigned, N> = {}) const {
-      return state_.hist_.axis(std::integral_constant<unsigned, N>()).bin(index(N));
+      BOOST_ASSERT(iter_.indices_.hist_);
+      return iter_.indices_.hist_->axis(std::integral_constant<unsigned, N>())
+          .bin(index(N));
     }
 
     /// Access current bin.
     /// @param d axis dimension.
-    decltype(auto) bin(unsigned d) const { return state_.hist_.axis(d).bin(index(d)); }
+    decltype(auto) bin(unsigned d) const {
+      BOOST_ASSERT(iter_.indices_.hist_);
+      return iter_.indices_.hist_->axis(d).bin(index(d));
+    }
 
     /** Computes density in current cell.
 
@@ -140,155 +154,237 @@ public:
       without bin widths, like axis::category, are treated as having unit bin with.
     */
     double density() const {
+      BOOST_ASSERT(iter_.indices_.hist_);
       double x = 1;
       unsigned d = 0;
-      state_.hist_.for_each_axis([&](const auto& a) {
+      iter_.indices_.hist_->for_each_axis([&](const auto& a) {
         const auto w = axis::traits::width_as<double>(a, this->index(d++));
         x *= w ? w : 1;
       });
       return get() / x;
     }
 
-  protected:
-    accessor(state_type& s, value_iterator i) noexcept : state_(s), iter_(i) {}
+    // forward all comparison operators to the value
 
-    state_type& state_;
-    value_iterator iter_;
-    friend class range_iterator;
-    friend class detached_accessor;
-  };
-
-  /// Accessor that owns a copy of the iterator state.
-  class detached_accessor : public accessor {
-  public:
-    detached_accessor(const accessor& x) : accessor(state_, x.iter_), state_(x.state_) {}
-    detached_accessor(const detached_accessor& x)
-        : detached_accessor(static_cast<const accessor&>(x)) {}
-    detached_accessor& operator=(const detached_accessor& x) {
-      state_ = x.state_;
-      accessor::iter_ = x.iter_;
-      return *this;
+    template <class U>
+    bool operator<(const U& o) const noexcept {
+      return get() < o;
     }
 
+    template <class U>
+    bool operator>(const U& o) const noexcept {
+      return get() > o;
+    }
+
+    template <class U>
+    bool operator==(const U& o) const noexcept {
+      return get() == o;
+    }
+
+    template <class U>
+    bool operator!=(const U& o) const noexcept {
+      return get() != o;
+    }
+
+    template <class U>
+    bool operator<=(const U& o) const noexcept {
+      return get() <= o;
+    }
+
+    template <class U>
+    bool operator>=(const U& o) const noexcept {
+      return get() >= o;
+    }
+
+    template <class U>
+    friend bool operator<(const U& x, const accessor& y) noexcept {
+      return y.operator>(x);
+    }
+
+    template <class U>
+    friend bool operator>(const U& x, const accessor& y) noexcept {
+      return y.operator<(x);
+    }
+
+    template <class U>
+    friend bool operator<=(const U& x, const accessor& y) noexcept {
+      return y.operator>=(x);
+    }
+
+    template <class U>
+    friend bool operator>=(const U& x, const accessor& y) noexcept {
+      return y.operator<=(x);
+    }
+
+    template <class U>
+    friend bool operator==(const U& x, const accessor& y) noexcept {
+      return y.operator==(x);
+    }
+
+    template <class U>
+    friend bool operator!=(const U& x, const accessor& y) noexcept {
+      return y.operator!=(x);
+    }
+
+    friend bool operator<(const accessor& x, const accessor& y) noexcept {
+      return x.operator<(y.get());
+    }
+
+    friend bool operator>(const accessor& x, const accessor& y) noexcept {
+      return x.operator>(y.get());
+    }
+
+    friend bool operator==(const accessor& x, const accessor& y) noexcept {
+      return x.operator==(y.get());
+    }
+
+    friend bool operator!=(const accessor& x, const accessor& y) noexcept {
+      return x.operator!=(y.get());
+    }
+
+    friend bool operator<=(const accessor& x, const accessor& y) noexcept {
+      return x.operator<=(y.get());
+    }
+
+    friend bool operator>=(const accessor& x, const accessor& y) noexcept {
+      return x.operator>=(y.get());
+    }
+
+    operator value_type() const noexcept { return get(); }
+
   private:
-    state_type state_;
+    accessor(iterator& i) noexcept : iter_(i) {}
+
+    accessor(const accessor&) = default; // only callable by indexed_range::iterator
+
+    iterator& iter_;
+
+    friend class iterator;
   };
 
   /// implementation detail
-  class range_iterator {
+  class iterator {
   public:
-    using value_type = accessor;
-    using reference = accessor&;
-    using pointer = accessor*;
-    using difference_type = void;
-    using iterator_category = std::input_iterator_tag;
+    using value_type = typename indexed_range::value_type;
+    using reference = accessor;
 
   private:
-    struct value_proxy {
-      detached_accessor operator*() { return ref; }
-      detached_accessor ref;
+    struct pointer_proxy {
+      reference* operator->() noexcept { return std::addressof(ref_); }
+      reference ref_;
     };
 
   public:
-    reference operator*() noexcept { return value_; }
-    pointer operator->() noexcept { return &value_; }
+    using pointer = pointer_proxy;
+    using difference_type = std::ptrdiff_t;
+    using iterator_category = std::forward_iterator_tag;
 
-    range_iterator& operator++() {
+    reference operator*() noexcept { return *this; }
+    pointer operator->() noexcept { return pointer_proxy{operator*()}; }
+
+    iterator& operator++() {
+      BOOST_ASSERT(indices_.hist_);
       std::size_t stride = 1;
-      auto c = value_.state_.indices_;
+      auto c = indices_.begin();
       ++c->idx;
-      ++value_.iter_;
-      while (c->idx == c->end &&
-             (c != (value_.state_.indices_ + value_.state_.hist_.rank() - 1))) {
+      ++iter_;
+      while (c->idx == c->end && (c != (indices_.end() - 1))) {
         c->idx = c->begin;
-        value_.iter_ -= (c->end - c->begin) * stride;
+        iter_ -= (c->end - c->begin) * stride;
         stride *= c->extent;
         ++c;
         ++c->idx;
-        value_.iter_ += stride;
+        iter_ += stride;
       }
       return *this;
     }
 
-    value_proxy operator++(int) {
-      value_proxy x{value_};
+    iterator operator++(int) {
+      auto prev = *this;
       operator++();
-      return x;
+      return prev;
     }
 
-    bool operator==(const range_iterator& x) const noexcept {
-      return value_.iter_ == x.value_.iter_;
-    }
-    bool operator!=(const range_iterator& x) const noexcept { return !operator==(x); }
+    bool operator==(const iterator& x) const noexcept { return iter_ == x.iter_; }
+    bool operator!=(const iterator& x) const noexcept { return !operator==(x); }
 
   private:
-    range_iterator(state_type& s, value_iterator i) noexcept : value_(s, i) {}
+    iterator(histogram_type* h) : iter_(h->begin()), indices_(h) {}
+
+    value_iterator iter_;
+
+    struct index_data {
+      axis::index_type idx, begin, end, extent;
+    };
+
+    struct indices_t : std::array<index_data, buffer_size> {
+      using base_type = std::array<index_data, buffer_size>;
+
+      indices_t(histogram_type* h) noexcept : hist_{h} {}
+
+      constexpr auto end() noexcept { return base_type::begin() + hist_->rank(); }
+      constexpr auto end() const noexcept { return base_type::begin() + hist_->rank(); }
+      histogram_type* hist_;
+    } indices_;
 
-    accessor value_;
     friend class indexed_range;
   };
 
-  indexed_range(Histogram& hist, coverage cov)
-      : state_(hist), begin_(hist.begin()), end_(begin_) {
+  indexed_range(histogram_type& hist, coverage cov) : begin_(&hist), end_(&hist) {
     std::size_t stride = 1;
-    auto ca = state_.indices_;
-    const auto clast = ca + state_.hist_.rank() - 1;
-    state_.hist_.for_each_axis([ca, clast, cov, &stride, this](const auto& a) mutable {
-      using opt = axis::traits::static_options<decltype(a)>;
-      constexpr int under = opt::test(axis::option::underflow);
-      constexpr int over = opt::test(axis::option::overflow);
-      const auto size = a.size();
+    auto ca = begin_.indices_.begin();
+    const auto clast = ca + begin_.indices_.hist_->rank() - 1;
+    begin_.indices_.hist_->for_each_axis(
+        [ca, clast, cov, &stride, this](const auto& a) mutable {
+          using opt = axis::traits::static_options<decltype(a)>;
+          constexpr int under = opt::test(axis::option::underflow);
+          constexpr int over = opt::test(axis::option::overflow);
+          const auto size = a.size();
 
-      ca->extent = size + under + over;
-      // -1 if underflow and cover all, else 0
-      ca->begin = cov == coverage::all ? -under : 0;
-      // size + 1 if overflow and cover all, else size
-      ca->end = cov == coverage::all ? size + over : size;
-      ca->idx = ca->begin;
+          ca->extent = size + under + over;
+          // -1 if underflow and cover all, else 0
+          ca->begin = cov == coverage::all ? -under : 0;
+          // size + 1 if overflow and cover all, else size
+          ca->end = cov == coverage::all ? size + over : size;
+          ca->idx = ca->begin;
 
-      begin_ += (ca->begin + under) * stride;
-      end_ += ((ca < clast ? ca->begin : ca->end) + under) * stride;
+          begin_.iter_ += (ca->begin + under) * stride;
+          end_.iter_ += ((ca < clast ? ca->begin : ca->end) + under) * stride;
 
-      stride *= ca->extent;
-      ++ca;
-    });
+          stride *= ca->extent;
+          ++ca;
+        });
   }
 
-  range_iterator begin() noexcept {
-    auto begin = begin_;
-    begin_ = end_;
-    return {state_, begin};
-  }
-  range_iterator end() noexcept { return {state_, end_}; }
+  iterator begin() noexcept { return begin_; }
+  iterator end() noexcept { return end_; }
 
 private:
-  state_type state_;
-  value_iterator begin_, end_;
+  iterator begin_, end_;
 };
 
-/** Generates an indexed iterator range over the histogram cells.
+/** Generates an indexed range of <a
+  href="https://en.cppreference.com/w/cpp/named_req/ForwardIterator">forward iterators</a>
+  over the histogram cells.
 
   Use this in a range-based for loop:
 
   ```
-  for (auto x : indexed(hist)) { ... }
+  for (auto&& x : indexed(hist)) { ... }
   ```
 
-  This highly optimized loop is at least comparable in speed to a hand-written loop over
-  the histogram cells and often much faster, depending on the histogram configuration. The
-  iterators dereference to an indexed_range::accessor, which has methods to query the
-  current indices and bins and acts like a pointer to the cell value. Accessors, like
-  pointers, are cheap to copy but get invalidated when the range iterator is incremented.
-  Likewise, any copies of a range iterator become invalid if one of them is incremented.
-  A indexed_range::detached_accessor can be stored for later use, but manually copying the
-  data of interest from the accessor is usually more efficient.
+  This generates an optimized loop which is nearly always faster than a hand-written loop
+  over the histogram cells. The iterators dereference to an indexed_range::accessor, which
+  has methods to query the current indices and bins and acts like a pointer to the cell
+  value. The returned iterators are forward iterators. They can be stored in a container,
+  but may not be used after the life-time of the histogram ends.
 
   @returns indexed_range
 
   @param hist Reference to the histogram.
   @param cov  Iterate over all or only inner bins (optional, default: inner).
  */
-template <typename Histogram>
+template <class Histogram>
 auto indexed(Histogram&& hist, coverage cov = coverage::inner) {
   return indexed_range<std::remove_reference_t<Histogram>>{std::forward<Histogram>(hist),
                                                            cov};
diff --git a/include/boost/histogram/unlimited_storage.hpp b/include/boost/histogram/unlimited_storage.hpp
index 78ebb745..d851c1e1 100644
--- a/include/boost/histogram/unlimited_storage.hpp
+++ b/include/boost/histogram/unlimited_storage.hpp
@@ -106,16 +106,16 @@ void buffer_destroy(Allocator& a, typename std::allocator_traits<Allocator>::poi
 /**
   Memory-efficient storage for integral counters which cannot overflow.
 
-  This storage provides a no-overflow-guarantee if it is filled with integral weights
-  only. This storage implementation keeps a contiguous array of elemental counters, one
-  for each cell. If an operation is requested, which would overflow a counter, the whole
-  array is replaced with another of a wider integral type, then the operation is executed.
-  The storage uses integers of 8, 16, 32, 64 bits, and then switches to a multiprecision
+  This storage provides a no-overflow-guarantee if the counters are incremented with
+  integer weights. It maintains a contiguous array of elemental counters, one for each
+  cell. If an operation is requested which would overflow a counter, the array is
+  replaced with another of a wider integral type, then the operation is executed. The
+  storage uses integers of 8, 16, 32, 64 bits, and then switches to a multiprecision
   integral type, similar to those in
   [Boost.Multiprecision](https://www.boost.org/doc/libs/develop/libs/multiprecision/doc/html/index.html).
 
-  A scaling operation or adding a floating point number turns the elements into doubles,
-  which voids the no-overflow-guarantee.
+  A scaling operation or adding a floating point number triggers a conversion of the
+  elemental counters into doubles, which voids the no-overflow-guarantee.
 */
 template <class Allocator>
 class unlimited_storage {
@@ -525,7 +525,7 @@ public:
   const_iterator begin() const noexcept { return {&buffer_, 0}; }
   const_iterator end() const noexcept { return {&buffer_, size()}; }
 
-  /// @private used by unit tests, not part of generic storage interface
+  /// implementation detail; used by unit tests, not part of generic storage interface
   template <class T>
   unlimited_storage(std::size_t s, const T* p, const allocator_type& a = {})
       : buffer_(std::move(a)) {
diff --git a/test/indexed_test.cpp b/test/indexed_test.cpp
index b9d2557a..34d86439 100644
--- a/test/indexed_test.cpp
+++ b/test/indexed_test.cpp
@@ -4,6 +4,7 @@
 // (See accompanying file LICENSE_1_0.txt
 // or copy at http://www.boost.org/LICENSE_1_0.txt)
 
+#include <algorithm>
 #include <boost/core/lightweight_test.hpp>
 #include <boost/histogram/axis/integer.hpp>
 #include <boost/histogram/axis/ostream.hpp>
@@ -15,6 +16,7 @@
 #include <boost/mp11/algorithm.hpp>
 #include <boost/mp11/list.hpp>
 #include <iterator>
+#include <ostream>
 #include <type_traits>
 #include <vector>
 #include "throw_exception.hpp"
@@ -35,8 +37,6 @@ void run_1d_tests(mp_list<IsDynamic, Coverage>) {
 
   auto ind = indexed(h, Coverage());
   auto it = ind.begin();
-  // calling begin second time yields end iterator but does not confuse the original range
-  BOOST_TEST(ind.begin() == ind.end());
   BOOST_TEST_EQ(it->indices().size(), 1);
   BOOST_TEST_EQ(it->indices()[0], Coverage() == coverage::all ? -1 : 0);
 
@@ -55,10 +55,10 @@ void run_1d_tests(mp_list<IsDynamic, Coverage>) {
   BOOST_TEST_EQ(it->bin(0), h.axis().bin(1));
   ++it;
   // check post-increment
-  auto prev = *it++;
-  BOOST_TEST_EQ(prev.index(0), 2);
-  BOOST_TEST_EQ(*prev, 4);
-  BOOST_TEST_EQ(prev.bin(0), h.axis().bin(2));
+  auto prev = it++;
+  BOOST_TEST_EQ(prev->index(0), 2);
+  BOOST_TEST_EQ(**prev, 4);
+  BOOST_TEST_EQ(prev->bin(0), h.axis().bin(2));
   if (Coverage() == coverage::all) {
     BOOST_TEST_EQ(it->index(0), 3);
     BOOST_TEST_EQ(**it, 5);
@@ -115,11 +115,47 @@ void run_density_tests(mp_list<IsDynamic, Coverage>) {
   // fill uniformly
   for (auto&& x : h) x = 1;
 
-  for (auto x : indexed(h, Coverage())) {
+  for (auto&& x : indexed(h, Coverage())) {
     BOOST_TEST_EQ(x.density(), *x / (x.bin(0).width() * x.bin(2).width()));
   }
 }
 
+template <class IsDynamic, class Coverage>
+void run_stdlib_tests(mp_list<IsDynamic, Coverage>) {
+  auto ax = axis::regular<>(3, 0, 1);
+  auto ay = axis::integer<>(0, 2);
+  auto h = make_s(IsDynamic(), std::vector<int>(), ax, ay);
+
+  struct generator {
+    int i = 0;
+    int operator()() { return ++i; }
+  };
+
+  auto ind = indexed(h, Coverage());
+  std::generate(ind.begin(), ind.end(), generator{});
+
+  {
+    int i = 0;
+    for (auto&& x : ind) BOOST_TEST_EQ(*x, ++i);
+  }
+
+  {
+    auto it = std::min_element(ind.begin(), ind.end());
+    BOOST_TEST(it == ind.begin());
+    BOOST_TEST(it != ind.end());
+  }
+
+  {
+    auto it = std::max_element(ind.begin(), ind.end());
+    // get last before end()
+    auto it2 = ind.begin();
+    auto it3 = it2;
+    while (it2 != ind.end()) it3 = it2++;
+    BOOST_TEST(it == it3);
+    BOOST_TEST(it != ind.begin());
+  }
+}
+
 int main() {
   mp_for_each<mp_product<mp_list, mp_list<mp_false, mp_true>,
                          mp_list<std::integral_constant<coverage, coverage::inner>,
@@ -128,6 +164,7 @@ int main() {
         run_1d_tests(x);
         run_3d_tests(x);
         run_density_tests(x);
+        run_stdlib_tests(x);
       });
   return boost::report_errors();
 }