diff --git a/.clang-tidy b/.clang-tidy
index e467891b616d..447cf868cead 100644
--- a/.clang-tidy
+++ b/.clang-tidy
@@ -96,6 +96,9 @@ CheckOptions:
     value: true
   - key: performance-move-const-arg.CheckTriviallyCopyableMove
     value: false
+  # The fix for this is not supported by GCC 9.
+  - key: modernize-loop-convert.UseCxx20ReverseRanges
+    value: false
 WarningsAsErrors: '*'
 # It is unclear if the header filter actually works or how to use it so
 # just include all headers
diff --git a/src/Evolution/Executables/GeneralizedHarmonic/EvolveGhBinaryBlackHole.hpp b/src/Evolution/Executables/GeneralizedHarmonic/EvolveGhBinaryBlackHole.hpp
index e23a1ff8ac46..655677896a1b 100644
--- a/src/Evolution/Executables/GeneralizedHarmonic/EvolveGhBinaryBlackHole.hpp
+++ b/src/Evolution/Executables/GeneralizedHarmonic/EvolveGhBinaryBlackHole.hpp
@@ -488,7 +488,9 @@ struct EvolutionMetavars {
         tmpl::pair<
             gh::gauges::GaugeCondition,
             tmpl::list<gh::gauges::DampedHarmonic, gh::gauges::Harmonic>>,
-        tmpl::pair<LtsTimeStepper, TimeSteppers::lts_time_steppers>,
+        // Restrict to monotonic time steppers in LTS to avoid control
+        // systems deadlocking.
+        tmpl::pair<LtsTimeStepper, TimeSteppers::monotonic_lts_time_steppers>,
         tmpl::pair<PhaseChange, PhaseControl::factory_creatable_classes>,
         tmpl::pair<StepChooser<StepChooserUse::LtsStep>,
                    StepChoosers::standard_step_choosers<system>>,
diff --git a/src/Evolution/Executables/GeneralizedHarmonic/EvolveGhSingleBlackHole.hpp b/src/Evolution/Executables/GeneralizedHarmonic/EvolveGhSingleBlackHole.hpp
index 41fecb87508c..15f7daa9e3a6 100644
--- a/src/Evolution/Executables/GeneralizedHarmonic/EvolveGhSingleBlackHole.hpp
+++ b/src/Evolution/Executables/GeneralizedHarmonic/EvolveGhSingleBlackHole.hpp
@@ -177,7 +177,13 @@ struct EvolutionMetavars : public GeneralizedHarmonicTemplateBase<3, UseLts> {
   struct factory_creation
       : tt::ConformsTo<Options::protocols::FactoryCreation> {
     using factory_classes = Options::add_factory_classes<
-        typename gh_base::factory_creation::factory_classes,
+        // Restrict to monotonic time steppers in LTS to avoid control
+        // systems deadlocking.
+        tmpl::insert<
+            tmpl::erase<typename gh_base::factory_creation::factory_classes,
+                        LtsTimeStepper>,
+            tmpl::pair<LtsTimeStepper,
+                       TimeSteppers::monotonic_lts_time_steppers>>,
         tmpl::pair<Event,
                    tmpl::flatten<tmpl::list<
                        intrp::Events::Interpolate<3, ApparentHorizon,
diff --git a/src/Evolution/Executables/GrMhd/GhValenciaDivClean/GhValenciaDivCleanBase.hpp b/src/Evolution/Executables/GrMhd/GhValenciaDivClean/GhValenciaDivCleanBase.hpp
index a6bef229d6a3..b2a7faef33ff 100644
--- a/src/Evolution/Executables/GrMhd/GhValenciaDivClean/GhValenciaDivCleanBase.hpp
+++ b/src/Evolution/Executables/GrMhd/GhValenciaDivClean/GhValenciaDivCleanBase.hpp
@@ -616,7 +616,13 @@ struct GhValenciaDivCleanTemplateBase<
             boundary_conditions>,
         tmpl::pair<gh::gauges::GaugeCondition, gh::gauges::all_gauges>,
         tmpl::pair<evolution::initial_data::InitialData, initial_data_list>,
-        tmpl::pair<LtsTimeStepper, TimeSteppers::lts_time_steppers>,
+        // Restrict to monotonic time steppers in LTS to avoid control
+        // systems deadlocking.
+        tmpl::pair<
+            LtsTimeStepper,
+            tmpl::conditional_t<use_control_systems,
+                                TimeSteppers::monotonic_lts_time_steppers,
+                                TimeSteppers::lts_time_steppers>>,
         tmpl::pair<PhaseChange, PhaseControl::factory_creatable_classes>,
         tmpl::pair<StepChooser<StepChooserUse::LtsStep>,
                    StepChoosers::standard_step_choosers<system>>,
diff --git a/src/Evolution/Executables/ScalarTensor/EvolveScalarTensorSingleBlackHole.hpp b/src/Evolution/Executables/ScalarTensor/EvolveScalarTensorSingleBlackHole.hpp
index 8cef0b370e03..af38b458ff5f 100644
--- a/src/Evolution/Executables/ScalarTensor/EvolveScalarTensorSingleBlackHole.hpp
+++ b/src/Evolution/Executables/ScalarTensor/EvolveScalarTensorSingleBlackHole.hpp
@@ -154,7 +154,13 @@ struct EvolutionMetavars : public ScalarTensorTemplateBase<EvolutionMetavars> {
   struct factory_creation
       : tt::ConformsTo<Options::protocols::FactoryCreation> {
     using factory_classes = Options::add_factory_classes<
-        typename st_base::factory_creation::factory_classes,
+        // Restrict to monotonic time steppers in LTS to avoid control
+        // systems deadlocking.
+        tmpl::insert<
+            tmpl::erase<typename st_base::factory_creation::factory_classes,
+                        LtsTimeStepper>,
+            tmpl::pair<LtsTimeStepper,
+                       TimeSteppers::monotonic_lts_time_steppers>>,
         tmpl::pair<Event,
                    tmpl::flatten<tmpl::list<
                        intrp::Events::Interpolate<volume_dim, AhA,
diff --git a/src/Time/Actions/ChangeStepSize.hpp b/src/Time/Actions/ChangeStepSize.hpp
index dbe0065d206e..ff8982ace95c 100644
--- a/src/Time/Actions/ChangeStepSize.hpp
+++ b/src/Time/Actions/ChangeStepSize.hpp
@@ -114,9 +114,8 @@ bool change_step_size(const gsl::not_null<db::DataBox<DbTags>*> box) {
         << time_step_id.substep_time() << ".");
   }
 
-  const auto& next_time_id = db::get<Tags::Next<Tags::TimeStepId>>(*box);
-  const auto new_step =
-      choose_lts_step_size(next_time_id.step_time(), desired_step);
+  const auto new_step = choose_lts_step_size(
+      time_step_id.step_time() + current_step, desired_step);
   db::mutate<Tags::Next<Tags::TimeStep>>(
       [&new_step](const gsl::not_null<TimeDelta*> next_step) {
         *next_step = new_step;
diff --git a/src/Time/BoundaryHistory.hpp b/src/Time/BoundaryHistory.hpp
index 11940e4ad64f..9963470ec657 100644
--- a/src/Time/BoundaryHistory.hpp
+++ b/src/Time/BoundaryHistory.hpp
@@ -15,7 +15,6 @@
 
 #include "DataStructures/CircularDeque.hpp"
 #include "DataStructures/MathWrapper.hpp"
-#include "DataStructures/StaticDeque.hpp"
 #include "Time/History.hpp"
 #include "Time/TimeStepId.hpp"
 #include "Utilities/Algorithm.hpp"
@@ -342,20 +341,16 @@ class BoundaryHistory {
   void clear_substeps_local(size_t n);
   void clear_substeps_remote(size_t n);
 
-  StaticDeque<StepData<LocalData>, history_max_past_steps + 2> local_data_{};
+  CircularDeque<StepData<LocalData>> local_data_{};
   CircularDeque<StepData<RemoteData>> remote_data_{};
 
   template <typename Data>
   using CouplingSubsteps =
       boost::container::static_vector<Data, history_max_substeps + 1>;
 
-  // Putting the CircularDeque outermost means that we are inserting
-  // and removing containers that do not allocate, so we don't have to
-  // worry about that.
   // NOLINTNEXTLINE(spectre-mutable)
   mutable CircularDeque<CouplingSubsteps<
-      StaticDeque<CouplingSubsteps<std::optional<CouplingResult>>,
-                  decltype(local_data_)::max_size()>>>
+      CircularDeque<CouplingSubsteps<std::optional<CouplingResult>>>>>
       couplings_;
 };
 
diff --git a/src/Time/TimeSteppers/AdamsBashforth.cpp b/src/Time/TimeSteppers/AdamsBashforth.cpp
index ea59e0b401b5..e7c05f0225a4 100644
--- a/src/Time/TimeSteppers/AdamsBashforth.cpp
+++ b/src/Time/TimeSteppers/AdamsBashforth.cpp
@@ -4,15 +4,12 @@
 #include "Time/TimeSteppers/AdamsBashforth.hpp"
 
 #include <algorithm>
-#include <boost/container/small_vector.hpp>
 #include <boost/iterator/transform_iterator.hpp>
 #include <cstddef>
+#include <cstdint>
 #include <iterator>
-#include <limits>
 #include <pup.h>
-#include <utility>
 
-#include "NumericalAlgorithms/Interpolation/LagrangePolynomial.hpp"
 #include "Time/ApproximateTime.hpp"
 #include "Time/BoundaryHistory.hpp"
 #include "Time/EvolutionOrdering.hpp"
@@ -21,6 +18,7 @@
 #include "Time/Time.hpp"
 #include "Time/TimeStepId.hpp"
 #include "Time/TimeSteppers/AdamsCoefficients.hpp"
+#include "Time/TimeSteppers/AdamsLts.hpp"
 #include "Utilities/Algorithm.hpp"
 #include "Utilities/ErrorHandling/Assert.hpp"
 #include "Utilities/ErrorHandling/Error.hpp"
@@ -223,8 +221,16 @@ void AdamsBashforth::add_boundary_delta_impl(
     const TimeSteppers::ConstBoundaryHistoryTimes& remote_times,
     const TimeSteppers::BoundaryHistoryEvaluator<T>& coupling,
     const TimeDelta& time_step) const {
-  boundary_impl(result, local_times, remote_times, coupling,
-                local_times.back().step_time() + time_step);
+  const auto step_start = local_times.back().step_time();
+  const size_t integration_order =
+      local_times.integration_order(local_times.size() - 1);
+
+  const adams_lts::AdamsScheme scheme{adams_lts::SchemeType::Explicit,
+                                      integration_order};
+  const auto lts_coefficients = adams_lts::lts_coefficients(
+      local_times, remote_times, step_start, step_start + time_step, scheme,
+      scheme, scheme);
+  adams_lts::apply_coefficients(result, lts_coefficients, coupling);
 }
 
 void AdamsBashforth::clean_boundary_history_impl(
@@ -256,362 +262,20 @@ void AdamsBashforth::boundary_dense_output_impl(
     // which is the input value of `result`.
     return;
   }
-  return boundary_impl(result, local_times, remote_times, coupling,
-                       ApproximateTime{time});
-}
-
-namespace {
-template <typename T>
-class SmallStepIterator {
- public:
-  using iterator_category = std::forward_iterator_tag;
-  using value_type = Time;
-  using pointer = const Time*;
-  using reference = const Time&;
-  using difference_type = std::ptrdiff_t;
-
-  enum class Side { Local, Remote, Both };
-
-  using history_iterator = typename ConstBoundaryHistoryTimes::const_iterator;
-
-  SmallStepIterator() = default;
-
-  SmallStepIterator(history_iterator local_begin, history_iterator remote_begin,
-                    history_iterator local_end, history_iterator remote_end)
-      : local_id_(std::move(local_begin)),
-        remote_id_(std::move(remote_begin)),
-        local_end_(std::move(local_end)),
-        remote_end_(std::move(remote_end)) {}
-
-  reference operator*() const {
-    return std::max(*local_id_, *remote_id_).step_time();
-  }
-  pointer operator->() const { return &**this; }
-
-  Side side() const {
-    if (*local_id_ < *remote_id_) {
-      return Side::Remote;
-    } else if (*remote_id_ < *local_id_) {
-      return Side::Local;
-    } else {
-      return Side::Both;
-    }
-  }
-
-  // These are the m^s(n) in the paper.
-  const history_iterator& local_iterator() const { return local_id_; }
-  const history_iterator& remote_iterator() const { return remote_id_; }
-
-  SmallStepIterator& operator++() {
-    ASSERT(local_id_ != local_end_ and remote_id_ != remote_end_,
-           "Overran iterator");
-    auto local_candidate = std::next(local_id_);
-    auto remote_candidate = std::next(remote_id_);
-
-    // NOLINTNEXTLINE(bugprone-branch-clone)
-    if (local_candidate == local_end_ and remote_candidate == remote_end_) {
-      local_id_ = std::move(local_candidate);
-      remote_id_ = std::move(remote_candidate);
-      // NOLINTNEXTLINE(bugprone-branch-clone)
-    } else if (local_candidate == local_end_) {
-      remote_id_ = std::move(remote_candidate);
-      // NOLINTNEXTLINE(bugprone-branch-clone)
-    } else if (remote_candidate == remote_end_) {
-      local_id_ = std::move(local_candidate);
-    } else if (*local_candidate < *remote_candidate) {
-      local_id_ = std::move(local_candidate);
-    } else if (*remote_candidate < *local_candidate) {
-      remote_id_ = std::move(remote_candidate);
-    } else {
-      local_id_ = std::move(local_candidate);
-      remote_id_ = std::move(remote_candidate);
-    }
-
-    return *this;
-  }
-
-  bool done() const {
-    return local_id_ == local_end_ and remote_id_ == remote_end_;
-  }
-
- private:
-  history_iterator local_id_{};
-  history_iterator remote_id_{};
-  history_iterator local_end_{};
-  history_iterator remote_end_{};
-};
-
-template <typename T>
-bool operator==(const SmallStepIterator<T>& a, const SmallStepIterator<T>& b) {
-  if (a.done() and b.done()) {
-    return true;
-  }
-  if (a.done() or b.done()) {
-    return false;
-  }
-  return *a == *b;
-}
-
-template <typename T>
-bool operator!=(const SmallStepIterator<T>& a, const SmallStepIterator<T>& b) {
-  return not(a == b);
-}
-
-template <typename T>
-bool operator<(const SmallStepIterator<T>& a, const SmallStepIterator<T>& b) {
-  return a.local_iterator() < b.local_iterator() or
-         a.remote_iterator() < b.remote_iterator();
-}
-
-template <typename T>
-bool operator>(const SmallStepIterator<T>& a, const SmallStepIterator<T>& b) {
-  return b < a;
-}
-
-template <typename It>
-It bounded_next_impl(std::input_iterator_tag /*meta*/, It it, const It& bound,
-                     typename std::iterator_traits<It>::difference_type n) {
-  ASSERT(n >= 0, "Can't advance an input iterator backwards.");
-  for (typename std::iterator_traits<It>::difference_type i = 0;
-       i < n and it != bound; ++i, ++it) {
-  }
-  return it;
-}
-
-template <typename It>
-It bounded_next_impl(std::random_access_iterator_tag /*meta*/, const It& it,
-                     typename std::iterator_traits<It>::difference_type n,
-                     const It& bound) {
-  if (bound - it < n) {
-    return bound;
-  } else {
-    return it + n;
-  }
-}
-
-template <typename It>
-It bounded_next(const It& it, const It& bound, const size_t n) {
-  return bounded_next_impl(
-      typename std::iterator_traits<It>::iterator_category{}, it, bound,
-      static_cast<typename std::iterator_traits<It>::difference_type>(n));
-}
-}  // namespace
-
-template <typename T, typename TimeType>
-void AdamsBashforth::boundary_impl(
-    const gsl::not_null<T*> result,
-    const ConstBoundaryHistoryTimes& local_times,
-    const ConstBoundaryHistoryTimes& remote_times,
-    const BoundaryHistoryEvaluator<T>& coupling,
-    const TimeType& end_time) const {
-  // Might be different from order_ during self-start.
   const auto current_order =
       local_times.integration_order(local_times.size() - 1);
-
-  ASSERT(current_order <= order_,
-         "Local history is too long for target order (" << current_order
-         << " should not exceed " << order_ << ")");
-  ASSERT(remote_times.size() >= current_order,
-         "Remote history is too short (" << remote_times.size()
-         << " should be at least " << current_order << ")");
-
-  // Avoid billions of casts
-  const auto order_s = static_cast<
-      typename ConstBoundaryHistoryTimes::iterator::difference_type>(
-      current_order);
-
-  // Start and end of the step we are trying to take
-  const Time start_time = local_times.back().step_time();
-  const auto time_step = end_time - start_time;
-
-  ASSERT(local_times.size() == current_order,
-         "Incorrect local data, have " << local_times.size() << " entry, need "
-         << current_order);
-
-  if (std::equal(local_times.begin(), local_times.end(),
-                 remote_times.end() - order_s)) {
-    // No local time-stepping going on.
-    OrderVector<Time> control_points{};
-    std::transform(local_times.begin(), local_times.end(),
-                   std::back_inserter(control_points),
-                   [](const TimeStepId& t) { return t.step_time(); });
-    const auto coefficients = adams_coefficients::coefficients(
-        control_points.begin(), control_points.end(), start_time, end_time);
-
-    auto local_it = local_times.begin();
-    auto remote_it = remote_times.end() - order_s;
-    for (auto coefficients_it = coefficients.begin();
-         coefficients_it != coefficients.end();
-         ++coefficients_it, ++local_it, ++remote_it) {
-      *result += *coefficients_it * *coupling(*local_it, *remote_it);
-    }
-    return;
-  }
-
-  ASSERT(current_order == order_,
-         "Cannot perform local time-stepping while self-starting.");
-
-  const evolution_less<> less{time_step.is_positive()};
-  const auto remote_begin =
-      std::upper_bound(remote_times.begin(), remote_times.end(),
-                       local_times.back()) -
-      order_s;
-
-  ASSERT(std::is_sorted(local_times.begin(), local_times.end()),
-         "Local history not in order");
-  ASSERT(std::is_sorted(remote_begin, remote_times.end()),
-         "Remote history not in order");
-  ASSERT(not less(start_time, (remote_begin + (order_s - 1))->step_time()),
-         "Remote history does not extend far enough back");
-  ASSERT(less((remote_times.end() - 1)->step_time(), end_time),
-         "Please supply only older data: "
-         << (remote_times.end() - 1)->step_time() << " is not before "
-         << end_time);
-
-  using difference_type = std::ptrdiff_t;
-
-  SmallStepIterator<T> contributing_small_step(
-      local_times.begin(), remote_begin, local_times.end(), remote_times.end());
-  SmallStepIterator<T> small_step_of_current_step = std::next(
-      contributing_small_step, static_cast<difference_type>(current_order - 1));
-  while (*small_step_of_current_step != start_time) {
-    ++contributing_small_step;
-    ++small_step_of_current_step;
-  }
-
-  // The size of this vector is the number of small steps in the
-  // current step, which will almost always be 1 or 2, but could be
-  // larger if two elements decide to do 4:1 stepping or something.
-  boost::container::small_vector<OrderVector<double>, 2>
-      small_step_coefficients{};
-  {
-    auto coefficient_eval_begin = contributing_small_step;
-    auto coefficient_eval_end = small_step_of_current_step;
-    while (coefficient_eval_end != SmallStepIterator<T>{}) {
-      auto next_end = std::next(coefficient_eval_end);
-      const double next_time = next_end == SmallStepIterator<T>{}
-                                   ? end_time.value()
-                                   : next_end->value();
-      small_step_coefficients.push_back(adams_coefficients::coefficients(
-          coefficient_eval_begin, next_end, *coefficient_eval_end,
-          ApproximateTime{next_time}));
-      ++coefficient_eval_begin;
-      coefficient_eval_end = std::move(next_end);
-    }
-  }
-
-  // Sum over the small steps that contribute to this step, doing the
-  // appropriate interpolation for each.
-  for (size_t contributing_step_index = 0;
-       contributing_small_step != SmallStepIterator<T>{};
-       ++contributing_small_step, ++contributing_step_index) {
-    if (contributing_small_step.side() != SmallStepIterator<T>::Side::Local) {
-      double overall_prefactor = 0.0;
-      auto small_step_within_current_step = static_cast<size_t>(
-          std::max(static_cast<difference_type>(contributing_step_index + 1 -
-                                                current_order),
-                   static_cast<difference_type>(0)));
-      const size_t small_step_within_current_step_end =
-          std::min(small_step_coefficients.size(), contributing_step_index + 1);
-      for (;
-           small_step_within_current_step < small_step_within_current_step_end;
-           ++small_step_within_current_step) {
-        overall_prefactor +=
-            small_step_coefficients[small_step_within_current_step]
-                                   [contributing_step_index -
-                                    small_step_within_current_step];
-      }
-      if (contributing_small_step.side() == SmallStepIterator<T>::Side::Both) {
-        *result += overall_prefactor *
-                   *coupling(*contributing_small_step.local_iterator(),
-                             *contributing_small_step.remote_iterator());
-      } else {
-        // Side::Remote
-        OrderVector<double> past_steps(current_order);
-        std::transform(
-            local_times.end() - static_cast<difference_type>(current_order),
-            local_times.end(), past_steps.begin(),
-            [](const TimeStepId& t) { return t.step_time().value(); });
-        size_t interpolation_index = 0;
-        for (auto interpolation_time =
-                 local_times.end() -
-                 static_cast<difference_type>(current_order);
-             interpolation_time != local_times.end();
-             ++interpolation_time, ++interpolation_index) {
-          const double coefficient =
-              overall_prefactor *
-              lagrange_polynomial(interpolation_index,
-                                  contributing_small_step->value(),
-                                  past_steps.begin(), past_steps.end());
-          *result += coefficient *
-                     *coupling(*interpolation_time,
-                               *contributing_small_step.remote_iterator());
-        }
-      }
-    } else {
-      // Side::Local
-      auto interpolation_time =
-          std::max(small_step_of_current_step.remote_iterator(),
-                   contributing_small_step.remote_iterator()) -
-          static_cast<difference_type>(current_order - 1);
-      const auto interpolation_time_end =
-          bounded_next(contributing_small_step, SmallStepIterator<T>{},
-                       current_order)
-              .remote_iterator();
-      for (; interpolation_time != interpolation_time_end;
-           ++interpolation_time) {
-        double coefficient = 0.0;
-        size_t small_step_within_current_step_index = 0;
-        auto small_step_within_current_step = small_step_of_current_step;
-        while (small_step_within_current_step.remote_iterator() <
-               interpolation_time) {
-          ++small_step_within_current_step;
-          ++small_step_within_current_step_index;
-        }
-        auto small_step_within_current_step_end = contributing_small_step;
-        const auto bound_from_interpolation_time =
-            bounded_next(interpolation_time, remote_times.end(), current_order);
-        for (size_t i = 0;
-             i < current_order and
-             small_step_within_current_step_end != SmallStepIterator<T>{} and
-             small_step_within_current_step_end.remote_iterator() <
-                 bound_from_interpolation_time;
-             ++i) {
-          ++small_step_within_current_step_end;
-        }
-        auto remote_steps_end = remote_times.begin();
-        double lagrange_factor = std::numeric_limits<double>::signaling_NaN();
-        for (; small_step_within_current_step !=
-               small_step_within_current_step_end;
-             ++small_step_within_current_step,
-             ++small_step_within_current_step_index) {
-          if (remote_steps_end !=
-              small_step_within_current_step.remote_iterator() + 1) {
-            remote_steps_end =
-                small_step_within_current_step.remote_iterator() + 1;
-            OrderVector<double> past_steps(current_order);
-            std::transform(
-                remote_steps_end - static_cast<difference_type>(current_order),
-                remote_steps_end, past_steps.begin(),
-                [](const TimeStepId& t) { return t.step_time().value(); });
-            lagrange_factor = lagrange_polynomial(
-                current_order -
-                    static_cast<size_t>(remote_steps_end - interpolation_time),
-                contributing_small_step->value(), past_steps.begin(),
-                past_steps.end());
-          }
-          coefficient +=
-              lagrange_factor *
-              small_step_coefficients[small_step_within_current_step_index]
-                                     [contributing_step_index -
-                                      small_step_within_current_step_index];
-        }
-        *result +=
-            coefficient * *coupling(*contributing_small_step.local_iterator(),
-                                    *interpolation_time);
-      }
-    }
-  }
+  const adams_lts::AdamsScheme scheme{adams_lts::SchemeType::Explicit,
+                                      current_order};
+  const auto small_step_start =
+      std::max(local_times.back(), remote_times.back()).step_time();
+  const auto lts_coefficients =
+      adams_lts::lts_coefficients(local_times, remote_times,
+                                  local_times.back().step_time(),
+                                  small_step_start, scheme, scheme, scheme) +
+      adams_lts::lts_coefficients(local_times, remote_times, small_step_start,
+                                  ApproximateTime{time}, scheme, scheme,
+                                  scheme);
+  adams_lts::apply_coefficients(result, lts_coefficients, coupling);
 }
 
 bool operator==(const AdamsBashforth& lhs, const AdamsBashforth& rhs) {
diff --git a/src/Time/TimeSteppers/AdamsBashforth.hpp b/src/Time/TimeSteppers/AdamsBashforth.hpp
index 705f858a89ba..8523395557ce 100644
--- a/src/Time/TimeSteppers/AdamsBashforth.hpp
+++ b/src/Time/TimeSteppers/AdamsBashforth.hpp
@@ -304,13 +304,6 @@ class AdamsBashforth : public LtsTimeStepper {
       const TimeSteppers::BoundaryHistoryEvaluator<T>& coupling,
       const double time) const;
 
-  template <typename T, typename TimeType>
-  void boundary_impl(gsl::not_null<T*> result,
-                     const ConstBoundaryHistoryTimes& local_times,
-                     const ConstBoundaryHistoryTimes& remote_times,
-                     const BoundaryHistoryEvaluator<T>& coupling,
-                     const TimeType& end_time) const;
-
   TIME_STEPPER_DECLARE_OVERLOADS
   LTS_TIME_STEPPER_DECLARE_OVERLOADS
 
diff --git a/src/Time/TimeSteppers/AdamsLts.cpp b/src/Time/TimeSteppers/AdamsLts.cpp
new file mode 100644
index 000000000000..de44d95ace67
--- /dev/null
+++ b/src/Time/TimeSteppers/AdamsLts.cpp
@@ -0,0 +1,449 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Time/TimeSteppers/AdamsLts.hpp"
+
+#include <algorithm>
+#include <cstddef>
+#include <iterator>
+#include <type_traits>
+#include <utility>
+
+#include "DataStructures/MathWrapper.hpp"
+#include "NumericalAlgorithms/Interpolation/LagrangePolynomial.hpp"
+#include "Time/ApproximateTime.hpp"
+#include "Time/BoundaryHistory.hpp"
+#include "Time/EvolutionOrdering.hpp"
+#include "Time/Time.hpp"
+#include "Time/TimeStepId.hpp"
+#include "Time/TimeSteppers/AdamsCoefficients.hpp"
+#include "Utilities/Algorithm.hpp"
+#include "Utilities/ErrorHandling/Assert.hpp"
+#include "Utilities/ErrorHandling/Error.hpp"
+#include "Utilities/GenerateInstantiations.hpp"
+#include "Utilities/Gsl.hpp"
+
+namespace TimeSteppers::adams_lts {
+Time exact_substep_time(const TimeStepId& id) {
+  const Time result =
+      id.substep() == 0 ? id.step_time() : id.step_time() + id.step_size();
+  ASSERT(result.value() == id.substep_time(), "Substep not at expected time.");
+  return result;
+}
+
+namespace {
+template <typename T>
+using OrderVector = adams_coefficients::OrderVector<T>;
+
+template <typename Op>
+LtsCoefficients& add_assign_impl(LtsCoefficients& a, const LtsCoefficients& b,
+                                 const Op& op) {
+  const auto key_equal = [](const LtsCoefficients::value_type& l,
+                            const LtsCoefficients::value_type& r) {
+    return get<0>(l) == get<0>(r) and get<1>(l) == get<1>(r);
+  };
+  const auto key_less = [](const LtsCoefficients::value_type& l,
+                           const LtsCoefficients::value_type& r) {
+    return get<0>(l) < get<0>(r) or
+           (get<0>(l) == get<0>(r) and get<1>(l) < get<1>(r));
+  };
+
+  // Two passes: first the common entries, and then the new entries.
+  size_t common_entries = 0;
+  {
+    auto a_it = a.begin();
+    auto b_it = b.begin();
+    while (a_it != a.end() and b_it != b.end()) {
+      if (key_less(*a_it, *b_it)) {
+        ++a_it;
+      } else if (key_less(*b_it, *a_it)) {
+        ++b_it;
+      } else {
+        ++common_entries;
+        get<2>(*a_it) += op(get<2>(*b_it));
+        ++a_it;
+        ++b_it;
+      }
+    }
+  }
+  a.resize(a.size() + b.size() - common_entries);
+  {
+    auto write = a.rbegin();
+    auto a_it = write + static_cast<LtsCoefficients::difference_type>(
+                            b.size() - common_entries);
+    auto b_it = b.rbegin();
+    while (b_it != b.rend()) {
+      if (a_it == a.rend() or key_less(*a_it, *b_it)) {
+        *write = *b_it;
+        get<2>(*write) = op(get<2>(*write));
+        ++b_it;
+      } else {
+        if (key_equal(*a_it, *b_it)) {
+          ++b_it;
+        }
+        *write = *a_it;
+        ++a_it;
+      }
+      ++write;
+    }
+  }
+  return a;
+}
+}  // namespace
+
+LtsCoefficients& operator+=(LtsCoefficients& a, const LtsCoefficients& b) {
+  return add_assign_impl(a, b, [](const double x) { return x; });
+}
+LtsCoefficients& operator-=(LtsCoefficients& a, const LtsCoefficients& b) {
+  return add_assign_impl(a, b, [](const double x) { return -x; });
+}
+
+LtsCoefficients operator+(LtsCoefficients&& a, LtsCoefficients&& b) {
+  return std::move(a += b);
+}
+LtsCoefficients operator+(LtsCoefficients&& a, const LtsCoefficients& b) {
+  return std::move(a += b);
+}
+LtsCoefficients operator+(const LtsCoefficients& a, LtsCoefficients&& b) {
+  return std::move(b += a);
+}
+LtsCoefficients operator+(const LtsCoefficients& a, const LtsCoefficients& b) {
+  auto a2 = a;
+  return std::move(a2 += b);
+}
+
+LtsCoefficients operator-(LtsCoefficients&& a, LtsCoefficients&& b) {
+  return std::move(a -= b);
+}
+LtsCoefficients operator-(LtsCoefficients&& a, const LtsCoefficients& b) {
+  return std::move(a -= b);
+}
+LtsCoefficients operator-(const LtsCoefficients& a, const LtsCoefficients& b) {
+  auto a2 = a;
+  return std::move(a2 -= b);
+}
+LtsCoefficients operator-(const LtsCoefficients& a, LtsCoefficients&& b) {
+  alg::for_each(b, [](auto& coef) { get<2>(coef) *= -1.0; });
+  return std::move(b += a);
+}
+
+template <typename T>
+void apply_coefficients(const gsl::not_null<T*> result,
+                        const LtsCoefficients& coefficients,
+                        const BoundaryHistoryEvaluator<T>& coupling) {
+  for (const auto& term : coefficients) {
+    *result += get<2>(term) * *coupling(get<0>(term), get<1>(term));
+  }
+}
+
+bool operator==(const AdamsScheme& a, const AdamsScheme& b) {
+  return a.type == b.type and a.order == b.order;
+}
+bool operator!=(const AdamsScheme& a, const AdamsScheme& b) {
+  return not(a == b);
+}
+
+namespace {
+// Collect the ids used for interpolating during a step to `end_time`
+// from `times`.
+//
+// For implicit schemes, the step containing (or ending at) `end_time`
+// must have predictor data available.
+template <typename TimeType>
+OrderVector<TimeStepId> find_relevant_ids(
+    const ConstBoundaryHistoryTimes& times, const TimeType& end_time,
+    const AdamsScheme& scheme) {
+  OrderVector<TimeStepId> ids{};
+  using difference_type = std::iterator_traits<
+      ConstBoundaryHistoryTimes::const_iterator>::difference_type;
+  const evolution_less<> less{times.front().time_runs_forward()};
+  // Can't do a binary search because times are not sorted during self-start.
+  auto used_range_end = times.end();
+  while (used_range_end != times.begin()) {
+    if (less((used_range_end - 1)->step_time(), end_time)) {
+      break;
+    }
+    --used_range_end;
+  }
+  const auto number_of_past_steps =
+      static_cast<difference_type>(scheme.order) -
+      (scheme.type == SchemeType::Implicit ? 1 : 0);
+  ASSERT(used_range_end - times.begin() >= number_of_past_steps,
+         "Insufficient past data.");
+  std::copy(used_range_end - number_of_past_steps, used_range_end,
+            std::back_inserter(ids));
+  if (scheme.type == SchemeType::Implicit) {
+    const auto last_step =
+        static_cast<size_t>(used_range_end - times.begin() - 1);
+    ASSERT(times.number_of_substeps(last_step) == 2,
+           "Must have substep data for implicit stepping.");
+    ids.push_back(times[{last_step, 1}]);
+  }
+  return ids;
+}
+
+// Choose the relevant times from `local` and `remote` for defining
+// small steps for `small_step_scheme`, using the specified schemes
+// for interpolation on the local and remote sides.
+//
+// This is the most recent values from the union of the `local` and
+// `remote` times with any values that should only be used for
+// interpolation removed.
+OrderVector<Time> merge_to_small_steps(const OrderVector<Time>& local,
+                                       const OrderVector<Time>& remote,
+                                       const evolution_less<Time>& less,
+                                       const AdamsScheme& local_scheme,
+                                       const AdamsScheme& remote_scheme,
+                                       const AdamsScheme& small_step_scheme) {
+  OrderVector<Time> small_steps(small_step_scheme.order);
+  auto local_it = local.rbegin();
+  auto remote_it = remote.rbegin();
+
+  ASSERT(not(local_scheme.type == SchemeType::Implicit and
+             remote_scheme.type == SchemeType::Explicit and
+             less(*local_it, *remote_it)),
+         "Explicit time " << *remote_it << " after implicit " << *local_it);
+  ASSERT(not(remote_scheme.type == SchemeType::Implicit and
+             local_scheme.type == SchemeType::Explicit and
+             less(*remote_it, *local_it)),
+         "Explicit time " << *local_it << " after implicit " << *remote_it);
+
+  if (small_step_scheme.type == SchemeType::Explicit) {
+    // Don't use implicit interpolation points for an explicit step
+    if (local_scheme.type == SchemeType::Implicit) {
+      ++local_it;
+    }
+    if (remote_scheme.type == SchemeType::Implicit) {
+      ++remote_it;
+    }
+  } else {
+    if (local_scheme.type == SchemeType::Implicit and
+        remote_scheme.type == SchemeType::Implicit) {
+      // If both the interpolation schemes are implicit, we will get
+      // two times after the small step we are working on, one from
+      // each.  One of them (if they are different) belongs to a later
+      // small step, and we should ignore it.  If they are the same,
+      // ignoring one of them is harmless.
+      if (less(*local_it, *remote_it)) {
+        ++remote_it;
+      } else {
+        ++local_it;
+      }
+    }
+  }
+
+  for (auto out = small_steps.rbegin(); out != small_steps.rend(); ++out) {
+    if (local_it == local.rend()) {
+      ASSERT(remote_it != remote.rend(), "Ran out of data");
+      *out = *remote_it;
+      ++remote_it;
+    } else if (remote_it == remote.rend()) {
+      *out = *local_it;
+      ++local_it;
+    } else {
+      *out = std::max(*local_it, *remote_it, less);
+      if (*local_it == *out) {
+        ++local_it;
+      }
+      if (*remote_it == *out) {
+        ++remote_it;
+      }
+    }
+  }
+  return small_steps;
+}
+
+// Evaluate the Lagrange interpolating polynomials with the given
+// `control_times` at `time`.  The returned vector contains the values
+// of all the Lagrange polynomials.
+template <typename TimeType>
+OrderVector<double> interpolation_coefficients(
+    const OrderVector<Time>& control_times, const TimeType& time) {
+  if constexpr (std::is_same_v<TimeType, Time>) {
+    // Skip the Lagrange polynomial calculations if we are evaluating
+    // at a control time.  This should be common.
+    for (size_t i = 0; i < control_times.size(); ++i) {
+      if (control_times[i] == time) {
+        OrderVector<double> coefficients(control_times.size(), 0.0);
+        coefficients[i] = 1.0;
+        return coefficients;
+      }
+    }
+  }
+
+  OrderVector<double> control_times_fp(control_times.size());
+  alg::transform(control_times, control_times_fp.begin(),
+                 [](const Time& t) { return t.value(); });
+
+  OrderVector<double> coefficients{};
+  for (size_t i = 0; i < control_times.size(); ++i) {
+    coefficients.push_back(lagrange_polynomial(
+        i, time.value(), control_times_fp.begin(), control_times_fp.end()));
+  }
+  return coefficients;
+}
+
+template <typename TimeType>
+LtsCoefficients lts_coefficients_for_gts(
+    const OrderVector<TimeStepId>& control_ids, const Time& start_time,
+    const TimeType& end_time) {
+  // The sides are stepping at the same rate, so no LTS is happening
+  // at this boundary.
+  OrderVector<Time> control_times(control_ids.size());
+  alg::transform(control_ids, control_times.begin(), exact_substep_time);
+
+  const OrderVector<double> gts_coefficients = adams_coefficients::coefficients(
+      control_times.begin(), control_times.end(), start_time, end_time);
+  LtsCoefficients lts_coefficients{};
+  for (size_t step = 0; step < gts_coefficients.size(); ++step) {
+    lts_coefficients.emplace_back(control_ids[step], control_ids[step],
+                                  gts_coefficients[step]);
+  }
+  return lts_coefficients;
+}
+}  // namespace
+
+template <typename TimeType>
+LtsCoefficients lts_coefficients(const ConstBoundaryHistoryTimes& local_times,
+                                 const ConstBoundaryHistoryTimes& remote_times,
+                                 const Time& start_time,
+                                 const TimeType& end_time,
+                                 const AdamsScheme& local_scheme,
+                                 const AdamsScheme& remote_scheme,
+                                 const AdamsScheme& small_step_scheme) {
+  if (start_time == end_time) {
+    return {};
+  }
+  const evolution_less<Time> time_less{local_times.front().time_runs_forward()};
+
+  LtsCoefficients step_coefficients{};
+
+  TimeType small_step_end = end_time;
+  for (;;) {
+    const OrderVector<TimeStepId> local_ids =
+        find_relevant_ids(local_times, small_step_end, local_scheme);
+    const OrderVector<TimeStepId> remote_ids =
+        find_relevant_ids(remote_times, small_step_end, remote_scheme);
+
+    // Check is the there is actually local time-stepping happening at
+    // this boundary.  Only check for the latest small step, before we
+    // have generated any coefficients.
+    if (step_coefficients.empty() and small_step_scheme == local_scheme and
+        small_step_scheme == remote_scheme and local_ids == remote_ids) {
+      return lts_coefficients_for_gts(local_ids, start_time, end_time);
+    }
+
+    OrderVector<Time> local_control_times(local_scheme.order);
+    alg::transform(local_ids, local_control_times.begin(), exact_substep_time);
+    OrderVector<Time> remote_control_times(remote_scheme.order);
+    alg::transform(remote_ids, remote_control_times.begin(),
+                   exact_substep_time);
+
+    const OrderVector<Time> small_step_times = merge_to_small_steps(
+        local_control_times, remote_control_times, time_less, local_scheme,
+        remote_scheme, small_step_scheme);
+    const Time current_small_step =
+        small_step_times[small_step_times.size() -
+                         (small_step_scheme.type == SchemeType::Implicit ? 2
+                                                                         : 1)];
+    ASSERT(not time_less(current_small_step, start_time),
+           "Reached time " << current_small_step
+           << " without hitting start time " << start_time
+           << " while iterating over small steps.  Most likely the supplied "
+              "start time was not a step boundary.");
+
+    const OrderVector<double> small_step_coefficients =
+        adams_coefficients::coefficients(small_step_times.begin(),
+                                         small_step_times.end(),
+                                         current_small_step, small_step_end);
+
+    for (size_t contributing_small_step = 0;
+         contributing_small_step < small_step_times.size();
+         ++contributing_small_step) {
+      const OrderVector<double> local_interpolation_coefficients =
+          interpolation_coefficients(local_control_times,
+                                     small_step_times[contributing_small_step]);
+      const OrderVector<double> remote_interpolation_coefficients =
+          interpolation_coefficients(remote_control_times,
+                                     small_step_times[contributing_small_step]);
+      for (size_t local_step_index = 0;
+           local_step_index < local_interpolation_coefficients.size();
+           ++local_step_index) {
+        if (local_interpolation_coefficients[local_step_index] == 0.0) {
+          continue;
+        }
+        for (size_t remote_step_index = 0;
+             remote_step_index < remote_interpolation_coefficients.size();
+             ++remote_step_index) {
+          if (remote_interpolation_coefficients[remote_step_index] == 0.0) {
+            continue;
+          }
+          step_coefficients.emplace_back(
+              local_ids[local_step_index], remote_ids[remote_step_index],
+              small_step_coefficients[contributing_small_step] *
+                  local_interpolation_coefficients[local_step_index] *
+                  remote_interpolation_coefficients[remote_step_index]);
+        }
+      }
+    }
+    if (current_small_step == start_time) {
+      break;
+    }
+
+    if constexpr (std::is_same_v<TimeType, Time>) {
+      // We're iterating backwards, so the next step temporally is the
+      // one we just did.
+      small_step_end = current_small_step;
+    } else {
+      ERROR(
+          "Multiple-small-step dense output is not supported.  Calculate "
+          "coefficients as the sum of the exact update for the complete small "
+          "steps and a dense update for a single partial small step.");
+    }
+  }
+
+  // Combine duplicate entries
+  ASSERT(not step_coefficients.empty(),
+         "Generated no coefficients.  This an algorithmic bug.  Please file "
+         "an issue.");
+  alg::sort(step_coefficients);
+  auto unique_entry = step_coefficients.begin();
+  for (auto generated_entry = std::next(step_coefficients.begin());
+       generated_entry != step_coefficients.end();
+       ++generated_entry) {
+    if (get<0>(*generated_entry) == get<0>(*unique_entry) and
+        get<1>(*generated_entry) == get<1>(*unique_entry)) {
+      get<2>(*unique_entry) += get<2>(*generated_entry);
+    } else {
+      *++unique_entry = *generated_entry;
+    }
+  }
+  step_coefficients.erase(std::next(unique_entry), step_coefficients.end());
+  return step_coefficients;
+}
+
+#define MATH_WRAPPER_TYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
+
+#define INSTANTIATE(_, data)                          \
+  template void apply_coefficients(                   \
+      gsl::not_null<MATH_WRAPPER_TYPE(data)*> result, \
+      const LtsCoefficients& coefficients,            \
+      const BoundaryHistoryEvaluator<MATH_WRAPPER_TYPE(data)>& coupling);
+
+GENERATE_INSTANTIATIONS(INSTANTIATE, (MATH_WRAPPER_TYPES))
+#undef INSTANTIATE
+#undef MATH_WRAPPER_TYPE
+
+#define TIME_TYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
+
+#define INSTANTIATE(_, data)                                                 \
+  template LtsCoefficients lts_coefficients(                                 \
+      const ConstBoundaryHistoryTimes& local_times,                          \
+      const ConstBoundaryHistoryTimes& remote_times, const Time& start_time, \
+      const TIME_TYPE(data) & end_time, const AdamsScheme& local_scheme,     \
+      const AdamsScheme& remote_scheme, const AdamsScheme& small_step_scheme);
+
+GENERATE_INSTANTIATIONS(INSTANTIATE, (Time, ApproximateTime))
+#undef INSTANTIATE
+#undef TIME_TYPE
+}  // namespace TimeSteppers::adams_lts
diff --git a/src/Time/TimeSteppers/AdamsLts.hpp b/src/Time/TimeSteppers/AdamsLts.hpp
new file mode 100644
index 000000000000..eb36c299f32d
--- /dev/null
+++ b/src/Time/TimeSteppers/AdamsLts.hpp
@@ -0,0 +1,128 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <boost/container/small_vector.hpp>
+#include <cstddef>
+#include <tuple>
+
+#include "Time/TimeStepId.hpp"
+#include "Time/TimeSteppers/AdamsCoefficients.hpp"
+
+/// \cond
+class Time;
+namespace TimeSteppers {
+template <typename T>
+class BoundaryHistoryEvaluator;
+class ConstBoundaryHistoryTimes;
+}  // namespace TimeSteppers
+namespace gsl {
+template <class T>
+class not_null;
+}  // namespace gsl
+/// \endcond
+
+/// Shared LTS implementation for the two Adams-based methods.
+namespace TimeSteppers::adams_lts {
+/// Get the time of a substep as a `Time`, assuming substeps are at
+/// the end of the step.
+Time exact_substep_time(const TimeStepId& id);
+
+// For order-k 2:1 stepping, in each small step, half the points will
+// require interpolation (k entries each) and the others will not (1
+// entry each).  So (2 small steps) * ((k/2 interpolations) * (from k
+// points) + (k/2 non-interpolations)) = k (k + 1).  (The small steps
+// round k/2 in different directions and the effect cancels out.)
+constexpr size_t lts_coefficients_static_size =
+    adams_coefficients::maximum_order * (adams_coefficients::maximum_order + 1);
+
+/// Storage for LTS coefficients that should not allocate in typical
+/// cases.  Each entry is a tuple of (local id, remote id,
+/// coefficient).  The contents should be kept sorted, as some
+/// functions assume that.
+struct LtsCoefficients
+    : boost::container::small_vector<std::tuple<TimeStepId, TimeStepId, double>,
+                                     lts_coefficients_static_size> {
+  using boost::container::small_vector<
+      std::tuple<TimeStepId, TimeStepId, double>,
+      lts_coefficients_static_size>::small_vector;
+};
+
+LtsCoefficients& operator+=(LtsCoefficients& a, const LtsCoefficients& b);
+LtsCoefficients& operator-=(LtsCoefficients& a, const LtsCoefficients& b);
+
+LtsCoefficients operator+(LtsCoefficients&& a, LtsCoefficients&& b);
+LtsCoefficients operator+(LtsCoefficients&& a, const LtsCoefficients& b);
+LtsCoefficients operator+(const LtsCoefficients& a, LtsCoefficients&& b);
+LtsCoefficients operator+(const LtsCoefficients& a, const LtsCoefficients& b);
+
+LtsCoefficients operator-(LtsCoefficients&& a, LtsCoefficients&& b);
+LtsCoefficients operator-(LtsCoefficients&& a, const LtsCoefficients& b);
+LtsCoefficients operator-(const LtsCoefficients& a, LtsCoefficients&& b);
+LtsCoefficients operator-(const LtsCoefficients& a, const LtsCoefficients& b);
+
+/// Add the LTS boundary terms for to \p result for the given set of
+/// coefficients.
+template <typename T>
+void apply_coefficients(gsl::not_null<T*> result,
+                        const LtsCoefficients& coefficients,
+                        const BoundaryHistoryEvaluator<T>& coupling);
+
+/// Type of coefficients used for an AdamsScheme
+enum class SchemeType { Explicit, Implicit };
+
+struct AdamsScheme {
+  SchemeType type;
+  size_t order;
+};
+
+bool operator==(const AdamsScheme& a, const AdamsScheme& b);
+bool operator!=(const AdamsScheme& a, const AdamsScheme& b);
+
+/*!
+ * Calculate the nonzero terms in an Adams LTS boundary contribution.
+ *
+ * The coefficients are generated for a step from \p start_time to \p
+ * end_time, integrating using the \p small_step_scheme.
+ * Interpolation from the elements is performed using polynomials
+ * appropriate for the given \p local_scheme and \p remote_scheme.
+ *
+ * This function returns the sum of the coefficients for the small
+ * steps, i.e., the steps between all the times either side updates
+ * between the time bounds.  Each small step is calculated as
+ *
+ * \f{equation}
+ * \Delta \tilde{y}_n = \Delta \tilde{t}_n \sum_{ij} D(y^L_{n,-i}, y^R_{n,-j})
+ * \sum_q \tilde{\alpha}_{nq}
+ * \ell_i(\tilde{t}_{n-q}; t^L_\cdots, \ldots)
+ * \ell_j(\tilde{t}_{n-q}; t^R_\cdots, \ldots),
+ * \f}
+ *
+ * where $y^L_{n,0}$ and $y^R_{n,0}$ are the values of the local and
+ * remote variables at the start of the step containing small step
+ * $n$, $y^L_{n,1}$ and $y^R_{n,1}$ are the predictor values for those
+ * steps, and other index values count backwards in the history of
+ * that side.  The range of $i$ and choice of the control times
+ * $t^L_\cdots$ are determined by \p local_scheme, $j$ and
+ * $t^R_\cdots$ are determined by \p remote_scheme, and the Adams
+ * coefficients $\tilde{\alpha}_{nq}$ correspond to \p
+ * small_step_scheme.
+ *
+ * When called for dense output, the arguments must represent a
+ * single small step, i.e., the step cannot cross any of the control
+ * times.  Any additional terms can be generated by a second call
+ * treating the remainder of the step as non-dense.
+ *
+ * \tparam TimeType The type `Time` for a step aligned with the
+ * control times or `ApproximateTime` for dense output.
+ */
+template <typename TimeType>
+LtsCoefficients lts_coefficients(const ConstBoundaryHistoryTimes& local_times,
+                                 const ConstBoundaryHistoryTimes& remote_times,
+                                 const Time& start_time,
+                                 const TimeType& end_time,
+                                 const AdamsScheme& local_scheme,
+                                 const AdamsScheme& remote_scheme,
+                                 const AdamsScheme& small_step_scheme);
+}  // namespace TimeSteppers::adams_lts
diff --git a/src/Time/TimeSteppers/AdamsMoultonPc.cpp b/src/Time/TimeSteppers/AdamsMoultonPc.cpp
index 57c70e9b2d5d..2b5f323e99a9 100644
--- a/src/Time/TimeSteppers/AdamsMoultonPc.cpp
+++ b/src/Time/TimeSteppers/AdamsMoultonPc.cpp
@@ -13,6 +13,7 @@
 #include "Time/History.hpp"
 #include "Time/SelfStart.hpp"
 #include "Time/TimeSteppers/AdamsCoefficients.hpp"
+#include "Time/TimeSteppers/AdamsLts.hpp"
 #include "Utilities/Algorithm.hpp"
 #include "Utilities/ErrorHandling/Assert.hpp"
 #include "Utilities/ErrorHandling/Error.hpp"
@@ -145,9 +146,53 @@ TimeStepId AdamsMoultonPc<Monotonic>::next_time_id_for_error(
   return next_time_id(current_id, time_step);
 }
 
+template <bool Monotonic>
+bool AdamsMoultonPc<Monotonic>::neighbor_data_required(
+    const TimeStepId& next_substep_id,
+    const TimeStepId& neighbor_data_id) const {
+  // Because of self-start, step times may not be monotonic across
+  // slabs, so check that first.
+  if (neighbor_data_id.slab_number() != next_substep_id.slab_number()) {
+    return neighbor_data_id.slab_number() < next_substep_id.slab_number();
+  }
+
+  const evolution_less<Time> before{neighbor_data_id.time_runs_forward()};
+
+  if constexpr (Monotonic) {
+    const auto next_time = adams_lts::exact_substep_time(next_substep_id);
+    const auto neighbor_time = adams_lts::exact_substep_time(neighbor_data_id);
+    return before(neighbor_time, next_time) or
+           (neighbor_time == next_time and neighbor_data_id.substep() == 1 and
+            next_substep_id.substep() == 0);
+  } else {
+    if (next_substep_id.substep() == 1) {
+      // predictor
+      return before(neighbor_data_id.step_time(),
+                    next_substep_id.step_time()) or
+             (neighbor_data_id.step_time() == next_substep_id.step_time() and
+              neighbor_data_id.substep() == 0);
+    } else {
+      // corrector
+      return before(neighbor_data_id.step_time(), next_substep_id.step_time());
+    }
+  }
+}
+
+template <bool Monotonic>
+bool AdamsMoultonPc<Monotonic>::neighbor_data_required(
+    const double dense_output_time, const TimeStepId& neighbor_data_id) const {
+  const evolution_less<double> before{neighbor_data_id.time_runs_forward()};
+  if constexpr (Monotonic) {
+    return not before(dense_output_time,
+                      adams_lts::exact_substep_time(neighbor_data_id).value());
+  } else {
+    return before(neighbor_data_id.step_time().value(), dense_output_time);
+  }
+}
+
 template <bool Monotonic>
 void AdamsMoultonPc<Monotonic>::pup(PUP::er& p) {
-  TimeStepper::pup(p);
+  LtsTimeStepper::pup(p);
   p | order_;
 }
 
@@ -245,6 +290,206 @@ bool AdamsMoultonPc<Monotonic>::can_change_step_size_impl(
          });
 }
 
+template <bool Monotonic>
+template <typename T>
+void AdamsMoultonPc<Monotonic>::add_boundary_delta_impl(
+    const gsl::not_null<T*> result,
+    const TimeSteppers::ConstBoundaryHistoryTimes& local_times,
+    const TimeSteppers::ConstBoundaryHistoryTimes& remote_times,
+    const TimeSteppers::BoundaryHistoryEvaluator<T>& coupling,
+    const TimeDelta& time_step) const {
+  ASSERT(not local_times.empty(), "No local data provided.");
+  ASSERT(not remote_times.empty(), "No remote data provided.");
+  const auto current_order =
+      local_times.integration_order(local_times.size() - 1);
+  if constexpr (Monotonic) {
+    const adams_lts::AdamsScheme predictor_scheme{
+        adams_lts::SchemeType::Explicit, current_order - 1};
+    const adams_lts::AdamsScheme corrector_scheme{
+        adams_lts::SchemeType::Implicit, current_order};
+    const auto step_start = local_times.back().step_time();
+    const auto step_end = step_start + time_step;
+    const auto small_step_start =
+        std::max(local_times.back(), remote_times.back()).step_time();
+    const auto synchronization_time =
+        std::min(local_times.back(), remote_times.back()).step_time();
+    const auto is_synchronization_time = [&](const TimeStepId& id) {
+      return id.step_time() == synchronization_time;
+    };
+    ASSERT(alg::any_of(local_times, is_synchronization_time) and
+               alg::any_of(remote_times, is_synchronization_time),
+           "Only nested step patterns (N:1) are supported.");
+
+    if (local_times.number_of_substeps(local_times.size() - 1) == 1) {
+      // Predictor
+      auto lts_coefficients = adams_lts::lts_coefficients(
+          local_times, remote_times, small_step_start, step_end,
+          predictor_scheme, predictor_scheme, predictor_scheme);
+      if (not is_synchronization_time(remote_times.back())) {
+        lts_coefficients += adams_lts::lts_coefficients(
+            local_times, remote_times, synchronization_time, small_step_start,
+            predictor_scheme, corrector_scheme, corrector_scheme);
+      }
+      adams_lts::apply_coefficients(result, lts_coefficients, coupling);
+    } else {
+      // Corrector
+      if (remote_times.number_of_substeps(remote_times.size() - 1) == 2) {
+        // Aligned corrector
+        ASSERT(adams_lts::exact_substep_time(
+                   remote_times[{remote_times.size() - 1, 1}]) == step_end,
+               "Have remote substep data, but it isn't aligned with the local "
+               "data.");
+        const auto lts_coefficients =
+            adams_lts::lts_coefficients(
+                local_times, remote_times, synchronization_time, step_end,
+                corrector_scheme, corrector_scheme, corrector_scheme) -
+            adams_lts::lts_coefficients(
+                local_times, remote_times, synchronization_time, step_start,
+                corrector_scheme, predictor_scheme, corrector_scheme);
+        adams_lts::apply_coefficients(result, lts_coefficients, coupling);
+      } else {
+        // Unaligned corrector
+        ASSERT(step_start == small_step_start,
+               "Trying to take unaligned step, but remote side is smaller.");
+        const auto lts_coefficients = adams_lts::lts_coefficients(
+            local_times, remote_times, step_start, step_end, corrector_scheme,
+            predictor_scheme, corrector_scheme);
+        adams_lts::apply_coefficients(result, lts_coefficients, coupling);
+      }
+    }
+  } else {
+    adams_lts::AdamsScheme scheme{adams_lts::SchemeType::Implicit,
+                                  current_order};
+    auto remote_scheme = scheme;
+
+    if (local_times.number_of_substeps(local_times.size() - 1) == 1) {
+      // Predictor
+      scheme = {adams_lts::SchemeType::Explicit, current_order - 1};
+      ASSERT(remote_times.back() <= local_times.back(),
+             "Unexpected remote values available.");
+      // If the sides are not aligned, we use the predictor data
+      // available from the neighbor.  If they are, that data has not
+      // been received.
+      remote_scheme = {remote_times.back() == local_times.back()
+                           ? adams_lts::SchemeType::Explicit
+                           : adams_lts::SchemeType::Implicit,
+                       current_order - 1};
+    }
+
+    const auto lts_coefficients = adams_lts::lts_coefficients(
+        local_times, remote_times, local_times.back().step_time(),
+        local_times.back().step_time() + time_step, scheme, remote_scheme,
+        scheme);
+    adams_lts::apply_coefficients(result, lts_coefficients, coupling);
+  }
+}
+
+template <bool Monotonic>
+void AdamsMoultonPc<Monotonic>::clean_boundary_history_impl(
+    const TimeSteppers::MutableBoundaryHistoryTimes& local_times,
+    const TimeSteppers::MutableBoundaryHistoryTimes& remote_times) const {
+  if (local_times.empty() or
+      local_times.number_of_substeps(local_times.size() - 1) != 2) {
+    return;
+  }
+  ASSERT(not remote_times.empty(), "No remote data available.");
+
+  const bool synchronized =
+      remote_times.number_of_substeps(remote_times.size() - 1) == 2 and
+      adams_lts::exact_substep_time(local_times[{local_times.size() - 1, 1}]) ==
+          adams_lts::exact_substep_time(
+              remote_times[{remote_times.size() - 1, 1}]);
+
+  if (Monotonic and not synchronized) {
+    return;
+  }
+
+  const auto required_points =
+      local_times.integration_order(local_times.size() - 1) - 2;
+
+  while (local_times.size() > required_points) {
+    local_times.pop_front();
+  }
+  for (size_t i = 0; i < local_times.size(); ++i) {
+    local_times.clear_substeps(i);
+  }
+
+  // If the sides are not aligned, then we are in the middle of the
+  // remote step, so still need its data.
+  if (synchronized) {
+    while (remote_times.size() > required_points) {
+      remote_times.pop_front();
+    }
+    for (size_t i = 0; i < remote_times.size(); ++i) {
+      remote_times.clear_substeps(i);
+    }
+  }
+}
+
+template <bool Monotonic>
+template <typename T>
+void AdamsMoultonPc<Monotonic>::boundary_dense_output_impl(
+    const gsl::not_null<T*> result,
+    const TimeSteppers::ConstBoundaryHistoryTimes& local_times,
+    const TimeSteppers::ConstBoundaryHistoryTimes& remote_times,
+    const TimeSteppers::BoundaryHistoryEvaluator<T>& coupling,
+    const double time) const {
+  if constexpr (Monotonic) {
+    ASSERT(local_times.number_of_substeps(local_times.size() - 1) == 1,
+           "Dense output must be done before predictor evaluation.");
+
+    const auto current_order =
+        local_times.integration_order(local_times.size() - 1);
+    const adams_lts::AdamsScheme predictor_scheme{
+        adams_lts::SchemeType::Explicit, current_order - 1};
+    const adams_lts::AdamsScheme corrector_scheme{
+        adams_lts::SchemeType::Implicit, current_order};
+    const auto small_step_start =
+        std::max(local_times.back(), remote_times.back()).step_time();
+    const auto synchronization_time =
+        std::min(local_times.back(), remote_times.back()).step_time();
+    const auto is_synchronization_time = [&](const TimeStepId& id) {
+      return id.step_time() == synchronization_time;
+    };
+    ASSERT(alg::any_of(local_times, is_synchronization_time) and
+               alg::any_of(remote_times, is_synchronization_time),
+           "Only nested step patterns (N:1) are supported.");
+
+    auto lts_coefficients = adams_lts::lts_coefficients(
+        local_times, remote_times, small_step_start, ApproximateTime{time},
+        predictor_scheme, predictor_scheme, predictor_scheme);
+    if (not is_synchronization_time(remote_times.back())) {
+      lts_coefficients += adams_lts::lts_coefficients(
+          local_times, remote_times, synchronization_time, small_step_start,
+          predictor_scheme, corrector_scheme, corrector_scheme);
+    }
+    adams_lts::apply_coefficients(result, lts_coefficients, coupling);
+  } else {
+    if (local_times.back().step_time().value() == time) {
+      // Nothing to do.  The requested time is the start of the step,
+      // which is the input value of `result`.
+      return;
+    }
+    ASSERT(local_times.number_of_substeps(local_times.size() - 1) == 2,
+           "Dense output must be done after predictor evaluation.");
+
+    const auto current_order =
+        local_times.integration_order(local_times.size() - 1);
+    const adams_lts::AdamsScheme scheme{adams_lts::SchemeType::Implicit,
+                                        current_order};
+    const auto small_step_start =
+        std::max(local_times.back(), remote_times.back()).step_time();
+    const auto lts_coefficients =
+        adams_lts::lts_coefficients(local_times, remote_times,
+                                    local_times.back().step_time(),
+                                    small_step_start, scheme, scheme, scheme) +
+        adams_lts::lts_coefficients(local_times, remote_times, small_step_start,
+                                    ApproximateTime{time}, scheme, scheme,
+                                    scheme);
+    adams_lts::apply_coefficients(result, lts_coefficients, coupling);
+  }
+}
+
 template <bool Monotonic>
 bool operator==(const AdamsMoultonPc<Monotonic>& lhs,
                 const AdamsMoultonPc<Monotonic>& rhs) {
@@ -259,6 +504,8 @@ bool operator!=(const AdamsMoultonPc<Monotonic>& lhs,
 
 TIME_STEPPER_DEFINE_OVERLOADS_TEMPLATED(AdamsMoultonPc<Monotonic>,
                                         bool Monotonic)
+LTS_TIME_STEPPER_DEFINE_OVERLOADS_TEMPLATED(AdamsMoultonPc<Monotonic>,
+                                            bool Monotonic)
 
 template <bool Monotonic>
 PUP::able::PUP_ID AdamsMoultonPc<Monotonic>::my_PUP_ID = 0;  // NOLINT
diff --git a/src/Time/TimeSteppers/AdamsMoultonPc.hpp b/src/Time/TimeSteppers/AdamsMoultonPc.hpp
index f59d24f91786..f55d58f8b964 100644
--- a/src/Time/TimeSteppers/AdamsMoultonPc.hpp
+++ b/src/Time/TimeSteppers/AdamsMoultonPc.hpp
@@ -5,8 +5,11 @@
 
 #include <cstddef>
 #include <cstdint>
+#include <pup.h>
+#include <string>
 
 #include "Options/String.hpp"
+#include "Time/TimeSteppers/LtsTimeStepper.hpp"
 #include "Time/TimeSteppers/TimeStepper.hpp"
 #include "Utilities/Serialization/CharmPupable.hpp"
 #include "Utilities/TMPL.hpp"
@@ -19,7 +22,11 @@ class er;
 }  // namespace PUP
 namespace TimeSteppers {
 template <typename T>
+class BoundaryHistoryEvaluator;
+class ConstBoundaryHistoryTimes;
+template <typename T>
 class ConstUntypedHistory;
+class MutableBoundaryHistoryTimes;
 template <typename T>
 class MutableUntypedHistory;
 }  // namespace TimeSteppers
@@ -81,7 +88,7 @@ namespace TimeSteppers {
  * </table>
  */
 template <bool Monotonic>
-class AdamsMoultonPc : public TimeStepper {
+class AdamsMoultonPc : public LtsTimeStepper {
  public:
   static std::string name() {
     return Monotonic ? "AdamsMoultonPcMonotonic" : "AdamsMoultonPc";
@@ -131,6 +138,14 @@ class AdamsMoultonPc : public TimeStepper {
   TimeStepId next_time_id_for_error(const TimeStepId& current_id,
                                     const TimeDelta& time_step) const override;
 
+  bool neighbor_data_required(
+      const TimeStepId& next_substep_id,
+      const TimeStepId& neighbor_data_id) const override;
+
+  bool neighbor_data_required(
+      double dense_output_time,
+      const TimeStepId& neighbor_data_id) const override;
+
   WRAPPED_PUPable_decl_template(AdamsMoultonPc);  // NOLINT
 
   explicit AdamsMoultonPc(CkMigrateMessage* /*unused*/) {}
@@ -159,7 +174,29 @@ class AdamsMoultonPc : public TimeStepper {
   bool can_change_step_size_impl(const TimeStepId& time_id,
                                  const ConstUntypedHistory<T>& history) const;
 
+  template <typename T>
+  void add_boundary_delta_impl(
+      gsl::not_null<T*> result,
+      const TimeSteppers::ConstBoundaryHistoryTimes& local_times,
+      const TimeSteppers::ConstBoundaryHistoryTimes& remote_times,
+      const TimeSteppers::BoundaryHistoryEvaluator<T>& coupling,
+      const TimeDelta& time_step) const;
+
+  void clean_boundary_history_impl(
+      const TimeSteppers::MutableBoundaryHistoryTimes& local_times,
+      const TimeSteppers::MutableBoundaryHistoryTimes& remote_times)
+      const override;
+
+  template <typename T>
+  void boundary_dense_output_impl(
+      gsl::not_null<T*> result,
+      const TimeSteppers::ConstBoundaryHistoryTimes& local_times,
+      const TimeSteppers::ConstBoundaryHistoryTimes& remote_times,
+      const TimeSteppers::BoundaryHistoryEvaluator<T>& coupling,
+      double time) const;
+
   TIME_STEPPER_DECLARE_OVERLOADS
+  LTS_TIME_STEPPER_DECLARE_OVERLOADS
 
   size_t order_{};
 };
diff --git a/src/Time/TimeSteppers/CMakeLists.txt b/src/Time/TimeSteppers/CMakeLists.txt
index 4b6c9845a62f..8d609a2dd1d1 100644
--- a/src/Time/TimeSteppers/CMakeLists.txt
+++ b/src/Time/TimeSteppers/CMakeLists.txt
@@ -6,6 +6,7 @@ spectre_target_sources(
   PRIVATE
   AdamsBashforth.cpp
   AdamsCoefficients.cpp
+  AdamsLts.cpp
   AdamsMoultonPc.cpp
   ClassicalRungeKutta4.cpp
   DormandPrince5.cpp
@@ -28,6 +29,7 @@ spectre_target_headers(
   HEADERS
   AdamsBashforth.hpp
   AdamsCoefficients.hpp
+  AdamsLts.hpp
   AdamsMoultonPc.hpp
   ClassicalRungeKutta4.hpp
   DormandPrince5.hpp
diff --git a/src/Time/TimeSteppers/Factory.hpp b/src/Time/TimeSteppers/Factory.hpp
index dc0cee19bc57..1dc3babcd6f0 100644
--- a/src/Time/TimeSteppers/Factory.hpp
+++ b/src/Time/TimeSteppers/Factory.hpp
@@ -27,7 +27,8 @@ using time_steppers =
                Rk5Owren, Rk5Tsitouras>;
 
 /// Typelist of available LtsTimeSteppers
-using lts_time_steppers = tmpl::list<AdamsBashforth>;
+using lts_time_steppers =
+    tmpl::list<AdamsBashforth, AdamsMoultonPc<false>, AdamsMoultonPc<true>>;
 
 /// Typelist of available ImexTimeSteppers
 using imex_time_steppers =
@@ -36,4 +37,9 @@ using imex_time_steppers =
 /// Typelist of TimeSteppers whose substep times are strictly increasing
 using increasing_substep_time_steppers =
     tmpl::list<AdamsBashforth, Rk3Owren, Rk4Owren>;
+
+/// Typelist of LtsTimeSteppers with monotonic() true, i.e., those
+/// that work with control systems.
+using monotonic_lts_time_steppers =
+    tmpl::list<AdamsBashforth, AdamsMoultonPc<true>>;
 }  // namespace TimeSteppers
diff --git a/src/Time/TimeSteppers/LtsTimeStepper.hpp b/src/Time/TimeSteppers/LtsTimeStepper.hpp
index 00237c24b3e1..e778c930f80b 100644
--- a/src/Time/TimeSteppers/LtsTimeStepper.hpp
+++ b/src/Time/TimeSteppers/LtsTimeStepper.hpp
@@ -22,6 +22,8 @@ class TimeStepId;
 /// \cond
 #define LTS_TIME_STEPPER_WRAPPED_TYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
 #define LTS_TIME_STEPPER_DERIVED_CLASS(data) BOOST_PP_TUPLE_ELEM(1, data)
+#define LTS_TIME_STEPPER_DERIVED_CLASS_TEMPLATE(data) \
+  BOOST_PP_TUPLE_ELEM(2, data)
 /// \endcond
 
 /// \ingroup TimeSteppersGroup
@@ -73,7 +75,7 @@ class LtsTimeStepper : public virtual TimeStepper {
       const TimeSteppers::ConstBoundaryHistoryTimes& remote_times, \
       const TimeSteppers::BoundaryHistoryEvaluator<                \
           LTS_TIME_STEPPER_WRAPPED_TYPE(data)>& coupling,          \
-      const double time) const = 0;
+      double time) const = 0;
 
   GENERATE_INSTANTIATIONS(LTS_TIME_STEPPER_DECLARE_VIRTUALS_IMPL,
                           (MATH_WRAPPER_TYPES))
@@ -162,7 +164,7 @@ class LtsTimeStepper : public virtual TimeStepper {
   ///     const TimeSteppers::ConstBoundaryHistoryTimes& local_times,
   ///     const TimeSteppers::ConstBoundaryHistoryTimes& remote_times,
   ///     const TimeSteppers::BoundaryHistoryEvaluator<T>& coupling,
-  ///     const double time) const;
+  ///     double time) const;
   /// ```
   template <typename LocalVars, typename RemoteVars, typename Coupling>
   void boundary_dense_output(
@@ -196,9 +198,10 @@ class LtsTimeStepper : public virtual TimeStepper {
       const TimeSteppers::ConstBoundaryHistoryTimes& remote_times, \
       const TimeSteppers::BoundaryHistoryEvaluator<                \
           LTS_TIME_STEPPER_WRAPPED_TYPE(data)>& coupling,          \
-      const double time) const override;
+      double time) const override;
 
 #define LTS_TIME_STEPPER_DEFINE_OVERLOADS_IMPL(_, data)                     \
+  LTS_TIME_STEPPER_DERIVED_CLASS_TEMPLATE(data)                             \
   void LTS_TIME_STEPPER_DERIVED_CLASS(data)::add_boundary_delta_forward(    \
       const gsl::not_null<LTS_TIME_STEPPER_WRAPPED_TYPE(data)*> result,     \
       const TimeSteppers::ConstBoundaryHistoryTimes& local_times,           \
@@ -209,6 +212,7 @@ class LtsTimeStepper : public virtual TimeStepper {
     return add_boundary_delta_impl(result, local_times, remote_times,       \
                                    coupling, time_step);                    \
   }                                                                         \
+  LTS_TIME_STEPPER_DERIVED_CLASS_TEMPLATE(data)                             \
   void LTS_TIME_STEPPER_DERIVED_CLASS(data)::boundary_dense_output_forward( \
       const gsl::not_null<LTS_TIME_STEPPER_WRAPPED_TYPE(data)*> result,     \
       const TimeSteppers::ConstBoundaryHistoryTimes& local_times,           \
@@ -233,6 +237,13 @@ class LtsTimeStepper : public virtual TimeStepper {
 /// Macro defining overloaded detail methods in classes derived from
 /// TimeStepper.  Must be placed in the cpp file for the derived
 /// class.
+/// @{
 #define LTS_TIME_STEPPER_DEFINE_OVERLOADS(derived_class)          \
   GENERATE_INSTANTIATIONS(LTS_TIME_STEPPER_DEFINE_OVERLOADS_IMPL, \
-                          (MATH_WRAPPER_TYPES), (derived_class))
+                          (MATH_WRAPPER_TYPES), (derived_class), ())
+#define LTS_TIME_STEPPER_DEFINE_OVERLOADS_TEMPLATED(derived_class, \
+                                                    template_args) \
+  GENERATE_INSTANTIATIONS(LTS_TIME_STEPPER_DEFINE_OVERLOADS_IMPL,  \
+                          (MATH_WRAPPER_TYPES), (derived_class),   \
+                          (template <template_args>))
+/// @}
diff --git a/src/Time/TimeSteppers/TimeStepper.hpp b/src/Time/TimeSteppers/TimeStepper.hpp
index 367ac4909c8a..ba5da1cf2f79 100644
--- a/src/Time/TimeSteppers/TimeStepper.hpp
+++ b/src/Time/TimeSteppers/TimeStepper.hpp
@@ -69,7 +69,7 @@ class TimeStepper : public PUP::able {
       const gsl::not_null<TIME_STEPPER_WRAPPED_TYPE(data)*> u,             \
       const TimeSteppers::ConstUntypedHistory<TIME_STEPPER_WRAPPED_TYPE(   \
           data)>& history,                                                 \
-      const double time) const = 0;                                        \
+      double time) const = 0;                                              \
   virtual bool can_change_step_size_forward(                               \
       const TimeStepId& time_id,                                           \
       const TimeSteppers::ConstUntypedHistory<TIME_STEPPER_WRAPPED_TYPE(   \
diff --git a/tests/Unit/Helpers/Time/TimeSteppers/LtsHelpers.cpp b/tests/Unit/Helpers/Time/TimeSteppers/LtsHelpers.cpp
index d9fadeb6907b..6d6d46a84dfc 100644
--- a/tests/Unit/Helpers/Time/TimeSteppers/LtsHelpers.cpp
+++ b/tests/Unit/Helpers/Time/TimeSteppers/LtsHelpers.cpp
@@ -6,169 +6,499 @@
 #include "Helpers/Time/TimeSteppers/LtsHelpers.hpp"
 
 #include <algorithm>
-#include <array>
 #include <cmath>
 #include <cstddef>
 #include <cstdint>
 #include <deque>
+#include <limits>
+#include <optional>
+#include <utility>
 
+#include "Helpers/Time/TimeSteppers/TimeStepperTestUtils.hpp"
 #include "Time/BoundaryHistory.hpp"
+#include "Time/EvolutionOrdering.hpp"
 #include "Time/History.hpp"
 #include "Time/Slab.hpp"
 #include "Time/Time.hpp"
 #include "Time/TimeStepId.hpp"
 #include "Time/TimeSteppers/LtsTimeStepper.hpp"
+#include "Utilities/Algorithm.hpp"
 #include "Utilities/ErrorHandling/Assert.hpp"
 #include "Utilities/Gsl.hpp"
+#include "Utilities/Rational.hpp"
 
 namespace TimeStepperTestUtils::lts {
-void test_equal_rate(const LtsTimeStepper& stepper, const size_t order,
-                     const size_t number_of_past_steps, const double epsilon,
-                     const bool forward) {
-  // This does an integral putting the entire derivative into the
-  // boundary term.
-  const double unused_local_deriv = 4444.;
-
-  auto analytic = [](double t) { return sin(t); };
-  auto driver = [](double t) { return cos(t); };
-  auto coupling = [=](const double local, const double remote) {
-    CHECK(local == unused_local_deriv);
-    return remote;
-  };
-
-  Approx approx = Approx::custom().epsilon(epsilon);
+namespace {
+void arbitrary_compare_to_volume(const LtsTimeStepper& stepper,
+                                 const bool time_runs_forward,
+                                 std::deque<Rational> step_pattern,
+                                 std::deque<Rational> neighbor_step_pattern,
+                                 const double tolerance) {
+  const auto order = stepper.order();
+  const auto local_approx = Approx::custom().epsilon(tolerance);
+  const auto step_sign = time_runs_forward ? 1 : -1;
 
-  const uint64_t num_steps = 100;
-  const Slab slab(0.875, 1.);
-  const TimeDelta step_size = (forward ? 1 : -1) * slab.duration() / num_steps;
+  // Arbitrary function that gives "generic" results.  We don't
+  // attempt to perform a meaningful integral.
+  const auto arbitrary_rhs = [](const TimeStepId& t) {
+    return 1.23 * static_cast<double>(t.substep()) + t.substep_time();
+  };
 
-  TimeStepId time_id(forward, 0, forward ? slab.start() : slab.end());
-  double y = analytic(time_id.substep_time());
   TimeSteppers::History<double> volume_history{order};
   TimeSteppers::BoundaryHistory<double, double, double> boundary_history{};
 
+  const Slab initial_slab(1.2, 3.4);
+  TimeStepId id(time_runs_forward, 0,
+                time_runs_forward ? initial_slab.start() : initial_slab.end());
   {
-    Time history_time = time_id.step_time();
-    TimeDelta history_step_size = step_size;
-    for (size_t j = 0; j < number_of_past_steps; ++j) {
-      ASSERT(history_time.slab() == history_step_size.slab(), "Slab mismatch");
-      if ((history_step_size.is_positive() and
-           history_time.is_at_slab_start()) or
-          (not history_step_size.is_positive() and
-           history_time.is_at_slab_end())) {
-        const Slab new_slab =
-            history_time.slab().advance_towards(-history_step_size);
-        history_time = history_time.with_slab(new_slab);
-        history_step_size = history_step_size.with_slab(new_slab);
-      }
-      history_time -= history_step_size;
-      const TimeStepId history_id(forward, 0, history_time);
-      volume_history.insert_initial(history_id, analytic(history_time.value()),
-                                    0.);
-      boundary_history.local().insert_initial(history_id, order,
-                                              unused_local_deriv);
-      boundary_history.remote().insert_initial(history_id, order,
-                                               driver(history_time.value()));
+    auto init_step = step_sign * initial_slab.duration();
+    auto init_time = id.step_time();
+    int64_t slab_number = 0;
+    while (volume_history.size() < stepper.number_of_past_steps()) {
+      --slab_number;
+      init_step =
+          init_step.with_slab(init_step.slab().advance_towards(-init_step));
+      init_time -= init_step;
+      const TimeStepId init_id(time_runs_forward, slab_number, init_time);
+      const double deriv = arbitrary_rhs(init_id);
+      volume_history.insert_initial(init_id, 0.0, deriv);
+      boundary_history.local().insert_initial(init_id, order, deriv);
+      // Currently unset and unused in the evolution code.
+      constexpr auto remote_order = std::numeric_limits<size_t>::max();
+      boundary_history.remote().insert_initial(init_id, remote_order, deriv);
     }
   }
 
-  for (uint64_t i = 0; i < num_steps; ++i) {
-    for (uint64_t substep = 0;
-         substep < stepper.number_of_substeps();
-         ++substep) {
-      volume_history.insert(time_id, y, 0.);
-      boundary_history.local().insert(time_id, order, unused_local_deriv);
-      boundary_history.remote().insert(time_id, order,
-                                       driver(time_id.substep_time()));
+  const bool equal_rate = step_pattern == neighbor_step_pattern;
+  const auto coupling = [&](const double local, const double remote) {
+    if (equal_rate) {
+      CHECK(local == remote);
+    }
+    return local;
+  };
+
+  auto neighbor_id = id;
+
+  const auto add_neighbor_entries = [&](const auto& limit) {
+    while (stepper.neighbor_data_required(limit, neighbor_id)) {
+      boundary_history.remote().insert(neighbor_id, order,
+                                       arbitrary_rhs(neighbor_id));
+      ASSERT(not neighbor_step_pattern.empty(), "Test logic error");
+      const TimeDelta neighbor_step = step_sign *
+                                      neighbor_step_pattern.front() *
+                                      neighbor_id.step_time().slab().duration();
+      neighbor_id = stepper.next_time_id(neighbor_id, neighbor_step);
+      if (neighbor_id.substep() == 0) {
+        neighbor_step_pattern.pop_front();
+      }
+    }
+  };
+
+  while (not step_pattern.empty()) {
+    const TimeDelta step =
+        step_sign * step_pattern.front() * id.step_time().slab().duration();
+    step_pattern.pop_front();
+
+    const double dense_test = (id.step_time() + step / 3).value();
+    bool dense_succeeded = false;
+    do {
+      {
+        const double deriv = arbitrary_rhs(id);
+        volume_history.insert(id, 0.0, deriv);
+        boundary_history.local().insert(id, order, deriv);
+      }
+
+      // Check dense output
+      {
+        double volume_value = 0.0;
+        if (stepper.dense_update_u(make_not_null(&volume_value), volume_history,
+                                   dense_test)) {
+          CHECK(not dense_succeeded);
+          dense_succeeded = true;
+          add_neighbor_entries(dense_test);
+          double boundary_value = 0.0;
+          stepper.boundary_dense_output(make_not_null(&boundary_value),
+                                        boundary_history, dense_test, coupling);
+          CHECK(volume_value == local_approx(boundary_value));
+        }
+      }
+
+      id = stepper.next_time_id(id, step);
+      add_neighbor_entries(id);
 
-      stepper.update_u(make_not_null(&y), volume_history, step_size);
+      // Check normal output
+      {
+        double volume_value = 0.0;
+        stepper.update_u(make_not_null(&volume_value), volume_history, step);
+        double boundary_value = 0.0;
+        stepper.add_boundary_delta(make_not_null(&boundary_value),
+                                   boundary_history, step, coupling);
+        CHECK(volume_value == local_approx(boundary_value));
+      }
       stepper.clean_history(make_not_null(&volume_history));
-      stepper.add_boundary_delta(&y, boundary_history, step_size, coupling);
       stepper.clean_boundary_history(make_not_null(&boundary_history));
-      time_id = stepper.next_time_id(time_id, step_size);
-    }
-    CHECK(y == approx(analytic(time_id.substep_time())));
+    } while (id.substep() != 0);
+    CHECK(dense_succeeded);
   }
-  // Make sure history is being cleaned up.  The limit of 20 is
-  // arbitrary, but much larger than the order of any integrators we
-  // care about and much smaller than the number of time steps in the
-  // test.
-  CHECK(boundary_history.local().size() < 20);
-  CHECK(boundary_history.remote().size() < 20);
 }
+}  // namespace
 
-void test_dense_output(const LtsTimeStepper& stepper) {
-  // We only support variable time-step, multistep LTS integration.
-  // Any multistep, variable time-step integrator must give the same
-  // results from dense output as from just taking a short step
-  // because we require dense output to be continuous.  A sufficient
-  // test is therefore to run with an LTS pattern and check that the
-  // dense output predicts the actual step result.
-  const Slab slab(0., 1.);
+void test_equal_rate(const LtsTimeStepper& stepper) {
+  INFO("test_equal_rate");
+  const std::deque<Rational> step_pattern{{1}, {1, 4}, {1, 4}, {1, 2}};
+  arbitrary_compare_to_volume(stepper, true, step_pattern, step_pattern,
+                              1.0e-15);
+  arbitrary_compare_to_volume(stepper, false, step_pattern, step_pattern,
+                              1.0e-15);
+}
 
-  // We don't use any meaningful values.  We only care that the dense
-  // output gives the same result as normal output.
-  // NOLINTNEXTLINE(spectre-mutable)
-  auto get_value = [value = 1.]() mutable { return value *= 1.1; };
+void test_uncoupled(const LtsTimeStepper& stepper, const double tolerance) {
+  INFO("test_uncoupled");
+  const std::deque<Rational> step_pattern{{1}, {1, 4}, {1, 4}, {1, 2}};
+  const std::deque<Rational> neighbor_step_pattern{{1, 4}, {1, 4}, {1, 2}, {1}};
+  arbitrary_compare_to_volume(stepper, true, step_pattern,
+                              neighbor_step_pattern, tolerance);
+  arbitrary_compare_to_volume(stepper, false, step_pattern,
+                              neighbor_step_pattern, tolerance);
+}
 
-  const auto coupling = [](const double a, const double b) { return a * b; };
+namespace {
+template <typename Coupling>
+class Element {
+ public:
+  template <typename AnalyticSelf, typename AnalyticNeighbor>
+  Element(const LtsTimeStepper& stepper, Coupling coupling,
+          std::deque<Rational> step_pattern, const TimeStepId& time_step_id,
+          const Rational& neighbor_first_step,
+          const AnalyticSelf& analytic_self,
+          const AnalyticNeighbor& analytic_neighbor)
+      : stepper_(&stepper),
+        coupling_(std::move(coupling)),
+        step_pattern_(std::move(step_pattern)),
+        time_step_id_(time_step_id),
+        value_(analytic_self(time_step_id_.substep_time())),
+        step_size_(next_step_size()),
+        next_time_step_id_(stepper_->next_time_id(time_step_id_, step_size_)),
+        volume_history_(stepper_->order()),
+        next_message_(stepper_->next_time_id(
+            time_step_id_, neighbor_first_step *
+                               time_step_id_.step_time().slab().duration())) {
+    auto init_step = (time_step_id_.time_runs_forward() ? 1 : -1) *
+                     time_step_id_.step_time().slab().duration();
+    auto init_time = time_step_id_.step_time();
+    int64_t slab_number = 0;
+    while (volume_history_.size() < stepper_->number_of_past_steps() + 1) {
+      const TimeStepId init_id(time_step_id_.time_runs_forward(), slab_number,
+                               init_time);
+      volume_history_.insert_initial(init_id, analytic_self(init_time.value()),
+                                     0.0);
+      boundary_history_.local().insert_initial(
+          init_id, stepper_->order(), analytic_self(init_time.value()));
+      // Currently unset and unused in the evolution code.
+      constexpr auto remote_order = std::numeric_limits<size_t>::max();
+      boundary_history_.remote().insert_initial(
+          init_id, remote_order, analytic_neighbor(init_time.value()));
+      --slab_number;
+      init_step =
+          init_step.with_slab(init_step.slab().advance_towards(-init_step));
+      init_time -= init_step;
+    }
+  }
 
-  const auto make_time_id = [](const Time& t) {
-    return TimeStepId(true, 0, t);
-  };
+  bool done() const { return step_pattern_.empty(); }
+  TimeStepId next_time_step_id() const { return next_time_step_id_; }
+  Time next_step_time() const { return time_step_id_.step_time() + step_size_; }
 
-  TimeSteppers::BoundaryHistory<double, double, double> history{};
-  {
-    const Slab init_slab = slab.retreat();
-    for (size_t i = 0; i < stepper.number_of_past_steps(); ++i) {
-      const Time init_time =
-          init_slab.end() -
-          init_slab.duration() * (i + 1) / stepper.number_of_past_steps();
-      history.local().insert_initial(make_time_id(init_time), stepper.order(),
-                                     get_value());
-      history.remote().insert_initial(make_time_id(init_time), stepper.order(),
-                                      get_value());
+  bool local_dense_output_ready(const double time) const {
+    const evolution_less<double> before{time_step_id_.time_runs_forward()};
+    if (not before(time, next_step_time().value())) {
+      return false;
+    }
+    double dummy_value = 0.0;
+    return stepper_->dense_update_u(make_not_null(&dummy_value),
+                                    volume_history_, time);
+  }
+
+  double dense_output(const double time) {
+    REQUIRE(process_messages(time));
+    double dense_result = *volume_history_.complete_step_start().value;
+    // Can skip the volume update because all derivatives are zero.
+    stepper_->boundary_dense_output(make_not_null(&dense_result),
+                                    boundary_history_, time, coupling_);
+    return dense_result;
+  }
+
+  std::optional<std::pair<TimeStepId, double>> step() {
+    if (not process_messages(next_time_step_id_)) {
+      return std::nullopt;
+    }
+    stepper_->update_u(make_not_null(&value_), volume_history_, step_size_);
+    stepper_->add_boundary_delta(make_not_null(&value_), boundary_history_,
+                                 step_size_, coupling_);
+    stepper_->clean_history(make_not_null(&volume_history_));
+    stepper_->clean_boundary_history(make_not_null(&boundary_history_));
+
+    time_step_id_ = next_time_step_id_;
+    if (time_step_id_.substep() == 0) {
+      step_pattern_.pop_front();
+      if (not done()) {
+        step_size_ = next_step_size();
+      }
+      step_size_ = step_size_.with_slab(time_step_id_.step_time().slab());
+    }
+    next_time_step_id_ = stepper_->next_time_id(time_step_id_, step_size_);
+    volume_history_.insert(time_step_id_, value_, 0.0);
+    boundary_history_.local().insert(
+        time_step_id_, volume_history_.integration_order(), value_);
+    return {{time_step_id_, value_}};
+  }
+
+  void receive_data(const TimeStepId& id, const double value,
+                    const TimeStepId& next_id) {
+    ASSERT(id == next_message_, "Test logic error.");
+    messages_.emplace_back(id, value);
+    next_message_ = next_id;
+  }
+
+ private:
+  TimeDelta next_step_size() const {
+    return step_pattern_.front() * time_step_id_.step_time().slab().duration();
+  }
+
+  template <typename T>
+  bool process_messages(const T& time) {
+    // Currently unset and unused in the evolution code.
+    constexpr auto integration_order = std::numeric_limits<size_t>::max();
+    while (not messages_.empty() and
+           stepper_->neighbor_data_required(time, messages_.front().first)) {
+      boundary_history_.remote().insert(
+          messages_.front().first, integration_order, messages_.front().second);
+      messages_.pop_front();
     }
+    return not(messages_.empty() and
+               stepper_->neighbor_data_required(time, next_message_));
   }
 
-  std::array<std::deque<TimeDelta>, 2> dt{
-      {{slab.duration() / 2, slab.duration() / 4, slab.duration() / 4},
-       {slab.duration() / 6, slab.duration() / 6, slab.duration() * 2 / 9,
-        slab.duration() * 4 / 9}}};
+  gsl::not_null<const LtsTimeStepper*> stepper_;
+  Coupling coupling_;
+
+  std::deque<Rational> step_pattern_;
+  TimeStepId time_step_id_;
+  double value_;
+  TimeDelta step_size_;
+  TimeStepId next_time_step_id_;
+  TimeSteppers::History<double> volume_history_{};
+  TimeSteppers::BoundaryHistory<double, double, double> boundary_history_{};
+  std::deque<std::pair<TimeStepId, double>> messages_{};
+  TimeStepId next_message_;
+};
+
+// Test system:
+// dx/dt = x y
+// dy/dt = - x y
+//
+// Solution:
+// x = c / [1 + exp(-c (t - d))]
+// y = c - x
+namespace product_system {
+double rhs_x(const double x, const double y) { return x * y; }
+double rhs_y(const double x, const double y) { return -x * y; }
+
+// Arbitrary
+constexpr double conserved_sum = 0.7;  // c
+constexpr double offset = 0.4;  // d
+
+double analytic_x(const double t) {
+  return conserved_sum / (1.0 + exp(-conserved_sum * (t - offset)));
+}
+double analytic_y(const double t) { return conserved_sum - analytic_x(t); }
+}  // namespace product_system
+
+void test_conservation_impl(const LtsTimeStepper& stepper,
+                            std::deque<Rational> step_pattern_x,
+                            std::deque<Rational> step_pattern_y) {
+  const bool time_runs_forward = step_pattern_x.front() > 0;
+  const evolution_less<Time> before{time_runs_forward};
+
+  const Slab initial_slab(1.2, 3.4);
+  const TimeStepId initial_id(
+      time_runs_forward, 0,
+      time_runs_forward ? initial_slab.start() : initial_slab.end());
+
+  const auto first_step_x = step_pattern_x.front();
+  Element element_x(
+      stepper,
+      [](const double local, const double remote) {
+        return product_system::rhs_x(local, remote);
+      },
+      std::move(step_pattern_x), initial_id, step_pattern_y.front(),
+      product_system::analytic_x, product_system::analytic_y);
+  Element element_y(
+      stepper,
+      [](const double local, const double remote) {
+        return product_system::rhs_y(remote, local);
+      },
+      std::move(step_pattern_y), initial_id, first_step_x,
+      product_system::analytic_y, product_system::analytic_x);
+
+  double test_time = initial_id.step_time().value();
 
-  Time t = slab.start();
-  Time next_check = t + dt[0][0];
-  std::array<Time, 2> next{{t, t}};
   for (;;) {
-    const size_t side = next[0] <= next[1] ? 0 : 1;
-
-    if (side == 0) {
-      history.local().insert(make_time_id(next[0]), stepper.order(),
-                             get_value());
-    } else {
-      history.remote().insert(make_time_id(next[1]), stepper.order(),
-                              get_value());
+    const bool need_more_x = not element_x.local_dense_output_ready(test_time);
+    const bool need_more_y = not element_y.local_dense_output_ready(test_time);
+    // If an element can produce dense output, it must have
+    // transmitted enough data for its neighbor to do the same.
+    if (not(need_more_x or need_more_y)) {
+      const double dense_x = element_x.dense_output(test_time);
+      const double dense_y = element_y.dense_output(test_time);
+      CHECK(dense_x + dense_y == approx(product_system::conserved_sum));
+      test_time = std::min(element_x.next_step_time(),
+                           element_y.next_step_time(), before)
+                      .value();
+      continue;
     }
 
-    const TimeDelta this_dt = gsl::at(dt, side).front();
-    gsl::at(dt, side).pop_front();
+    if (element_x.done() and element_y.done()) {
+      break;
+    }
+
+    if (need_more_x) {
+      const auto step_result = element_x.step();
+      if (step_result.has_value()) {
+        element_y.receive_data(step_result->first, step_result->second,
+                               element_x.next_time_step_id());
+        continue;
+      }
+    }
+
+    if (need_more_y) {
+      const auto step_result = element_y.step();
+      if (step_result.has_value()) {
+        element_x.receive_data(step_result->first, step_result->second,
+                               element_y.next_time_step_id());
+        continue;
+      }
+    }
+
+    INFO("Deadlocked");
+    REQUIRE(false);
+  }
+}
 
-    gsl::at(next, side) += this_dt;
+double test_convergence_error(const LtsTimeStepper& stepper,
+                              const std::deque<Rational>& step_pattern_x,
+                              const std::deque<Rational>& step_pattern_y,
+                              const int32_t repeats, const bool test_dense) {
+  const bool time_runs_forward = step_pattern_x.front() > 0;
 
-    if (std::min(next[0], next[1]) == next_check) {
-      double dense_result = 0.0;
-      stepper.boundary_dense_output(&dense_result, history, next_check.value(),
-                                    coupling);
-      double delta = 0.0;
-      stepper.add_boundary_delta(&delta, history, next_check - t, coupling);
-      stepper.clean_boundary_history(make_not_null(&history));
-      CHECK(dense_result == approx(delta));
-      if (next_check.is_at_slab_boundary()) {
+  const auto initial_slab =
+      time_runs_forward ? Slab::with_duration_from_start(-10.0, 5.0 / repeats)
+                        : Slab::with_duration_to_end(-5.0, 5.0 / repeats);
+  // Arbitrary but not a nice rational portion of a slab.
+  const double dense_time = -7.91237;
+
+  const TimeStepId initial_id(
+      time_runs_forward, 0,
+      time_runs_forward ? initial_slab.start() : initial_slab.end());
+
+  std::deque<Rational> full_pattern_x{};
+  std::deque<Rational> full_pattern_y{};
+  for (int32_t i = 0; i < repeats; ++i) {
+    full_pattern_x.insert(full_pattern_x.end(), step_pattern_x.begin(),
+                          step_pattern_x.end());
+    full_pattern_y.insert(full_pattern_y.end(), step_pattern_y.begin(),
+                          step_pattern_y.end());
+  }
+
+  Element element_x(
+      stepper,
+      [](const double local, const double remote) {
+        return product_system::rhs_x(local, remote);
+      },
+      std::move(full_pattern_x), initial_id, step_pattern_y.front(),
+      product_system::analytic_x, product_system::analytic_y);
+  Element element_y(
+      stepper,
+      [](const double local, const double remote) {
+        return product_system::rhs_y(remote, local);
+      },
+      std::move(full_pattern_y), initial_id, step_pattern_x.front(),
+      product_system::analytic_y, product_system::analytic_x);
+
+  for (;;) {
+    // Step y as much as possible so we don't have to worry about
+    // having enough data on the x side for dense output.
+    while (not element_y.done()) {
+      const auto step_result = element_y.step();
+      if (not step_result.has_value()) {
         break;
       }
-      t = next_check;
-      next_check += dt[0].front();
+      element_x.receive_data(step_result->first, step_result->second,
+                             element_y.next_time_step_id());
+    }
+
+    if (test_dense and element_x.local_dense_output_ready(dense_time)) {
+      return element_x.dense_output(dense_time) -
+             product_system::analytic_x(dense_time);
+    }
+
+    {
+      REQUIRE(not element_x.done());
+      const auto step_result = element_x.step();
+      if (step_result.has_value()) {
+        if (not test_dense and element_x.done()) {
+          return step_result->second -
+                 product_system::analytic_x(step_result->first.substep_time());
+        }
+        element_y.receive_data(step_result->first, step_result->second,
+                               element_x.next_time_step_id());
+      }
     }
   }
 }
+
+void test_convergence_impl(
+    const LtsTimeStepper& stepper,
+    const std::pair<int32_t, int32_t>& number_of_steps_range,
+    const int32_t stride, const bool dense) {
+  CAPTURE(dense);
+  std::deque<Rational> step_pattern_x{{1, 2}, {1, 8}, {1, 8}, {1, 4}};
+  std::deque<Rational> step_pattern_y{{1, 8}, {1, 8}, {1, 4}, {1, 2}};
+  CHECK(convergence_rate(
+            number_of_steps_range, stride, [&](const int32_t repeats) {
+              return test_convergence_error(stepper, step_pattern_x,
+                                            step_pattern_y, repeats, dense);
+            }) == approx(stepper.order()).margin(0.4));
+  alg::for_each(step_pattern_x, [](Rational& x) { x *= -1; });
+  alg::for_each(step_pattern_y, [](Rational& x) { x *= -1; });
+  CHECK(convergence_rate(
+            number_of_steps_range, stride, [&](const int32_t repeats) {
+              return test_convergence_error(stepper, step_pattern_x,
+                                            step_pattern_y, repeats, dense);
+            }) == approx(stepper.order()).margin(0.4));
+}
+}  // namespace
+
+void test_conservation(const LtsTimeStepper& stepper) {
+  std::deque<Rational> step_pattern_x{{1}, {1, 4}, {1, 4}, {1, 2}};
+  std::deque<Rational> step_pattern_y{{1, 4}, {1, 4}, {1, 2}, {1}};
+  test_conservation_impl(stepper, step_pattern_x, step_pattern_y);
+  alg::for_each(step_pattern_x, [](Rational& x) { x *= -1; });
+  alg::for_each(step_pattern_y, [](Rational& x) { x *= -1; });
+  test_conservation_impl(stepper, std::move(step_pattern_x),
+                         std::move(step_pattern_y));
+}
+
+void test_convergence(const LtsTimeStepper& stepper,
+                      const std::pair<int32_t, int32_t>& number_of_steps_range,
+                      const int32_t stride) {
+  test_convergence_impl(stepper, number_of_steps_range, stride, false);
+}
+
+void test_dense_convergence(
+    const LtsTimeStepper& stepper,
+    const std::pair<int32_t, int32_t>& number_of_steps_range,
+    const int32_t stride) {
+  test_convergence_impl(stepper, number_of_steps_range, stride, true);
+}
 }  // namespace TimeStepperTestUtils::lts
diff --git a/tests/Unit/Helpers/Time/TimeSteppers/LtsHelpers.hpp b/tests/Unit/Helpers/Time/TimeSteppers/LtsHelpers.hpp
index eed00c5aa4f2..cdb554b6fca5 100644
--- a/tests/Unit/Helpers/Time/TimeSteppers/LtsHelpers.hpp
+++ b/tests/Unit/Helpers/Time/TimeSteppers/LtsHelpers.hpp
@@ -4,6 +4,8 @@
 #pragma once
 
 #include <cstddef>
+#include <cstdint>
+#include <utility>
 
 /// \cond
 class LtsTimeStepper;
@@ -11,11 +13,23 @@ class LtsTimeStepper;
 
 namespace TimeStepperTestUtils::lts {
 /// Test boundary computations with the same step size on both
-/// neighbors.
-void test_equal_rate(const LtsTimeStepper& stepper, size_t order,
-                     size_t number_of_past_steps, double epsilon, bool forward);
+/// neighbors against the volume computation.
+void test_equal_rate(const LtsTimeStepper& stepper);
 
+/// Test uncoupled boundary computations against the volume
+/// computation.
+void test_uncoupled(const LtsTimeStepper& stepper, double tolerance);
 
-/// Test the accuracy of boundary dense output.
-void test_dense_output(const LtsTimeStepper& stepper);
+/// Test conservation of boundary dense output.
+void test_conservation(const LtsTimeStepper& stepper);
+
+// Test convergence rate of boundary integration.
+void test_convergence(const LtsTimeStepper& stepper,
+                      const std::pair<int32_t, int32_t>& number_of_steps_range,
+                      int32_t stride);
+
+// Test convergence rate of boundary dense output.
+void test_dense_convergence(
+    const LtsTimeStepper& stepper,
+    const std::pair<int32_t, int32_t>& number_of_steps_range, int32_t stride);
 }  // namespace TimeStepperTestUtils::lts
diff --git a/tests/Unit/Time/TimeSteppers/CMakeLists.txt b/tests/Unit/Time/TimeSteppers/CMakeLists.txt
index ed755c90932c..81e0e97886fa 100644
--- a/tests/Unit/Time/TimeSteppers/CMakeLists.txt
+++ b/tests/Unit/Time/TimeSteppers/CMakeLists.txt
@@ -5,6 +5,7 @@ set(LIBRARY_SOURCES
   ${LIBRARY_SOURCES}
   TimeSteppers/Test_AdamsBashforth.cpp
   TimeSteppers/Test_AdamsCoefficients.cpp
+  TimeSteppers/Test_AdamsLts.cpp
   TimeSteppers/Test_AdamsMoultonPc.cpp
   TimeSteppers/Test_ClassicalRungeKutta4.cpp
   TimeSteppers/Test_DormandPrince5.cpp
diff --git a/tests/Unit/Time/TimeSteppers/Test_AdamsBashforth.cpp b/tests/Unit/Time/TimeSteppers/Test_AdamsBashforth.cpp
index 1e945677b528..1f71b7b6b526 100644
--- a/tests/Unit/Time/TimeSteppers/Test_AdamsBashforth.cpp
+++ b/tests/Unit/Time/TimeSteppers/Test_AdamsBashforth.cpp
@@ -4,14 +4,10 @@
 #include "Framework/TestingFramework.hpp"
 
 #include <algorithm>
-#include <array>
 #include <cmath>
 #include <cstddef>
 #include <cstdint>
-#include <deque>
-#include <initializer_list>
 
-#include "DataStructures/MathWrapper.hpp"
 #include "Framework/TestCreation.hpp"
 #include "Framework/TestHelpers.hpp"
 #include "Helpers/Time/TimeSteppers/LtsHelpers.hpp"
@@ -25,10 +21,8 @@
 #include "Time/TimeSteppers/LtsTimeStepper.hpp"
 #include "Time/TimeSteppers/TimeStepper.hpp"
 #include "Utilities/ErrorHandling/Assert.hpp"
-#include "Utilities/ForceInline.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/Literals.hpp"
-#include "Utilities/MakeWithValue.hpp"
 
 SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.AdamsBashforth", "[Unit][Time]") {
   for (size_t order = 1; order < 9; ++order) {
@@ -153,193 +147,6 @@ SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.AdamsBashforth.Backwards",
 }
 
 namespace {
-// Non-copyable double to verify that the boundary code is not making
-// internal copies.
-class NCd {
- public:
-  NCd() = default;
-  explicit NCd(double x) : x_(x) {}
-  NCd(const NCd&) = delete;
-  NCd(NCd&&) = default;
-  NCd& operator=(const NCd&) = delete;
-  NCd& operator=(NCd&&) = default;
-  ~NCd() = default;
-
-  const double& operator()() const { return x_; }
-  double& operator()() { return x_; }
-
- private:
-  double x_;
-};
-
-auto make_math_wrapper(const gsl::not_null<NCd*> x) {
-  return ::make_math_wrapper(&(*x)());
-}
-
-auto make_math_wrapper(const NCd& x) { return ::make_math_wrapper(x()); }
-
-// Random numbers
-constexpr double c10 = 0.949716728952811;
-constexpr double c11 = 0.190663110072823;
-constexpr double c20 = 0.932407227651314;
-constexpr double c21 = 0.805454101952822;
-constexpr double c22 = 0.825876851406978;
-
-// Test coupling for integrating using two drivers (multiplied together)
-NCd quartic_coupling(const NCd& local, const NCd& remote) {
-  return NCd(local() * remote());
-}
-
-// Test functions for integrating a quartic using the above coupling.
-// The derivative of quartic_answer is the product of the other two.
-double quartic_side1(double x) { return c10 + x * c11; }
-double quartic_side2(double x) { return c20 + x * (c21 + x * c22); }
-double quartic_answer(double x) {
-  return x * (c10 * c20
-              + x * ((c10 * c21 + c11 * c20) / 2
-                     + x * ((c10 * c22 + c11 * c21) / 3
-                            + x * (c11 * c22 / 4))));
-}
-
-void do_lts_test(const std::array<TimeDelta, 2>& dt) {
-  // For general time steppers the boundary stepper cannot be run
-  // without simultaneously running the volume stepper.  For
-  // Adams-Bashforth methods, however, the volume contribution is zero
-  // if all the derivative contributions are from the boundary, so we
-  // can leave it out.
-
-  const bool forward_in_time = dt[0].is_positive();
-  const auto simulation_less = [forward_in_time](const Time& a, const Time& b) {
-    return forward_in_time ? a < b : b < a;
-  };
-
-  const auto make_time_id = [forward_in_time](const Time& t) {
-    return TimeStepId(forward_in_time, 0, t);
-  };
-
-  const Slab slab = dt[0].slab();
-  Time t = forward_in_time ? slab.start() : slab.end();
-
-  const size_t order = 4;
-  TimeSteppers::AdamsBashforth ab4(order);
-
-  TimeSteppers::BoundaryHistory<NCd, NCd, NCd> history{};
-  {
-    const Slab init_slab = slab.advance_towards(-dt[0]);
-
-    for (int32_t step = 1; step <= 3; ++step) {
-      {
-        const Time now = t - step * dt[0].with_slab(init_slab);
-        history.local().insert_initial(make_time_id(now), order,
-                                       NCd(quartic_side1(now.value())));
-      }
-      {
-        const Time now = t - step * dt[1].with_slab(init_slab);
-        history.remote().insert_initial(make_time_id(now), order,
-                                        NCd(quartic_side2(now.value())));
-      }
-    }
-  }
-
-  NCd y(quartic_answer(t.value()));
-  Time next_check = t + dt[0];
-  std::array<Time, 2> next{{t, t}};
-  for (;;) {
-    const auto side = static_cast<size_t>(
-        std::min_element(next.cbegin(), next.cend(), simulation_less)
-        - next.cbegin());
-
-    if (side == 0) {
-      history.local().insert(make_time_id(t), order,
-                             NCd(quartic_side1(t.value())));
-    } else {
-      history.remote().insert(make_time_id(t), order,
-                              NCd(quartic_side2(t.value())));
-    }
-
-    gsl::at(next, side) += gsl::at(dt, side);
-
-    t = *std::min_element(next.cbegin(), next.cend(), simulation_less);
-
-    ASSERT(not simulation_less(next_check, t), "Screwed up arithmetic");
-    if (t == next_check) {
-      ab4.add_boundary_delta(&y, history, dt[0], quartic_coupling);
-      ab4.clean_boundary_history(make_not_null(&history));
-      CHECK(y() == approx(quartic_answer(t.value())));
-      if (t.is_at_slab_boundary()) {
-        break;
-      }
-      next_check += dt[0];
-    }
-  }
-}
-
-void check_lts_vts() {
-  const Slab slab(0., 1.);
-
-  const auto make_time_id = [](const Time& t) {
-    return TimeStepId(true, 0, t);
-  };
-
-  Time t = slab.start();
-
-  const size_t order = 4;
-  TimeSteppers::AdamsBashforth ab4(order);
-
-  TimeSteppers::BoundaryHistory<NCd, NCd, NCd> history{};
-  {
-    const Slab init_slab = slab.retreat();
-    const TimeDelta init_dt = init_slab.duration() / 4;
-
-    // clang-tidy misfeature: warns about boost internals here
-    for (int32_t step = 1; step <= 3; ++step) {  // NOLINT
-      // clang-tidy misfeature: warns about boost internals here
-      const Time now = t - step * init_dt;  // NOLINT
-      history.local().insert_initial(make_time_id(now), order,
-                                     NCd(quartic_side1(now.value())));
-      history.remote().insert_initial(make_time_id(now), order,
-                                      NCd(quartic_side2(now.value())));
-    }
-  }
-
-  std::array<std::deque<TimeDelta>, 2> dt{{
-      {slab.duration() / 2, slab.duration() / 4, slab.duration() / 4},
-      {slab.duration() / 6, slab.duration() / 6, slab.duration() * 2 / 9,
-            slab.duration() * 4 / 9}}};
-
-  NCd y(quartic_answer(t.value()));
-  Time next_check = t + dt[0][0];
-  std::array<Time, 2> next{{t, t}};
-  for (;;) {
-    const auto side = static_cast<size_t>(
-        std::min_element(next.cbegin(), next.cend()) - next.cbegin());
-
-    if (side == 0) {
-      history.local().insert(make_time_id(next[0]), order,
-                             NCd(quartic_side1(next[0].value())));
-    } else {
-      history.remote().insert(make_time_id(next[1]), order,
-                              NCd(quartic_side2(next[1].value())));
-    }
-
-    const TimeDelta this_dt = gsl::at(dt, side).front();
-    gsl::at(dt, side).pop_front();
-
-    gsl::at(next, side) += this_dt;
-
-    if (*std::min_element(next.cbegin(), next.cend()) == next_check) {
-      ab4.add_boundary_delta(&y, history, next_check - t, quartic_coupling);
-      ab4.clean_boundary_history(make_not_null(&history));
-      CHECK(y() == approx(quartic_answer(next_check.value())));
-      if (next_check.is_at_slab_boundary()) {
-        break;
-      }
-      t = next_check;
-      next_check += dt[0].front();
-    }
-  }
-}
-
 void test_neighbor_data_required() {
   // Test is order-independent
   const TimeSteppers::AdamsBashforth stepper(4);
@@ -360,53 +167,25 @@ void test_neighbor_data_required() {
 }
 }  // namespace
 
+// [[Timeout, 10]]
 SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.AdamsBashforth.Boundary",
                   "[Unit][Time]") {
   test_neighbor_data_required();
 
-  // No local stepping
   for (size_t order = 1; order < 9; ++order) {
-    INFO(order);
+    CAPTURE(order);
     const TimeSteppers::AdamsBashforth stepper(order);
-    for (size_t start_points = 0; start_points < order; ++start_points) {
-      INFO(start_points);
-      const double epsilon = std::max(std::pow(1e-3, start_points + 1), 1e-14);
-      TimeStepperTestUtils::lts::test_equal_rate(stepper, start_points + 1,
-                                                 start_points, epsilon, true);
-      TimeStepperTestUtils::lts::test_equal_rate(stepper, start_points + 1,
-                                                 start_points, epsilon, false);
+    TimeStepperTestUtils::lts::test_equal_rate(stepper);
+    TimeStepperTestUtils::lts::test_uncoupled(stepper, 1e-12);
+    TimeStepperTestUtils::lts::test_conservation(stepper);
+    // Only test convergence for low-order methods, since it's hard to
+    // find parameters where the high-order ones are in the convergent
+    // limit but not roundoff-dominated.
+    if (order < 5) {
+      TimeStepperTestUtils::lts::test_convergence(stepper, {20, 100}, 20);
+      TimeStepperTestUtils::lts::test_dense_convergence(stepper, {40, 200}, 40);
     }
   }
-
-  // Local stepping with constant step sizes
-  const Slab slab(0., 1.);
-  for (const auto& full : {slab.duration(), -slab.duration()}) {
-    do_lts_test({{full / 4, full / 4}});
-    do_lts_test({{full / 4, full / 8}});
-    do_lts_test({{full / 8, full / 4}});
-    do_lts_test({{full / 16, full / 4}});
-    do_lts_test({{full / 4, full / 16}});
-    do_lts_test({{full / 32, full / 4}});
-    do_lts_test({{full / 4, full / 32}});
-
-    // Non-nesting cases
-    do_lts_test({{full / 4, full / 6}});
-    do_lts_test({{full / 6, full / 4}});
-    do_lts_test({{full / 5, full / 7}});
-    do_lts_test({{full / 7, full / 5}});
-    do_lts_test({{full / 5, full / 13}});
-    do_lts_test({{full / 13, full / 5}});
-  }
-
-  // Local stepping with varying time steps
-  check_lts_vts();
-
-  // Dense output
-  for (size_t order = 1; order < 9; ++order) {
-    INFO(order);
-    TimeStepperTestUtils::lts::test_dense_output(
-        TimeSteppers::AdamsBashforth(order));
-  }
 }
 
 SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.AdamsBashforth.Reversal",
diff --git a/tests/Unit/Time/TimeSteppers/Test_AdamsLts.cpp b/tests/Unit/Time/TimeSteppers/Test_AdamsLts.cpp
new file mode 100644
index 000000000000..ad1227181325
--- /dev/null
+++ b/tests/Unit/Time/TimeSteppers/Test_AdamsLts.cpp
@@ -0,0 +1,759 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Utilities/StdHelpers.hpp"
+
+#include "Framework/TestingFramework.hpp"
+
+#include <array>
+#include <cstddef>
+#include <cstdint>
+#include <limits>
+#include <map>
+#include <optional>
+#include <ostream>
+#include <utility>
+#include <vector>
+
+#include "DataStructures/DataVector.hpp"
+#include "Time/ApproximateTime.hpp"
+#include "Time/BoundaryHistory.hpp"
+#include "Time/Slab.hpp"
+#include "Time/Time.hpp"
+#include "Time/TimeStepId.hpp"
+#include "Time/TimeSteppers/AdamsLts.hpp"
+#include "Utilities/ErrorHandling/Assert.hpp"
+#include "Utilities/Gsl.hpp"
+#include "Utilities/MakeWithValue.hpp"
+#include "Utilities/Rational.hpp"
+
+namespace {
+namespace adams_lts = TimeSteppers::adams_lts;
+
+void test_exact_substep_time() {
+  const Slab slab(1.2, 3.4);
+  const Time quarter = slab.start() + slab.duration() / 4;
+  const Time middle = slab.start() + slab.duration() / 2;
+  CHECK(adams_lts::exact_substep_time(TimeStepId(true, 1, middle)) == middle);
+  CHECK(adams_lts::exact_substep_time(TimeStepId(false, 1, middle)) == middle);
+  CHECK(adams_lts::exact_substep_time(TimeStepId(
+            true, 1, quarter, 1, slab.duration() / 4, middle.value())) ==
+        middle);
+  CHECK(adams_lts::exact_substep_time(TimeStepId(
+            false, 1, middle, 1, -slab.duration() / 4, quarter.value())) ==
+        quarter);
+#ifdef SPECTRE_DEBUG
+  CHECK_THROWS_WITH(
+      adams_lts::exact_substep_time(TimeStepId(
+          true, 1, slab.start(), 1, slab.duration(), middle.value())),
+      Catch::Matchers::ContainsSubstring("Substep not at expected time"));
+#endif  // SPECTRE_DEBUG
+}
+
+void test_lts_coefficients_struct() {
+  const Slab slab(0.0, 1.0);
+  const TimeStepId id1(true, 0, slab.start());
+  const TimeStepId id2 = id1.next_substep(slab.duration() / 2, 1.0);
+  const TimeStepId id3 = id2.next_step(slab.duration() / 2);
+
+  const adams_lts::LtsCoefficients coefs1{
+      {id1, id1, 2.0}, {id1, id2, 3.0}, {id1, id3, 4.0}};
+  const adams_lts::LtsCoefficients coefs2{
+      {id1, id1, 5.0}, {id1, id3, 7.0}, {id2, id1, 9.0}};
+
+  const adams_lts::LtsCoefficients sum{
+      {id1, id1, 7.0}, {id1, id2, 3.0}, {id1, id3, 11.0}, {id2, id1, 9.0}};
+  const adams_lts::LtsCoefficients difference{
+      {id1, id1, -3.0}, {id1, id2, 3.0}, {id1, id3, -3.0}, {id2, id1, -9.0}};
+
+  {
+    auto s = coefs1;
+    CHECK(&(s += coefs2) == &s);
+    CHECK(s == sum);
+  }
+  {
+    auto d = coefs1;
+    CHECK(&(d -= coefs2) == &d);
+    CHECK(d == difference);
+  }
+
+  CHECK(coefs1 + coefs2 == sum);
+  CHECK(coefs1 - coefs2 == difference);
+  CHECK(adams_lts::LtsCoefficients(coefs1) + coefs2 == sum);
+  CHECK(adams_lts::LtsCoefficients(coefs1) - coefs2 == difference);
+  CHECK(coefs1 + adams_lts::LtsCoefficients(coefs2) == sum);
+  CHECK(coefs1 - adams_lts::LtsCoefficients(coefs2) == difference);
+  CHECK(adams_lts::LtsCoefficients(coefs1) +
+            adams_lts::LtsCoefficients(coefs2) ==
+        sum);
+  CHECK(adams_lts::LtsCoefficients(coefs1) -
+            adams_lts::LtsCoefficients(coefs2) ==
+        difference);
+
+  adams_lts::LtsCoefficients allocated_coefs{};
+  for (size_t i = 0; i < adams_lts::lts_coefficients_static_size + 1; ++i) {
+    allocated_coefs.emplace_back(
+        id1,
+        TimeStepId(
+            true, static_cast<int64_t>(i),
+            Slab(static_cast<double>(i), static_cast<double>(i + 1)).start()),
+        i);
+  }
+
+  {
+    auto a = allocated_coefs;
+    const auto* const data = a.data();
+    const auto res = std::move(a) + allocated_coefs;
+    CHECK(res.data() == data);
+    CHECK(res == allocated_coefs + allocated_coefs);
+  }
+  {
+    auto a = allocated_coefs;
+    const auto* const data = a.data();
+    const auto res = std::move(a) - allocated_coefs;
+    CHECK(res.data() == data);
+    CHECK(res == allocated_coefs - allocated_coefs);
+  }
+  {
+    auto a = allocated_coefs;
+    const auto* const data = a.data();
+    const auto res = allocated_coefs + std::move(a);
+    CHECK(res.data() == data);
+    CHECK(res == allocated_coefs + allocated_coefs);
+  }
+  {
+    auto a = allocated_coefs;
+    const auto* const data = a.data();
+    const auto res = allocated_coefs - std::move(a);
+    CHECK(res.data() == data);
+    CHECK(res == allocated_coefs - allocated_coefs);
+  }
+  {
+    auto a = allocated_coefs;
+    auto b = allocated_coefs;
+    const auto* const a_data = a.data();
+    const auto* const b_data = a.data();
+    const auto res = std::move(a) + std::move(b);
+    CHECK((res.data() == a_data or res.data() == b_data));
+    CHECK(res == allocated_coefs + allocated_coefs);
+  }
+  {
+    auto a = allocated_coefs;
+    auto b = allocated_coefs;
+    const auto* const a_data = a.data();
+    const auto* const b_data = a.data();
+    const auto res = std::move(a) - std::move(b);
+    CHECK((res.data() == a_data or res.data() == b_data));
+    CHECK(res == allocated_coefs - allocated_coefs);
+  }
+}
+
+template <typename T>
+void test_apply_coefficients(const T& used_for_size) {
+  const Slab slab(1.2, 3.4);
+
+  TimeSteppers::BoundaryHistory<double, double, T> history{};
+  history.local().insert(TimeStepId(true, 0, slab.start()), 1, 1.0);
+  history.local().insert(TimeStepId(true, 0, slab.end()), 1, 10.0);
+  history.remote().insert(TimeStepId(true, 0, slab.start()), 1, 100.0);
+  history.remote().insert(TimeStepId(true, 0, slab.end()), 1, 10000.0);
+
+  const adams_lts::LtsCoefficients coefficients{
+      {history.local()[0], history.remote()[0], 1.0},
+      {history.local()[0], history.remote()[1], 2.0},
+      {history.local()[1], history.remote()[1], 3.0}};
+
+  const auto coupling = [&used_for_size](const double local,
+                                         const double remote) {
+    return make_with_value<T>(used_for_size, local * remote);
+  };
+
+  auto result = make_with_value<T>(used_for_size, 2.0);
+  adams_lts::apply_coefficients(make_not_null(&result), coefficients,
+                                history.evaluator(coupling));
+
+  CHECK(result == make_with_value<T>(used_for_size, 320102.0));
+}
+
+struct FakeId {
+  int step_time;
+  std::optional<int> step_size_if_substep{};
+  int64_t slab = 0;
+};
+
+bool operator<(const FakeId& a, const FakeId& b) {
+  return a.slab < b.slab or
+         (a.slab == b.slab and (a.step_time < b.step_time or
+                                (a.step_time == b.step_time and
+                                 not a.step_size_if_substep.has_value() and
+                                 b.step_size_if_substep.has_value())));
+}
+
+std::ostream& operator<<(std::ostream& s, const FakeId& id) {
+  using ::operator<<;
+  return s << "[" << id.step_time << " " << id.step_size_if_substep << " "
+           << id.slab << "]";
+}
+
+using ExpectedCoefficients = std::map<std::pair<FakeId, FakeId>, double>;
+
+ExpectedCoefficients flip_sides(const ExpectedCoefficients& coefs) {
+  ExpectedCoefficients flipped{};
+  for (const auto& coef : coefs) {
+    flipped.insert({{coef.first.second, coef.first.first}, coef.second});
+  }
+  return flipped;
+}
+
+namespace step_coefficients_detail {
+constexpr int max_time = 16;
+
+Time make_time(const int time) {
+  return {Slab(-max_time, max_time), Rational(time + max_time, 2 * max_time)};
+}
+
+ApproximateTime make_time(const double time) { return {time}; }
+
+TimeStepId make_id(const FakeId& fake_id) {
+  if (fake_id.step_size_if_substep.has_value()) {
+    ASSERT(*fake_id.step_size_if_substep > 0, "Zero step size");
+    const auto unit_step =
+        Slab(-max_time, max_time).duration() / (2 * max_time);
+    return {
+        true,
+        fake_id.slab,
+        make_time(fake_id.step_time),
+        1,
+        *fake_id.step_size_if_substep * unit_step,
+        static_cast<double>(fake_id.step_time + *fake_id.step_size_if_substep)};
+  } else {
+    return {true, fake_id.slab, make_time(fake_id.step_time)};
+  }
+}
+
+FakeId fake_from_id(const TimeStepId& id) {
+  return {static_cast<int>(id.step_time().value()),
+          id.substep() == 1
+              ? std::optional(static_cast<int>(id.step_size().value()))
+              : std::nullopt,
+          id.slab_number()};
+}
+
+// Map t -> -t/2.  As our test slab is symmetrical about 0, this maps
+// the fraction as [0, 1] -> [3/4, 1/4].
+Time alternate_time(const Time& time) {
+  return {time.slab(), Rational(3, 4) - time.fraction() / 2};
+}
+
+ApproximateTime alternate_time(const ApproximateTime& time) {
+  return {-0.5 * time.value()};
+}
+
+TimeStepId alternate_id(const TimeStepId& id) {
+  if (id.substep() == 0) {
+    return {not id.time_runs_forward(), id.slab_number(),
+            alternate_time(id.step_time())};
+  } else {
+    return {not id.time_runs_forward(),
+            id.slab_number(),
+            alternate_time(id.step_time()),
+            id.substep(),
+            id.step_size() / -2,
+            id.substep_time() / -2.0};
+  }
+}
+
+template <typename T>
+ExpectedCoefficients step_coefficients_impl(
+    const std::vector<FakeId>& local_steps,
+    const std::vector<FakeId>& remote_steps,
+    const adams_lts::AdamsScheme& local_scheme,
+    const adams_lts::AdamsScheme& remote_scheme,
+    const adams_lts::AdamsScheme& small_step_scheme, const int step_start,
+    const T step_end) {
+  const auto history_order = std::numeric_limits<size_t>::max();  // unused
+
+  TimeSteppers::BoundaryHistory<double, double, double> history{};
+  TimeSteppers::BoundaryHistory<double, double, double> alt_history{};
+
+  for (const auto& step : local_steps) {
+    const auto id = make_id(step);
+    history.local().insert(id, history_order, 0.0);
+    alt_history.local().insert(alternate_id(id), history_order, 0.0);
+  }
+  for (const auto& step : remote_steps) {
+    const auto id = make_id(step);
+    history.remote().insert(id, history_order, 0.0);
+    alt_history.remote().insert(alternate_id(id), history_order, 0.0);
+  }
+
+  const auto coefficients = adams_lts::lts_coefficients(
+      history.local(), history.remote(), make_time(step_start),
+      make_time(step_end), local_scheme, remote_scheme, small_step_scheme);
+
+  // Compare with step scaled by -1/2.
+  const auto alt_coefficients = adams_lts::lts_coefficients(
+      alt_history.local(), alt_history.remote(),
+      alternate_time(make_time(step_start)),
+      alternate_time(make_time(step_end)), local_scheme, remote_scheme,
+      small_step_scheme);
+  REQUIRE(coefficients.size() == alt_coefficients.size());
+  for (size_t i = 0; i < coefficients.size(); ++i) {
+    const auto& coef = coefficients[i];
+    const auto& alt_coef = alt_coefficients[i];
+    CHECK(get<0>(alt_coef) == alternate_id(get<0>(coef)));
+    CHECK(get<1>(alt_coef) == alternate_id(get<1>(coef)));
+    CHECK(get<2>(alt_coef) == approx(get<2>(coef) / -2.0));
+  }
+
+  ExpectedCoefficients coefficients_map{};
+  for (const auto& entry : coefficients) {
+    const auto insert_success = coefficients_map.insert(
+        {{fake_from_id(get<0>(entry)), fake_from_id(get<1>(entry))},
+         get<2>(entry)});
+    if (not insert_success.second) {
+      // Duplicate entry.
+      CAPTURE(entry);
+      CHECK(false);
+      insert_success.first->second += get<2>(entry);
+    }
+  }
+  return coefficients_map;
+}
+}  // namespace step_coefficients_detail
+
+ExpectedCoefficients dense_coefficients(
+    const std::vector<FakeId>& local_steps,
+    const std::vector<FakeId>& remote_steps,
+    const adams_lts::AdamsScheme& local_scheme,
+    const adams_lts::AdamsScheme& remote_scheme,
+    const adams_lts::AdamsScheme& small_step_scheme, const int step_start,
+    const double step_end) {
+  auto result = step_coefficients_detail::step_coefficients_impl(
+      local_steps, remote_steps, local_scheme, remote_scheme, small_step_scheme,
+      step_start, step_end);
+  // NOLINTNEXTLINE(readability-suspicious-call-argument)
+  const auto reversed = step_coefficients_detail::step_coefficients_impl(
+      remote_steps, local_steps, remote_scheme, local_scheme, small_step_scheme,
+      step_start, step_end);
+  CHECK_ITERABLE_APPROX(result, flip_sides(reversed));
+  return result;
+}
+
+ExpectedCoefficients step_coefficients(
+    const std::vector<FakeId>& local_steps,
+    const std::vector<FakeId>& remote_steps,
+    const adams_lts::AdamsScheme& local_scheme,
+    const adams_lts::AdamsScheme& remote_scheme,
+    const adams_lts::AdamsScheme& small_step_scheme, const int step_start,
+    const int step_end, const bool also_check_dense = true) {
+  auto result = step_coefficients_detail::step_coefficients_impl(
+      local_steps, remote_steps, local_scheme, remote_scheme, small_step_scheme,
+      step_start, step_end);
+  // NOLINTNEXTLINE(readability-suspicious-call-argument)
+  const auto reversed = step_coefficients_detail::step_coefficients_impl(
+      remote_steps, local_steps, remote_scheme, local_scheme, small_step_scheme,
+      step_start, step_end);
+  CHECK_ITERABLE_APPROX(result, flip_sides(reversed));
+  if (also_check_dense) {
+    const auto dense = dense_coefficients(
+        local_steps, remote_steps, local_scheme, remote_scheme,
+        small_step_scheme, step_start, static_cast<double>(step_end));
+    CHECK_ITERABLE_APPROX(result, dense);
+  }
+  return result;
+}
+
+void test_lts_coefficients() {
+  // AB2 [0 1] step to 2
+  const std::array coefs_ab2{-1.0 / 2.0, 3.0 / 2.0};
+  // AB2 [0 1] step to 3/2
+  const std::array coefs_ab2_32{-1.0 / 8.0, 5.0 / 8.0};
+  // AB3 [0 1 2] step to 3
+  const std::array coefs_ab3{5.0 / 12.0, -4.0 / 3.0, 23.0 / 12.0};
+  // AB3 [0 1 2] step to 5/2
+  const std::array coefs_ab3_52{1.0 / 12.0, -7.0 / 24.0, 17.0 / 24.0};
+  // AM3 [0 1 (2)] step to 2
+  const std::array coefs_am3{-1.0 / 12.0, 2.0 / 3.0, 5.0 / 12.0};
+  // AM3 [0 1 (2)] dense to 3/2
+  const std::array coefs_am3_32{-1.0 / 24.0, 11.0 / 24.0, 1.0 / 12.0};
+
+  {
+    INFO("AB GTS order 1");
+    const adams_lts::AdamsScheme ab1{adams_lts::SchemeType::Explicit, 1};
+    const std::vector<FakeId> steps{{0}};
+    // clang-format off
+    const ExpectedCoefficients expected{
+        {{{0}, {0}}, 1.0}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps, steps, ab1, ab1, ab1, 0, 1), expected);
+    const ExpectedCoefficients expected_dense_12{
+        {{{0}, {0}}, 0.5}};
+    CHECK_ITERABLE_APPROX(
+        dense_coefficients(steps, steps, ab1, ab1, ab1, 0, 0.5),
+        expected_dense_12);
+    // clang-format on
+  }
+
+  {
+    INFO("AB GTS order 3");
+    const adams_lts::AdamsScheme ab3{adams_lts::SchemeType::Explicit, 3};
+    const std::vector<FakeId> steps{{0}, {1}, {2}};
+    // clang-format off
+    const ExpectedCoefficients expected{
+        {{{0}, {0}}, coefs_ab3[0]},
+        {{{1}, {1}}, coefs_ab3[1]},
+        {{{2}, {2}}, coefs_ab3[2]}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps, steps, ab3, ab3, ab3, 2, 3), expected);
+    const ExpectedCoefficients expected_dense_52{
+        {{{0}, {0}}, coefs_ab3_52[0]},
+        {{{1}, {1}}, coefs_ab3_52[1]},
+        {{{2}, {2}}, coefs_ab3_52[2]}};
+    CHECK_ITERABLE_APPROX(
+        dense_coefficients(steps, steps, ab3, ab3, ab3, 2, 2.5),
+        expected_dense_52);
+    // clang-format on
+  }
+
+  {
+    INFO("AM GTS order 3");
+    const adams_lts::AdamsScheme ab2{adams_lts::SchemeType::Explicit, 2};
+    const adams_lts::AdamsScheme am3{adams_lts::SchemeType::Implicit, 3};
+    const std::vector<FakeId> steps{{0}, {1}, {1, 1}};
+    // clang-format off
+    const ExpectedCoefficients expected{
+        {{{0}, {0}}, coefs_am3[0]},
+        {{{1}, {1}}, coefs_am3[1]},
+        {{{1, 1}, {1, 1}}, coefs_am3[2]}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps, steps, am3, am3, am3, 1, 2), expected);
+    const ExpectedCoefficients expected_predictor{
+        {{{0}, {0}}, coefs_ab2[0]},
+        {{{1}, {1}}, coefs_ab2[1]}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps, steps, ab2, ab2, ab2, 1, 2),
+        expected_predictor);
+    const ExpectedCoefficients expected_dense_32_nonmonotonic{
+        {{{0}, {0}}, coefs_am3_32[0]},
+        {{{1}, {1}}, coefs_am3_32[1]},
+        {{{1, 1}, {1, 1}}, coefs_am3_32[2]}};
+    CHECK_ITERABLE_APPROX(
+        dense_coefficients(steps, steps, am3, am3, am3, 1, 1.5),
+        expected_dense_32_nonmonotonic);
+    const ExpectedCoefficients expected_dense_32_monotonic{
+        {{{0}, {0}}, coefs_ab2_32[0]},
+        {{{1}, {1}}, coefs_ab2_32[1]}};
+    CHECK_ITERABLE_APPROX(
+        dense_coefficients(steps, steps, ab2, ab2, ab2, 1, 1.5),
+        expected_dense_32_monotonic);
+    // clang-format on
+  }
+
+  {
+    INFO("AB single-side order 3");
+    const adams_lts::AdamsScheme ab1{adams_lts::SchemeType::Explicit, 1};
+    const adams_lts::AdamsScheme ab3{adams_lts::SchemeType::Explicit, 3};
+    const std::vector<FakeId> steps{{0}, {1}, {2}};
+    const std::vector<FakeId> single_step{{0}};
+    // clang-format off
+    const ExpectedCoefficients expected{
+        {{{0}, {0}}, coefs_ab3[0]},
+        {{{1}, {0}}, coefs_ab3[1]},
+        {{{2}, {0}}, coefs_ab3[2]}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps, single_step, ab3, ab1, ab3, 2, 3), expected);
+    const ExpectedCoefficients expected_dense_52{
+        {{{0}, {0}}, coefs_ab3_52[0]},
+        {{{1}, {0}}, coefs_ab3_52[1]},
+        {{{2}, {0}}, coefs_ab3_52[2]}};
+    CHECK_ITERABLE_APPROX(
+        dense_coefficients(steps, single_step, ab3, ab1, ab3, 2, 2.5),
+        expected_dense_52);
+    // clang-format on
+  }
+
+  {
+    INFO("AM single-side order 3");
+    const adams_lts::AdamsScheme ab1{adams_lts::SchemeType::Explicit, 1};
+    const adams_lts::AdamsScheme ab2{adams_lts::SchemeType::Explicit, 2};
+    const adams_lts::AdamsScheme am3{adams_lts::SchemeType::Implicit, 3};
+    const std::vector<FakeId> steps{{0}, {1}, {1, 1}};
+    const std::vector<FakeId> single_step{{0}};
+    // clang-format off
+    const ExpectedCoefficients expected{
+        {{{0}, {0}}, coefs_am3[0]},
+        {{{1}, {0}}, coefs_am3[1]},
+        {{{1, 1}, {0}}, coefs_am3[2]}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps, single_step, am3, ab1, am3, 1, 2),
+        expected);
+    const ExpectedCoefficients expected_predictor{
+        {{{0}, {0}}, coefs_ab2[0]},
+        {{{1}, {0}}, coefs_ab2[1]}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps, single_step, ab2, ab1, ab2, 1, 2),
+        expected_predictor);
+    // Nonmonotonic requires additional future data.
+    const ExpectedCoefficients expected_dense_32_monotonic{
+        {{{0}, {0}}, coefs_ab2_32[0]},
+        {{{1}, {0}}, coefs_ab2_32[1]}};
+    CHECK_ITERABLE_APPROX(
+        dense_coefficients(steps, single_step, ab2, ab1, ab2, 1, 1.5),
+        expected_dense_32_monotonic);
+    // clang-format on
+  }
+
+  {
+    INFO("AB 2:1 order 3");
+    // -8          -4           0           4
+    //             -4    -2     0     2     4
+    const adams_lts::AdamsScheme ab3{adams_lts::SchemeType::Explicit, 3};
+    const std::vector<FakeId> steps_large{{-8}, {-4}, {0}};
+    const std::vector<FakeId> steps_small{{-4}, {-2}, {0}, {2}};
+    // clang-format off
+    const ExpectedCoefficients expected_large{
+        {{{0}, {2}}, 115.0 / 16.0},
+        {{{0}, {0}}, 7.0 / 6.0},
+        {{{0}, {-2}}, -11.0 / 16.0},
+        {{{-4}, {2}}, -115.0 / 24.0},
+        {{{-4}, {-2}}, -11.0 / 8.0},
+        {{{-4}, {-4}}, 5.0 / 6.0},
+        {{{-8}, {2}}, 23.0 / 16.0},
+        {{{-8}, {-2}}, 11.0 / 48.0}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps_large, steps_small, ab3, ab3, ab3, 0, 4, false),
+        expected_large);
+    const ExpectedCoefficients expected_small_1{
+        {{{0}, {0}}, 23.0 / 6.0},
+        {{{-2}, {0}}, -1.0},
+        {{{-2}, {-4}}, -2.0},
+        {{{-4}, {-4}}, 5.0 / 6.0},
+        {{{-2}, {-8}}, 1.0 / 3.0}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps_small, steps_large, ab3, ab3, ab3, 0, 2),
+        expected_small_1);
+    const ExpectedCoefficients expected_small_2{
+        {{{2}, {0}}, 115.0 / 16.0},
+        {{{0}, {0}}, -8.0 / 3.0},
+        {{{-2}, {0}}, 5.0 / 16.0},
+        {{{2}, {-4}}, -115.0 / 24.0},
+        {{{-2}, {-4}}, 5.0 / 8.0},
+        {{{2}, {-8}}, 23.0 / 16.0},
+        {{{-2}, {-8}}, -5.0 / 48.0}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps_small, steps_large, ab3, ab3, ab3, 2, 4),
+        expected_small_2);
+    // clang-format on
+  }
+
+  {
+    INFO("AB LTS -> GTS order 2");
+    // -2     0  1
+    //    -1  0  1
+    const adams_lts::AdamsScheme ab2{adams_lts::SchemeType::Explicit, 2};
+    const std::vector<FakeId> steps_large{{-2}, {0}};
+    const std::vector<FakeId> steps_small{{-1}, {0}};
+    // clang-format off
+    const ExpectedCoefficients expected{
+        {{{0}, {0}}, 3.0 / 2.0},
+        {{{0}, {-1}}, -1.0 / 4.0},
+        {{{-2}, {-1}}, -1.0 / 4.0}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps_large, steps_small, ab2, ab2, ab2, 0, 1),
+        expected);
+    // clang-format on
+  }
+
+  {
+    INFO("AB 3:1 order 2");
+    // -3        0        3
+    //       -1  0  1  2  3
+    const adams_lts::AdamsScheme ab2{adams_lts::SchemeType::Explicit, 2};
+    const std::vector<FakeId> steps_large{{-3}, {0}};
+    const std::vector<FakeId> steps_small{{-1}, {0}, {1}, {2}};
+    // clang-format off
+    const ExpectedCoefficients expected_large{
+        {{{-3}, {-1}}, -1.0 / 6.0},
+        {{{-3}, {1}}, -1.0 / 3.0},
+        {{{-3}, {2}}, -1.0},
+        {{{0}, {-1}}, -1.0 / 3.0},
+        {{{0}, {0}}, 1.0},
+        {{{0}, {1}}, 4.0 / 3.0},
+        {{{0}, {2}}, 5.0 / 2.0}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps_large, steps_small, ab2, ab2, ab2, 0, 3, false),
+        expected_large);
+    const ExpectedCoefficients expected_small_1{
+        {{{-1}, {-3}}, -1.0 / 6.0},
+        {{{-1}, {0}}, -1.0 / 3.0},
+        {{{0}, {0}}, 3.0 / 2.0}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps_small, steps_large, ab2, ab2, ab2, 0, 1),
+        expected_small_1);
+    const ExpectedCoefficients expected_small_2{
+        {{{0}, {0}}, -1.0 / 2.0},
+        {{{1}, {-3}}, -1.0 / 2.0},
+        {{{1}, {0}}, 2.0}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps_small, steps_large, ab2, ab2, ab2, 1, 2),
+        expected_small_2);
+    const ExpectedCoefficients expected_small_3{
+        {{{1}, {-3}}, 1.0 / 6.0},
+        {{{1}, {0}}, -2.0 / 3.0},
+        {{{2}, {-3}}, -1.0},
+        {{{2}, {0}}, 5.0 / 2.0}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps_small, steps_large, ab2, ab2, ab2, 2, 3),
+        expected_small_3);
+
+    const ExpectedCoefficients expected_dense_12{
+        {{{-3}, {-1}}, -1.0 / 24.0},
+        {{{0}, {-1}}, -1.0 / 12.0},
+        {{{0}, {0}}, 5.0 / 8.0}};
+    CHECK_ITERABLE_APPROX(
+        dense_coefficients(steps_large, steps_small, ab2, ab2, ab2, 0, 0.5),
+        expected_dense_12);
+    const ExpectedCoefficients expected_dense_32{
+        {{{0}, {0}}, -1.0 / 8.0},
+        {{{-3}, {1}}, -5.0 / 24.0},
+        {{{0}, {1}}, 5.0 / 6.0}};
+    CHECK_ITERABLE_APPROX(
+        dense_coefficients(steps_large, steps_small, ab2, ab2, ab2, 1, 1.5),
+        expected_dense_32);
+    const ExpectedCoefficients expected_dense_52{
+        {{{-3}, {1}}, 1.0 / 24.0},
+        {{{0}, {1}}, -1.0 / 6.0},
+        {{{-3}, {2}}, -5.0 / 12.0},
+        {{{0}, {2}}, 25.0 / 24.0}};
+    CHECK_ITERABLE_APPROX(
+        dense_coefficients(steps_large, steps_small, ab2, ab2, ab2, 2, 2.5),
+        expected_dense_52);
+    // clang-format on
+  }
+
+  {
+    INFO("AB unaligned order 2");
+    // 0  1     3  4     6
+    // 0     2  3     5  6
+    const adams_lts::AdamsScheme ab2{adams_lts::SchemeType::Explicit, 2};
+    const std::vector<FakeId> steps_a{{1}, {3}, {4}};
+    const std::vector<FakeId> steps_b{{2}, {3}, {5}};
+    // clang-format off
+    const ExpectedCoefficients expected_a_1{
+        {{{1}, {2}}, -1.0 / 4.0},
+        {{{3}, {2}}, -1.0 / 4.0},
+        {{{3}, {3}}, 3.0 / 2.0}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps_a, steps_b, ab2, ab2, ab2, 3, 4),
+        expected_a_1);
+    const ExpectedCoefficients expected_a_2{
+        {{{3}, {3}}, -1.0 / 2.0},
+        {{{3}, {5}}, -3.0 / 2.0},
+        {{{4}, {2}}, -3.0 / 2.0},
+        {{{4}, {3}}, 11.0 / 4.0},
+        {{{4}, {5}}, 11.0 / 4.0}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps_a, steps_b, ab2, ab2, ab2, 4, 6, false),
+        expected_a_2);
+    const ExpectedCoefficients expected_b_1{
+        {{{2}, {1}}, -1.0 / 4.0},
+        {{{2}, {3}}, -1.0 / 4.0},
+        {{{2}, {4}}, -3.0 / 2.0},
+        {{{3}, {3}}, 1.0},
+        {{{3}, {4}}, 3.0}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps_b, steps_a, ab2, ab2, ab2, 3, 5, false),
+        expected_b_1);
+    const ExpectedCoefficients expected_b_2{
+        {{{3}, {4}}, -1.0 / 4.0},
+        {{{5}, {3}}, -3.0 / 2.0},
+        {{{5}, {4}}, 11.0 / 4.0}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps_b, steps_a, ab2, ab2, ab2, 5, 6),
+        expected_b_2);
+    // clang-format on
+  }
+
+  {
+    INFO("AM 2:1 order 2");
+    // 0     2
+    // 0  1  2
+    const adams_lts::AdamsScheme ab1{adams_lts::SchemeType::Explicit, 1};
+    const adams_lts::AdamsScheme am2{adams_lts::SchemeType::Implicit, 2};
+    const std::vector<FakeId> steps_large{{0}, {0, 2}};
+    const std::vector<FakeId> steps_small{{0}, {0, 1}, {1}, {1, 1}};
+    const std::vector<FakeId> steps_small_for_nonmonotonic_predictor{{0}};
+    // clang-format off
+    // Monotonic predictor requires two calls, both of which are
+    // small-side tests below.
+    const ExpectedCoefficients expected_large_predictor_nonmonotonic{
+        {{{0}, {0}}, 2.0}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps_large, steps_small_for_nonmonotonic_predictor,
+                          ab1, ab1, ab1, 0, 2),
+        expected_large_predictor_nonmonotonic);
+    const ExpectedCoefficients expected_large_corrector{
+        {{{0}, {0}}, 1.0 / 2.0},
+        {{{0}, {0, 1}}, 1.0 / 4.0},
+        {{{0, 2}, {0, 1}}, 1.0 / 4.0},
+        {{{0}, {1}}, 1.0 / 4.0},
+        {{{0, 2}, {1}}, 1.0 / 4.0},
+        {{{0, 2}, {1, 1}}, 1.0 / 2.0}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps_large, steps_small, am2, am2, am2, 0, 2, false),
+        expected_large_corrector);
+    const ExpectedCoefficients expected_small_predictor_1{
+        {{{0}, {0}}, 1.0}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps_small, steps_large, ab1, ab1, ab1, 0, 1),
+        expected_small_predictor_1);
+    const ExpectedCoefficients expected_small_corrector_1{
+        {{{0}, {0}}, 1.0 / 2.0},
+        {{{0, 1}, {0}}, 1.0 / 4.0},
+        {{{0, 1}, {0, 2}}, 1.0 / 4.0}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps_small, steps_large, am2, am2, am2, 0, 1),
+        expected_small_corrector_1);
+    const ExpectedCoefficients expected_small_predictor_2{
+        {{{1}, {0}}, 1.0}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps_small, steps_large, ab1, ab1, ab1, 1, 2),
+        expected_small_predictor_2);
+    const ExpectedCoefficients expected_small_corrector_2{
+        {{{1}, {0}}, 1.0 / 4.0},
+        {{{1}, {0, 2}}, 1.0 / 4.0},
+        {{{1, 1}, {0, 2}}, 1.0 / 2.0}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps_small, steps_large, am2, am2, am2, 1, 2),
+        expected_small_corrector_2);
+    // clang-format on
+  }
+
+  {
+    INFO("AB self-start order 4");
+    // 0 1 2 3 -- 0 1
+    // 0 1 2 3 -- 0 1
+    const adams_lts::AdamsScheme ab4{adams_lts::SchemeType::Explicit, 4};
+    const std::vector<FakeId> steps{{2, {}, -1}, {3, {}, -1}, {0}, {1}};
+    // clang-format off
+    const ExpectedCoefficients expected{
+        {{{2, {}, -1}, {2, {}, -1}}, 13.0 / 24.0},
+        {{{3, {}, -1}, {3, {}, -1}}, -1.0 / 24.0},
+        {{{0}, {0}}, -1.0 / 24.0},
+        {{{1}, {1}}, 13.0 / 24.0}};
+    CHECK_ITERABLE_APPROX(
+        step_coefficients(steps, steps, ab4, ab4, ab4, 1, 2),
+        expected);
+    // clang-format on
+  }
+}
+
+SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.AdamsLts", "[Unit][Time]") {
+  test_exact_substep_time();
+  test_lts_coefficients_struct();
+  test_apply_coefficients(0.0);
+  test_apply_coefficients(DataVector(5, 0.0));
+  test_lts_coefficients();
+}
+}  // namespace
diff --git a/tests/Unit/Time/TimeSteppers/Test_AdamsMoultonPc.cpp b/tests/Unit/Time/TimeSteppers/Test_AdamsMoultonPc.cpp
index 6c4b63ab5d86..1c10d4a0821d 100644
--- a/tests/Unit/Time/TimeSteppers/Test_AdamsMoultonPc.cpp
+++ b/tests/Unit/Time/TimeSteppers/Test_AdamsMoultonPc.cpp
@@ -6,7 +6,9 @@
 #include <algorithm>
 #include <cmath>
 #include <cstddef>
+#include <cstdint>
 #include <utility>
+#include <vector>
 
 #include "Framework/TestCreation.hpp"
 #include "Framework/TestHelpers.hpp"
@@ -18,7 +20,9 @@
 #include "Time/TimeStepId.hpp"
 #include "Time/TimeSteppers/AdamsMoultonPc.hpp"
 #include "Time/TimeSteppers/TimeStepper.hpp"
+#include "Utilities/Gsl.hpp"
 #include "Utilities/Literals.hpp"
+#include "Utilities/Rational.hpp"
 
 namespace {
 template <bool Monotonic>
@@ -139,14 +143,264 @@ void test_am() {
   }
 }
 
-// [[Timeout, 10]]
+template <bool Monotonic>
+void test_neighbor_data_required(const bool time_runs_forward) {
+  // Test is order-independent
+  const TimeSteppers::AdamsMoultonPc<Monotonic> stepper(4);
+  const Slab slab(0.0, 1.0);
+  const auto long_step = (time_runs_forward ? 1 : -1) * slab.duration() / 2;
+  const auto short_step = long_step / 2;
+
+  const TimeStepId start(time_runs_forward, 1,
+                         time_runs_forward ? slab.start() : slab.end());
+  const TimeStepId long_sub = start.next_substep(long_step, 1.0);
+  const TimeStepId short_sub1 = start.next_substep(short_step, 1.0);
+  const TimeStepId middle(time_runs_forward, 1, start.step_time() + short_step);
+  const TimeStepId short_sub2 = middle.next_substep(short_step, 1.0);
+  const TimeStepId end(time_runs_forward, 1, start.step_time() + long_step);
+
+  const double middle_time = middle.step_time().value();
+
+  if constexpr (Monotonic) {
+    CHECK(not stepper.neighbor_data_required(start, start));
+
+    CHECK(stepper.neighbor_data_required(long_sub, start));
+    CHECK(not stepper.neighbor_data_required(long_sub, long_sub));
+    CHECK(stepper.neighbor_data_required(long_sub, short_sub1));
+    CHECK(stepper.neighbor_data_required(long_sub, middle));
+    CHECK(not stepper.neighbor_data_required(long_sub, short_sub2));
+
+    CHECK(stepper.neighbor_data_required(short_sub1, start));
+    CHECK(not stepper.neighbor_data_required(short_sub1, long_sub));
+    CHECK(not stepper.neighbor_data_required(short_sub1, short_sub1));
+    CHECK(not stepper.neighbor_data_required(short_sub1, middle));
+
+    CHECK(stepper.neighbor_data_required(middle, start));
+    CHECK(not stepper.neighbor_data_required(middle, long_sub));
+    CHECK(stepper.neighbor_data_required(middle, short_sub1));
+    CHECK(not stepper.neighbor_data_required(middle, middle));
+
+    CHECK(stepper.neighbor_data_required(middle_time, start));
+    CHECK(not stepper.neighbor_data_required(middle_time, long_sub));
+    CHECK(stepper.neighbor_data_required(middle_time, short_sub1));
+    CHECK(stepper.neighbor_data_required(middle_time, middle));
+    CHECK(not stepper.neighbor_data_required(middle_time, short_sub2));
+
+    CHECK(stepper.neighbor_data_required(short_sub2, start));
+    CHECK(not stepper.neighbor_data_required(short_sub2, long_sub));
+    CHECK(stepper.neighbor_data_required(short_sub2, short_sub1));
+    CHECK(stepper.neighbor_data_required(short_sub2, middle));
+    CHECK(not stepper.neighbor_data_required(short_sub2, short_sub2));
+
+    CHECK(stepper.neighbor_data_required(end, start));
+    CHECK(stepper.neighbor_data_required(end, long_sub));
+    CHECK(stepper.neighbor_data_required(end, short_sub1));
+    CHECK(stepper.neighbor_data_required(end, middle));
+    CHECK(stepper.neighbor_data_required(end, short_sub2));
+    CHECK(not stepper.neighbor_data_required(end, end));
+  } else {
+    CHECK(not stepper.neighbor_data_required(start, start));
+
+    CHECK(stepper.neighbor_data_required(long_sub, start));
+    CHECK(not stepper.neighbor_data_required(long_sub, long_sub));
+    CHECK(not stepper.neighbor_data_required(long_sub, short_sub1));
+    CHECK(not stepper.neighbor_data_required(long_sub, middle));
+    CHECK(not stepper.neighbor_data_required(long_sub, short_sub2));
+
+    CHECK(stepper.neighbor_data_required(short_sub1, start));
+    CHECK(not stepper.neighbor_data_required(short_sub1, long_sub));
+    CHECK(not stepper.neighbor_data_required(short_sub1, short_sub1));
+    CHECK(not stepper.neighbor_data_required(short_sub1, middle));
+
+    CHECK(stepper.neighbor_data_required(middle, start));
+    CHECK(stepper.neighbor_data_required(middle, long_sub));
+    CHECK(stepper.neighbor_data_required(middle, short_sub1));
+    CHECK(not stepper.neighbor_data_required(middle, middle));
+
+    CHECK(stepper.neighbor_data_required(middle_time, start));
+    CHECK(stepper.neighbor_data_required(middle_time, long_sub));
+    CHECK(stepper.neighbor_data_required(middle_time, short_sub1));
+    CHECK(not stepper.neighbor_data_required(middle_time, middle));
+
+    CHECK(stepper.neighbor_data_required(short_sub2, start));
+    CHECK(stepper.neighbor_data_required(short_sub2, long_sub));
+    CHECK(stepper.neighbor_data_required(short_sub2, short_sub1));
+    CHECK(stepper.neighbor_data_required(short_sub2, middle));
+    CHECK(not stepper.neighbor_data_required(short_sub2, short_sub2));
+
+    CHECK(stepper.neighbor_data_required(end, start));
+    CHECK(stepper.neighbor_data_required(end, long_sub));
+    CHECK(stepper.neighbor_data_required(end, short_sub1));
+    CHECK(stepper.neighbor_data_required(end, middle));
+    CHECK(stepper.neighbor_data_required(end, short_sub2));
+    CHECK(not stepper.neighbor_data_required(end, end));
+  }
+}
+
+template <bool Monotonic>
+void test_boundary() {
+  test_neighbor_data_required<Monotonic>(true);
+  test_neighbor_data_required<Monotonic>(false);
+
+  for (size_t order = 2; order < 9; ++order) {
+    CAPTURE(order);
+    const TimeSteppers::AdamsMoultonPc<Monotonic> stepper(order);
+    TimeStepperTestUtils::lts::test_equal_rate(stepper);
+    TimeStepperTestUtils::lts::test_uncoupled(stepper, 1e-12);
+    TimeStepperTestUtils::lts::test_conservation(stepper);
+    // Only test convergence for low-order methods, since it's hard to
+    // find parameters where the high-order ones are in the convergent
+    // limit but not roundoff-dominated.
+    if (order < 5) {
+      TimeStepperTestUtils::lts::test_convergence(stepper, {20, 100}, 20);
+      TimeStepperTestUtils::lts::test_dense_convergence(stepper, {50, 170}, 15);
+    }
+  }
+}
+
+TimeStepId make_id(const Rational& time_fraction, const uint64_t substep,
+                   const Rational& step_fraction = 0) {
+  ASSERT(substep == 0 or step_fraction != 0, "Need step size for a substep.");
+  const Slab slab(0.0, 1.0);
+  const TimeStepId step_id(true, 0, Time(slab, time_fraction));
+  if (substep == 0) {
+    return step_id;
+  } else {
+    return step_id.next_substep(TimeDelta(slab, step_fraction), 1.0);
+  }
+}
+
+using Counts = std::vector<std::pair<Rational, size_t>>;
+
+Counts substep_counts(const TimeSteppers::ConstBoundaryHistoryTimes& times) {
+  std::vector<std::pair<Rational, size_t>> counts{};
+  counts.reserve(times.size());
+  for (const auto& step : times) {
+    counts.emplace_back(step.step_time().fraction(),
+                        times.number_of_substeps(step));
+  }
+  return counts;
+}
+
+void test_non_monotonic_boundary_cleaning() {
+  const TimeSteppers::AdamsMoultonPc<false> am3(3);
+  TimeSteppers::BoundaryHistory<double, double, double> history{};
+
+  // Equal rate
+  history.local().insert(make_id({0, 8}, 0), 3, 0.0);
+  history.local().insert(make_id({1, 8}, 0), 3, 0.0);
+  history.remote().insert(make_id({0, 8}, 0), 3, 0.0);
+  history.remote().insert(make_id({1, 8}, 0), 3, 0.0);
+  am3.clean_boundary_history(make_not_null(&history));
+  CHECK(substep_counts(history.local()) == Counts{{{0, 8}, 1}, {{1, 8}, 1}});
+  CHECK(substep_counts(history.remote()) == Counts{{{0, 8}, 1}, {{1, 8}, 1}});
+  history.local().insert(make_id({1, 8}, 1, {1, 8}), 3, 0.0);
+  history.remote().insert(make_id({1, 8}, 1, {1, 8}), 3, 0.0);
+  am3.clean_boundary_history(make_not_null(&history));
+  CHECK(substep_counts(history.local()) == Counts{{{1, 8}, 1}});
+  CHECK(substep_counts(history.remote()) == Counts{{{1, 8}, 1}});
+
+  // Local finer
+  history.local().insert(make_id({2, 8}, 0), 3, 0.0);
+  history.remote().insert(make_id({2, 8}, 0), 3, 0.0);
+  am3.clean_boundary_history(make_not_null(&history));
+  CHECK(substep_counts(history.local()) == Counts{{{1, 8}, 1}, {{2, 8}, 1}});
+  CHECK(substep_counts(history.remote()) == Counts{{{1, 8}, 1}, {{2, 8}, 1}});
+  history.local().insert(make_id({2, 8}, 1, {1, 8}), 3, 0.0);
+  history.remote().insert(make_id({2, 8}, 1, {2, 8}), 3, 0.0);
+  am3.clean_boundary_history(make_not_null(&history));
+  CHECK(substep_counts(history.local()) == Counts{{{2, 8}, 1}});
+  CHECK(substep_counts(history.remote()) == Counts{{{1, 8}, 1}, {{2, 8}, 2}});
+  history.local().insert(make_id({3, 8}, 0), 3, 0.0);
+  am3.clean_boundary_history(make_not_null(&history));
+  CHECK(substep_counts(history.local()) == Counts{{{2, 8}, 1}, {{3, 8}, 1}});
+  CHECK(substep_counts(history.remote()) == Counts{{{1, 8}, 1}, {{2, 8}, 2}});
+  history.local().insert(make_id({3, 8}, 1, {1, 8}), 3, 0.0);
+  am3.clean_boundary_history(make_not_null(&history));
+  CHECK(substep_counts(history.local()) == Counts{{{3, 8}, 1}});
+  CHECK(substep_counts(history.remote()) == Counts{{{2, 8}, 1}});
+
+  // Local coarser
+  history.local().insert(make_id({4, 8}, 0), 3, 0.0);
+  history.remote().insert(make_id({4, 8}, 0), 3, 0.0);
+  am3.clean_boundary_history(make_not_null(&history));
+  CHECK(substep_counts(history.local()) == Counts{{{3, 8}, 1}, {{4, 8}, 1}});
+  CHECK(substep_counts(history.remote()) == Counts{{{2, 8}, 1}, {{4, 8}, 1}});
+  history.local().insert(make_id({4, 8}, 1, {2, 8}), 3, 0.0);
+  history.remote().insert(make_id({4, 8}, 1, {1, 8}), 3, 0.0);
+  history.remote().insert(make_id({5, 8}, 0), 3, 0.0);
+  history.remote().insert(make_id({5, 8}, 1, {1, 8}), 3, 0.0);
+  am3.clean_boundary_history(make_not_null(&history));
+  CHECK(substep_counts(history.local()) == Counts{{{4, 8}, 1}});
+  CHECK(substep_counts(history.remote()) == Counts{{{5, 8}, 1}});
+}
+
+void test_monotonic_boundary_cleaning() {
+  const TimeSteppers::AdamsMoultonPc<true> am3(3);
+  TimeSteppers::BoundaryHistory<double, double, double> history{};
+
+  // Equal rate
+  history.local().insert(make_id({0, 8}, 0), 3, 0.0);
+  history.local().insert(make_id({1, 8}, 0), 3, 0.0);
+  history.remote().insert(make_id({0, 8}, 0), 3, 0.0);
+  history.remote().insert(make_id({1, 8}, 0), 3, 0.0);
+  am3.clean_boundary_history(make_not_null(&history));
+  CHECK(substep_counts(history.local()) == Counts{{{0, 8}, 1}, {{1, 8}, 1}});
+  CHECK(substep_counts(history.remote()) == Counts{{{0, 8}, 1}, {{1, 8}, 1}});
+  history.local().insert(make_id({1, 8}, 1, {1, 8}), 3, 0.0);
+  history.remote().insert(make_id({1, 8}, 1, {1, 8}), 3, 0.0);
+  am3.clean_boundary_history(make_not_null(&history));
+  CHECK(substep_counts(history.local()) == Counts{{{1, 8}, 1}});
+  CHECK(substep_counts(history.remote()) == Counts{{{1, 8}, 1}});
+
+  // Local finer
+  history.local().insert(make_id({2, 8}, 0), 3, 0.0);
+  history.remote().insert(make_id({2, 8}, 0), 3, 0.0);
+  am3.clean_boundary_history(make_not_null(&history));
+  CHECK(substep_counts(history.local()) == Counts{{{1, 8}, 1}, {{2, 8}, 1}});
+  CHECK(substep_counts(history.remote()) == Counts{{{1, 8}, 1}, {{2, 8}, 1}});
+  history.local().insert(make_id({2, 8}, 1, {1, 8}), 3, 0.0);
+  history.remote().insert(make_id({2, 8}, 1, {2, 8}), 3, 0.0);
+  am3.clean_boundary_history(make_not_null(&history));
+  CHECK(substep_counts(history.local()) == Counts{{{1, 8}, 1}, {{2, 8}, 2}});
+  CHECK(substep_counts(history.remote()) == Counts{{{1, 8}, 1}, {{2, 8}, 2}});
+  history.local().insert(make_id({3, 8}, 0), 3, 0.0);
+  am3.clean_boundary_history(make_not_null(&history));
+  CHECK(substep_counts(history.local()) ==
+        Counts{{{1, 8}, 1}, {{2, 8}, 2}, {{3, 8}, 1}});
+  CHECK(substep_counts(history.remote()) == Counts{{{1, 8}, 1}, {{2, 8}, 2}});
+  history.local().insert(make_id({3, 8}, 1, {1, 8}), 3, 0.0);
+  am3.clean_boundary_history(make_not_null(&history));
+  CHECK(substep_counts(history.local()) == Counts{{{3, 8}, 1}});
+  CHECK(substep_counts(history.remote()) == Counts{{{2, 8}, 1}});
+
+  // Local coarser
+  history.local().insert(make_id({4, 8}, 0), 3, 0.0);
+  history.remote().insert(make_id({4, 8}, 0), 3, 0.0);
+  am3.clean_boundary_history(make_not_null(&history));
+  CHECK(substep_counts(history.local()) == Counts{{{3, 8}, 1}, {{4, 8}, 1}});
+  CHECK(substep_counts(history.remote()) == Counts{{{2, 8}, 1}, {{4, 8}, 1}});
+  history.local().insert(make_id({4, 8}, 1, {2, 8}), 3, 0.0);
+  history.remote().insert(make_id({4, 8}, 1, {1, 8}), 3, 0.0);
+  history.remote().insert(make_id({5, 8}, 0), 3, 0.0);
+  history.remote().insert(make_id({5, 8}, 1, {1, 8}), 3, 0.0);
+  am3.clean_boundary_history(make_not_null(&history));
+  CHECK(substep_counts(history.local()) == Counts{{{4, 8}, 1}});
+  CHECK(substep_counts(history.remote()) == Counts{{{5, 8}, 1}});
+}
+
+// [[Timeout, 30]]
 SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.AdamsMoultonPc", "[Unit][Time]") {
   test_am<false>();
+  test_boundary<false>();
+  test_non_monotonic_boundary_cleaning();
 }
 
-// [[Timeout, 10]]
+// [[Timeout, 30]]
 SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.AdamsMoultonPcMonotonic",
                   "[Unit][Time]") {
   test_am<true>();
+  test_boundary<true>();
+  test_monotonic_boundary_cleaning();
 }
 }  // namespace
diff --git a/tests/Unit/Time/TimeSteppers/Test_Factory.cpp b/tests/Unit/Time/TimeSteppers/Test_Factory.cpp
index 073b9889cc03..63e012da2a29 100644
--- a/tests/Unit/Time/TimeSteppers/Test_Factory.cpp
+++ b/tests/Unit/Time/TimeSteppers/Test_Factory.cpp
@@ -82,5 +82,11 @@ SPECTRE_TEST_CASE("Unit.Time.TimeSteppers.Factory", "[Unit][Time]") {
         }
         return true;
       });
+  check_list_contents<TimeSteppers::monotonic_lts_time_steppers>(
+      "monotonic_lts_time_steppers",
+      []<typename Stepper>(tmpl::type_<Stepper> /*meta*/) {
+        return std::is_convertible_v<Stepper*, LtsTimeStepper*> and
+               create_example(tmpl::type_<Stepper>{}).monotonic();
+      });
 }
 }  // namespace