diff --git a/src/include/serif/operators/concepts.hpp b/src/include/serif/operators/concepts.hpp
new file mode 100644
index 00000000..a8dba3a3
--- /dev/null
+++ b/src/include/serif/operators/concepts.hpp
@@ -0,0 +1,100 @@
+#pragma once
+
+#include <concepts>
+#include <type_traits>
+
+#include "serif/discretization/blocks/traits.hpp"
+#include "serif/discretization/functional_space/space.hpp"
+#include "serif/equations/concepts.hpp"
+
+namespace serif::operators {
+    // These restrictions describe the existing volume implementation,
+    // not which physical equations may be declared in serif::equations.
+    template <typename TermT>
+    concept ScalarValueTerm =
+        discretization::blocks::GlobalScalarTerm<TermT> ||
+        (discretization::blocks::SpatialTerm<TermT> &&
+        discretization::blocks::SpatialTermTraits<TermT>::rank == 0 &&
+        (
+            discretization::blocks::SpatialTermTraits<TermT>::family == discretization::functional_space::Family::h1 ||
+            discretization::blocks::SpatialTermTraits<TermT>::family == discretization::functional_space::Family::l2
+        ));
+
+    template <typename OperandT>
+    concept ScalarValueOperand =
+        discretization::forms::IsOperand<OperandT> &&
+        ScalarValueTerm<typename OperandT::term_type> &&
+        std::same_as<typename OperandT::operation_type, discretization::forms::Operation::Value>;
+
+    /// A conforming scalar gradient; discontinuous fields also require face terms.
+    template <typename OperandT>
+    concept ScalarGradientOperand =
+        discretization::forms::IsOperand<OperandT> &&
+        discretization::blocks::SpatialTerm<typename OperandT::term_type> &&
+        discretization::blocks::SpatialTermTraits<typename OperandT::term_type>::rank == 0 &&
+        discretization::blocks::SpatialTermTraits<typename OperandT::term_type>::family == discretization::functional_space::Family::h1 &&
+        std::same_as<typename OperandT::operation_type, discretization::forms::Operation::Gradient>;
+
+    /** @brief Full vector gradients are supported for componentwise H1 fields. */
+    template <typename OperandT>
+    concept VectorGradientOperand =
+        discretization::forms::IsOperand<OperandT> &&
+        discretization::blocks::SpatialTerm<typename OperandT::term_type> &&
+        discretization::blocks::SpatialTermTraits<typename OperandT::term_type>::rank == 1 &&
+        discretization::blocks::SpatialTermTraits<typename OperandT::term_type>::family == discretization::functional_space::Family::h1 &&
+        std::same_as<typename OperandT::operation_type, discretization::forms::Operation::Gradient>;
+
+    template <typename EquationT, typename OperandsT>
+    struct HasScalarVolumeOperands : std::false_type {};
+
+    template <typename EquationT, typename... Operands>
+    struct HasScalarVolumeOperands<EquationT, utils::types::TypeList<Operands...>> :
+        std::bool_constant<
+            sizeof...(Operands) != 0 &&
+            (ScalarValueOperand<Operands> && ...) &&
+            (requires(const EquationT& equation, Operands operand) {
+                {equation.additional_order(operand, 1)}->std::same_as<int>;
+            } && ...)> {};
+
+    /// Equations executable by the current scalar-value volume operator.
+    template <typename EquationT>
+    concept ScalarVolumeEquation =
+        equations::Equation<EquationT> &&
+        ScalarValueTerm<typename EquationT::Test> &&
+        HasScalarVolumeOperands<EquationT, typename EquationT::Operands>::value;
+    /// Mathematical operand combinations implemented by the current volume backend.
+    template <typename OperandT>
+    concept MappedOperand =
+        discretization::forms::IsOperand<OperandT> &&
+        (
+            ScalarValueOperand<OperandT> || ScalarGradientOperand<OperandT> || VectorGradientOperand<OperandT> ||
+            (
+                discretization::blocks::SpatialTerm<typename OperandT::term_type> &&
+                discretization::blocks::SpatialTermTraits<typename OperandT::term_type>::rank == 1 &&
+                discretization::blocks::SpatialTermTraits<typename OperandT::term_type>::family == discretization::functional_space::Family::h1 &&
+                (
+                    std::same_as<typename OperandT::operation_type, discretization::forms::Operation::Value> ||
+                    std::same_as<typename OperandT::operation_type, discretization::forms::Operation::Divergence>
+                )
+            )
+        );
+
+    template <typename C, typename ListT> struct VolumeOperands : std::false_type {};
+    template <typename C, typename... O>
+    struct VolumeOperands<C, utils::types::TypeList<O...>> : std::bool_constant<
+        sizeof...(O) != 0 && (MappedOperand<O> && ...) &&
+        (requires(const C& contribution, O operand) {
+        {contribution.additional_order(operand, 1)}->std::same_as<int>; } && ...)> {};
+
+    template <typename TupleT> struct VolumeContributions : std::false_type {};
+    template <typename... C>
+    struct VolumeContributions<std::tuple<C...>>
+        : std::bool_constant<((discretization::domain::IsDomainOrSet<equations::SupportOf<C>> &&
+                               MappedOperand<equations::TestOperandOf<C>> && VolumeOperands<C, typename C::Operands>::value) &&
+                              ...)> {};
+
+    /// One or more supported weak volume contributions assembled into a single residual row.
+    template <typename T>
+    concept VolumeEquation = equations::Equation<T> && VolumeContributions<equations::ContributionsOf<T>>::value;
+
+} // namespace serif::operators
diff --git a/src/include/serif/operators/contributions.hpp b/src/include/serif/operators/contributions.hpp
new file mode 100644
index 00000000..e35b7328
--- /dev/null
+++ b/src/include/serif/operators/contributions.hpp
@@ -0,0 +1,85 @@
+#pragma once
+
+#include <cstddef>
+#include <memory>
+#include <type_traits>
+#include <utility>
+#include <tuple>
+
+#include "serif/equations/concepts.hpp"
+
+namespace serif::operators {
+    /// Retain one leaf equation and the constant coefficient of its arithmetic expression.
+    template <typename EquationT>
+    struct ContributionBinding {
+        std::shared_ptr<const EquationT> equation;
+        double coefficient;
+    };
+
+    namespace detail {
+        /// A returned tuple may own generated equations or borrow an equation's existing children.
+        /// Retaining its parent keeps reference tuples valid without copying move-only equations.
+        template <typename EquationT>
+        struct ContributionChildren {
+            std::shared_ptr<const EquationT> parent;
+            decltype(std::declval<const EquationT&>().contributions()) values;
+
+            /// Evaluate contributions once, retaining either its returned value or its returned reference.
+            explicit ContributionChildren(std::shared_ptr<const EquationT> equation) :
+                parent(std::move(equation)),
+                values(parent->contributions()) {
+            }
+        };
+
+        /// Flatten equation arithmetic while preserving every generated tuple and borrowed parent.
+        template <typename EquationT>
+        [[nodiscard]] auto bind_contributions(std::shared_ptr<const EquationT> equation, double coefficient) {
+            if constexpr (requires { equation->contributions(); }) {
+                const auto children = std::make_shared<const ContributionChildren<EquationT>>(std::move(equation));
+                return std::apply([&](const auto&... child) {
+                    // Aliasing pointers address each child but retain the tuple and its parent together.
+                    return std::tuple_cat(bind_contributions(
+                        std::shared_ptr<const std::remove_cvref_t<decltype(child)>>(children, &child), coefficient
+                    )...);
+                }, children->values);
+            } else {
+                return std::tuple{ContributionBinding<EquationT>{
+                    .equation = std::move(equation),
+                    .coefficient = coefficient
+                }};
+            }
+        }
+
+        /// Visit arithmetic leaves without retaining them beyond the supplied callback.
+        template <typename EquationT, typename FunctionT>
+        void visit_contributions(
+            const EquationT& equation,
+            FunctionT& function,
+            std::size_t& index,
+            double coefficient
+        ) {
+            if constexpr (requires { equation.contributions(); }) {
+                std::apply([&](const auto&... child) {(visit_contributions(child, function, index, coefficient), ...);},equation.contributions());
+            } else {
+                function(equation, index++, coefficient);
+            }
+        }
+    } // namespace detail
+
+    /// Give each retained leaf a compile-time tuple position for independently cached contribution results.
+    template <equations::Equation EquationT>
+    [[nodiscard]] auto bind_contributions(std::shared_ptr<const EquationT> equation) {
+        return detail::bind_contributions(std::move(equation), 1.0);
+    }
+
+    /// Visit leaves in declaration order, multiplying constant coefficients along each path.
+    /// The callback receives (leaf equation, unique contribution index, coefficient).
+    template <equations::Equation EquationT, typename FunctionT>
+    void for_each_contribution(
+        const EquationT& equation,
+        FunctionT&& function
+    ) {
+        std::size_t index = 0;
+        detail::visit_contributions(equation, function, index, 1.0);
+    }
+} // namespace serif::operators
diff --git a/src/include/serif/operators/linearized_block.hpp b/src/include/serif/operators/linearized_block.hpp
new file mode 100644
index 00000000..60428d36
--- /dev/null
+++ b/src/include/serif/operators/linearized_block.hpp
@@ -0,0 +1,8 @@
+#pragma once
+
+#include "serif/backend/mfem/operators/linearized_block.hpp"
+
+namespace serif::operators {
+    using backend::mfem::LinearizedBlock;
+    using backend::mfem::linearized_block;
+}
diff --git a/src/include/serif/operators/linearized_system_view.hpp b/src/include/serif/operators/linearized_system_view.hpp
new file mode 100644
index 00000000..d8a4a42b
--- /dev/null
+++ b/src/include/serif/operators/linearized_system_view.hpp
@@ -0,0 +1,64 @@
+#pragma once
+
+#include <type_traits>
+
+#include "serif/discretization/blocks/system_layout.hpp"
+#include "serif/utils/dependencies/dependency_context.hpp"
+
+namespace serif::operators
+{
+    /** @brief Identify a physical Jacobian block by its equation rows and unknown coordinates. */
+    template <typename RowsT, typename ValuesT>
+    struct PhysicalBlock final
+    {
+        using Rows = RowsT;
+        using Values = ValuesT;
+    };
+
+    /** @brief Borrow an accepted collection of physical blocks and its independently ordered system layout. */
+    template <discretization::blocks::SystemLayoutType Layout, utils::dependencies::IsDependencyContext Blocks>
+    class LinearizedSystemView final
+    {
+    public:
+        using LayoutType = std::remove_cvref_t<Layout>;
+        using BlockSources = std::remove_cvref_t<Blocks>;
+
+        /** @brief Both owners must outlive the view and remain at this accepted linearization while it is used. */
+        LinearizedSystemView(const LayoutType& layout, const BlockSources& blocks) : m_layout(layout), m_blocks(blocks)
+        {
+        }
+
+        /** @brief A borrowed system cannot retain temporary coordinate metadata. */
+        LinearizedSystemView(const LayoutType&&, const BlockSources&) = delete;
+        /** @brief A borrowed system cannot retain temporary block sources. */
+        LinearizedSystemView(const LayoutType&, const BlockSources&&) = delete;
+
+        /** @brief Access a declared physical coupling; absent blocks are compile-time errors. */
+        template <typename Rows, typename Values>
+            requires utils::dependencies::BindingsWithBinding<typename BlockSources::Bindings, PhysicalBlock<Rows, Values>>
+        [[nodiscard]] const auto& block() const
+        {
+            return m_blocks.template get<PhysicalBlock<Rows, Values>>();
+        }
+
+        /** @brief Access the shared coordinate description without taking ownership. */
+        [[nodiscard]] const LayoutType& layout() const noexcept
+        {
+            return m_layout;
+        }
+
+        /** @brief Supply accepted block sources to ordinary dependency preparation. */
+        [[nodiscard]] const BlockSources& blocks() const noexcept
+        {
+            return m_blocks;
+        }
+
+    private:
+        const LayoutType& m_layout;
+        const BlockSources& m_blocks;
+    };
+
+    /** @brief Infer canonical owner types while retaining the explicit borrowing lifetime checks. */
+    template <discretization::blocks::SystemLayoutType Layout, utils::dependencies::IsDependencyContext Blocks>
+    LinearizedSystemView(const Layout&, const Blocks&) -> LinearizedSystemView<Layout, Blocks>;
+} // namespace serif::operators
diff --git a/src/include/serif/operators/plan.hpp b/src/include/serif/operators/plan.hpp
new file mode 100644
index 00000000..4ff219f0
--- /dev/null
+++ b/src/include/serif/operators/plan.hpp
@@ -0,0 +1,85 @@
+#pragma once
+
+#include <tuple>
+#include "serif/utils/dependencies/key_list_operations.hpp"
+#include "serif/discretization/forms/forms.hpp"
+#include "serif/operators/concepts.hpp"
+
+namespace serif::operators {
+    namespace detail {
+        template <typename ListT> struct AsKeyList;
+        template <typename... T> struct AsKeyList<utils::types::TypeList<T...>> {
+            using type = utils::dependencies::CombinedKeyList<utils::dependencies::KeyList<T>...>;
+        };
+        template <typename... Lists> struct Join {
+            using type = typename utils::dependencies::CombinedKeyList<typename AsKeyList<Lists>::type...>::Keys;
+        };
+        template <typename TupleT> struct ContributionOperands;
+        template <typename... C> struct ContributionOperands<std::tuple<C...>> : Join<typename C::Operands...> {};
+        template <typename ListT> struct Columns;
+        template <typename... O>
+        struct Columns<utils::types::TypeList<O...>> {
+            using type = typename utils::dependencies::CombinedKeyList<utils::dependencies::KeyList<discretization::forms::Operand<typename O::term_type>>...>::Keys;
+        };
+
+        template <typename Input, typename TestOperand> struct CouplingForm {
+            using Trial = typename Input::term_type;
+            using Test = typename TestOperand::term_type;
+            using type = std::conditional_t<
+                std::same_as<typename Input::operation_type, discretization::forms::Operation::Divergence>,
+                discretization::forms::divergence<Trial, Test>,
+                std::conditional_t<std::same_as<typename TestOperand::operation_type, discretization::forms::Operation::Divergence>,
+                                   discretization::forms::divergence<Test, Trial>,
+                                   std::conditional_t<std::same_as<Trial, Test>, discretization::forms::mass<Trial, Test>,
+                                                      discretization::forms::source<Trial, Test>>>>;
+        };
+    } // namespace detail
+
+    template <VolumeEquation EquationT,
+              typename OperandsT = typename detail::ContributionOperands<equations::ContributionsOf<EquationT>>::type>
+    struct EquationPlan;
+
+    template <VolumeEquation EquationT, typename... O> struct EquationPlan<EquationT, utils::types::TypeList<O...>> {
+        using Test = typename EquationT::Test;
+        using Support = discretization::blocks::SupportOf<Test>;
+        using Operands = utils::types::TypeList<O...>;
+        using Columns = typename detail::Columns<Operands>::type;
+        static constexpr std::size_t operand_count = sizeof...(O);
+        static constexpr std::size_t column_count = equations::OperandValues<Columns>::count;
+        static constexpr int component_count = (discretization::forms::component_count<discretization::forms::OperandValue<O>> + ... + 0);
+
+        template <typename OperandT, typename ContributionT = EquationT>
+        using Form = typename detail::CouplingForm<OperandT, equations::TestOperandOf<ContributionT>>::type;
+
+        /// Locate the unique coefficient column shared by every operation on a term.
+        template <typename OperandT> static consteval std::size_t column_index() {
+            return equations::OperandValues<Columns>::template index<typename OperandT::term_type>();
+        }
+        /// Locate the first scalar component of an operand in pointwise derivative storage.
+        template <typename OperandT> static consteval int component_offset() {
+            constexpr std::array counts{discretization::forms::component_count<discretization::forms::OperandValue<O>>...};
+            constexpr auto index = equations::OperandValues<Operands>::template index<OperandT>();
+            int offset = 0;
+            for (std::size_t i = 0; i < index; ++i) {
+                offset += counts[i];
+            }
+            return offset;
+        }
+        /// Visit all sampled operations in stable first-appearance order.
+        template <typename FunctionT> static void for_each_operand(FunctionT&& function) {
+            (function(O{}, equations::OperandValues<Operands>::template index<O>()), ...);
+        }
+        /// Visit only the operands needed by this contribution, retaining global sample indices.
+        template <typename ContributionT, typename FunctionT> static void for_each_input(FunctionT&& function) {
+            [&]<typename... Input>(utils::types::TypeList<Input...>) {
+                (function(Input{}, equations::OperandValues<Operands>::template index<Input>()), ...);
+            }(typename ContributionT::Operands{});
+        }
+        /// Visit each physical coefficient vector once for binding and packed actions.
+        template <typename FunctionT> static void for_each_column(FunctionT&& function) {
+            [&]<typename... Column>(utils::types::TypeList<Column...>) {
+                (function(Column{}, equations::OperandValues<Columns>::template index<Column>()), ...);
+            }(Columns{});
+        }
+    };
+} // namespace serif::operators