diff --git a/utilities/CMakeLists.txt b/utilities/CMakeLists.txt
index 2dd0cd05fb..76b6aad338 100644
--- a/utilities/CMakeLists.txt
+++ b/utilities/CMakeLists.txt
@@ -12,3 +12,4 @@ foreach(current IN LISTS EXECUTABLE_SOURCES)
 endforeach()
 
 add_subdirectory(external-interface)
+add_subdirectory(cost_analysis)
diff --git a/utilities/cost_analysis/CMakeLists.txt b/utilities/cost_analysis/CMakeLists.txt
new file mode 100644
index 0000000000..131767ac7c
--- /dev/null
+++ b/utilities/cost_analysis/CMakeLists.txt
@@ -0,0 +1,35 @@
+include(FindOrFetchCLI11)
+include(FindOrFetchJSON)
+
+add_executable(cost_analysis cost_analysis.cpp report.cpp)
+set_target_properties(cost_analysis PROPERTIES
+    CXX_SCAN_FOR_MODULES OFF
+    RUNTIME_OUTPUT_DIRECTORY "${PROJECT_BINARY_DIR}")
+target_link_libraries(cost_analysis PRIVATE
+    SeQuant::mbpt nlohmann_json::nlohmann_json CLI11::CLI11)
+target_set_warning_flags(cost_analysis)
+
+foreach(NAME IN ITEMS "ccsd_r2" "df_r1")
+    # --omit-revision drops the volatile git-revision line so a plain
+    # compare_files suffices; the tool writes the report next to the driver.
+    add_test(
+        NAME "sequant/cost_analysis/${NAME}/generate"
+        COMMAND cost_analysis --driver "${CMAKE_CURRENT_LIST_DIR}/examples/${NAME}.json" --omit-revision)
+    set_tests_properties("sequant/cost_analysis/${NAME}/generate"
+        PROPERTIES FIXTURES_SETUP "COST_ANALYSIS_${NAME}")
+
+    add_test(
+        NAME "sequant/cost_analysis/${NAME}/verify"
+        COMMAND "${CMAKE_COMMAND}" -E compare_files "${NAME}.md" "${NAME}.md.expected"
+        WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/examples")
+    set_tests_properties("sequant/cost_analysis/${NAME}/verify"
+        PROPERTIES FIXTURES_REQUIRED "COST_ANALYSIS_${NAME}")
+endforeach()
+
+# ccsd_r2 exercises --dump_tree; diff the tree dump against its fixture.
+add_test(
+    NAME "sequant/cost_analysis/ccsd_r2/dump_tree"
+    COMMAND "${CMAKE_COMMAND}" -E compare_files "R2.tree.txt" "R2.tree.txt.expected"
+    WORKING_DIRECTORY "${CMAKE_CURRENT_LIST_DIR}/examples")
+set_tests_properties("sequant/cost_analysis/ccsd_r2/dump_tree"
+    PROPERTIES FIXTURES_REQUIRED "COST_ANALYSIS_ccsd_r2")
diff --git a/utilities/cost_analysis/README.md b/utilities/cost_analysis/README.md
new file mode 100644
index 0000000000..7c495f78f6
--- /dev/null
+++ b/utilities/cost_analysis/README.md
@@ -0,0 +1,118 @@
+# Cost analysis
+
+A standalone utility that estimates the cost of evaluating a tensor equation
+without running it. It reads a serialized equation, optimizes and binarizes it,
+then writes a Markdown report on the largest intermediates, the most expensive
+contractions, peak storage, total operation count, and (optionally) cache reuse.
+Figures are symbolic `AsyCost` polynomials in the index-space sizes; memory and
+storage are additionally evaluated to megabytes at the driver's sizes (element
+width follows the field: 8 bytes real, 16 bytes complex). Operation counts are
+reported symbolically only.
+
+The optimize/binarize step mirrors the choices MPQC's `SeQuantEngine` makes
+(symmetrizer stripping, per-summand factorization, volatile-leaf weighting,
+cache gating), so the estimates track what a real MPQC evaluation would do.
+
+## Build & run
+
+Built as part of SeQuant (target `cost_analysis`). It takes a single argument, a
+[JSON](https://www.json.org/) _driver_:
+
+```bash
+cost_analysis --driver examples/ccsd_r2.json
+```
+
+Paths inside the driver (`equation_file`, `output.path`) — and the `dump_tree`
+`<name>.tree.txt` files — are resolved/written relative to the driver's own
+directory, not the invocation directory. The final "Report written to …" line
+prints the absolute path so the output is easy to locate.
+
+## Driver
+
+| Block | Key | Meaning | Default |
+|---|---|---|---|
+| `context` | `spbasis` | `spinor` or `spinfree` | `spinor` |
+| | `field` | `real` or `complex` | `real` |
+| | `convention` | registry: `min_sr`/`sr`/`mr`/`f12` | `sr` |
+| | `aux` | factorization spaces to register: `["df"]` (Κ) and/or `["thc"]` (L) | none |
+| `sizes` | `<label>` | approximate size per index space, e.g. `{"i": 10, "a": 38}` | — |
+| `optimize` | `objective` | `dense_flops` / `dense_size` / `dense_peak_size` | `dense_flops` |
+| | `volatile_leaf` | amplitude label being solved (e.g. `t`, `R`); empty disables volatile weighting | `R` |
+| | `reorder`, `cse_subnet`, `volatile_weight`, `machine_balance`, `fast_mem_elems` | forwarded to `sequant::OptimizeOptions` | see `OptSpec` |
+| `cache` | `enabled` | run the cache-reuse simulation | `true` |
+| | `min_repeats` | reuse threshold to cache an intermediate | `2` |
+| | `max_footprint` | cache budget, in **elements** (not MB); `0` = unlimited | `0` |
+| `output` | `path` | report filename | `cost_analysis.md` |
+| | `top_n` | rows in the largest/expensive tables | `20` |
+| | `dump_tree` | also write each result's binarized tree to `<name>.tree.txt` | `false` |
+| `equations` | `[{name, equation_file}]` | equations to analyze | — |
+
+> [!IMPORTANT]
+> Every index space that appears in an equation should be given a size. Any
+> space omitted from `sizes` falls back to its standard-registry default, so
+> the reported MB/operation figures stay well-defined but reflect that default
+> rather than your problem — size all used spaces (including any `aux` space) to
+> get meaningful numbers.
+
+Each `equation_file` holds one `<head> = <rhs>` in SeQuant serialization V1, e.g.
+
+```
+R{a1,a2;i1,i2} = g{a1,a2;i1,i2} + g{a3,a4;i1,i2} t{i3,i4;a3,a4} R{a1,a2;i3,i4}
+```
+
+> [!NOTE]
+> CSV/PNO expressions are not supported: `sequant::AsyCost` is not proto-index aware, so
+> it would size such an index as a bare one and report wrong costs.
+
+## Workflow
+
+Per result equation (`process()` in `cost_analysis.cpp`):
+
+1. **Deserialize** the `<head> = <rhs>` string into a `ResultExpr`.
+2. **Strip** a leading (anti)symmetrizer — it is re-expressed by the head, not
+   evaluated.
+3. **Flatten** nested sums/products (mirrors `SeQuantEngine`; a no-op for
+   already-flat serialized input).
+4. Per summand: **`optimize`** (choose a binary contraction order for the
+   objective) then **`binarize`** into an evaluation tree — *once*; the tree is
+   the single source of truth for every reported number.
+5. **Catalog** the tree's internal nodes, keyed by structural hash/equality so
+   equal intermediates across terms collapse into one entry (with a `uses`
+   count). Record each node's memory, local operation count, and O/V/X space
+   signature.
+
+A final **cache simulation** (`cache_manager`) runs over all results' trees to
+count intermediates that would be cached / persist across iterations.
+
+## Report
+
+One `##` section per result — a one-line summary and the **Largest
+intermediates**, **Most expensive contractions**, and **Shape census** tables —
+followed by a single **Cache** section pooled across all results. Memory is the
+tensor's element count (product of its index-space sizes) times the field's
+element width (8 B real / 16 B complex), in MB; operation counts stay symbolic.
+
+The quantities that aren't self-evident:
+- **distinct intermediates** — structurally-unique internal contraction nodes;
+  equal intermediates across terms collapse to one, scalars excluded.
+- **peak storage** — the transient working set of the *heaviest single
+  contraction* (its two operands plus result, via `min_storage`), maxed over the
+  terms. **Not** the cache footprint — the **Cache** table covers that.
+- **Shape census `Size`** — the largest instance of each O/V/X signature.
+- **Cache** — a simulation over the volatile-gated `cache_manager` (as MPQC
+  builds it), pooled across all results. **Cached** / **Persistent** count
+  intermediates reused ≥ `min_repeats` (volatile-leaf gated) / those that persist
+  across iterations; **footprints** are their summed element counts in MB — a
+  cache-content total, distinct from the per-result transient **peak storage**.
+
+## Tests
+
+Reference tests (`ctest -R sequant/cost_analysis`) run the tool on `examples/`
+with `--omit-revision` (which drops the volatile git-revision line, making the
+report byte-stable) and compare it against a frozen `*.md.expected` via
+`cmake -E compare_files`:
+
+- `.../ccsd_r2/*` — spin-orbital CCSD R2 (a `Sum` of terms) whose two ladder
+  terms share one `g*t` intermediate, exercising the reuse census; also diffs
+  the `--dump_tree` output against `R2.tree.txt.expected`.
+- `.../df_r1/*` — a single density-fitted product (the non-`Sum` path).
diff --git a/utilities/cost_analysis/cost_analysis.cpp b/utilities/cost_analysis/cost_analysis.cpp
new file mode 100644
index 0000000000..57a616fb99
--- /dev/null
+++ b/utilities/cost_analysis/cost_analysis.cpp
@@ -0,0 +1,392 @@
+// SeQuant expression cost/factorization analysis.
+//
+// Reads a JSON driver (space sizes + equation files), reads each equation in
+// as-is, then optimizes -> binarizes it and writes a Markdown report on the
+// largest intermediates, shape census, and (gated) cache footprint. Shared
+// data types live in cost_analysis.hpp; report formatting lives in
+// report.{hpp,cpp}.
+
+#include "cost_analysis.hpp"
+#include "report.hpp"
+
+#include <CLI/CLI.hpp>
+#include <nlohmann/json.hpp>
+
+#include <SeQuant/core/context.hpp>
+#include <SeQuant/core/eval/cache_manager.hpp>
+#include <SeQuant/core/index_space_registry.hpp>
+#include <SeQuant/core/io/shorthands.hpp>  // deserialize
+#include <SeQuant/core/optimize/optimize.hpp>
+#include <SeQuant/core/reserved.hpp>
+#include <SeQuant/core/runtime.hpp>  // set_locale
+#include <SeQuant/core/tensor_canonicalizer.hpp>
+#include <SeQuant/core/utility/expr.hpp>
+#include <SeQuant/core/utility/string.hpp>     // toUtf16
+#include <SeQuant/domain/mbpt/convention.hpp>  // make_min_sr_spaces
+
+#include <filesystem>
+#include <fstream>
+#include <iostream>
+#include <memory>
+#include <stdexcept>
+#include <string>
+#include <utility>
+#include <vector>
+
+using nlohmann::json;
+using namespace sequant;
+using namespace sequant::mbpt;
+
+namespace {
+
+Config load_config(const json& d) {
+  Config c;
+
+  if (d.contains("context")) {
+    const auto& ctx = d.at("context");
+    c.spinor = ctx.value("spbasis", std::string("spinor")) == "spinor";
+    const auto field = ctx.value("field", std::string("real"));
+    if (field != "real" && field != "complex")
+      throw std::runtime_error("unknown field (want real|complex): " + field);
+    c.real_field = field == "real";
+    c.convention = ctx.value("convention", std::string("sr"));
+    for (const auto& a : ctx.value("aux", json::array())) {
+      const auto s = a.get<std::string>();
+      if (s != "df" && s != "thc")
+        throw std::runtime_error("unknown aux (want df|thc): " + s);
+      c.aux.push_back(s);
+    }
+  }
+
+  for (const auto& [label, size] : d.at("sizes").items()) {
+    const auto sz = size.get<long long>();
+    if (sz <= 0)
+      throw std::runtime_error("index-space size must be > 0: " + label);
+    c.sizes[label] = static_cast<std::size_t>(sz);
+  }
+
+  if (d.contains("optimize")) {
+    const auto& o = d.at("optimize");
+    c.optimize.objective = o.value("objective", std::string("dense_flops"));
+    c.optimize.reorder = o.value("reorder", true);
+    c.optimize.cse_subnet = o.value("cse_subnet", false);
+    c.optimize.volatile_leaf = o.value("volatile_leaf", std::string("R"));
+    c.optimize.volatile_weight = o.value("volatile_weight", 10.0);
+    c.optimize.machine_balance = o.value("machine_balance", 0.0);
+    c.optimize.fast_mem_elems = o.value("fast_mem_elems", 0.0);
+  }
+  if (d.contains("cache")) {
+    const auto& ca = d.at("cache");
+    c.cache.enabled = ca.value("enabled", true);
+    // A negative min_repeats wraps to a huge size_t (JSON ints are signed) and
+    // silently disables caching, so reject it.
+    const long mr = ca.value("min_repeats", 2L);
+    if (mr < 0) throw std::runtime_error("cache.min_repeats must be >= 0");
+    c.cache.min_repeats = static_cast<std::size_t>(mr);
+    c.cache.max_footprint = ca.value("max_footprint", 0.0);
+    if (c.cache.max_footprint < 0.0)
+      throw std::runtime_error("cache.max_footprint must be >= 0");
+  }
+  if (d.contains("output")) {
+    const auto& ou = d.at("output");
+    c.out.path = ou.value("path", std::string("cost_analysis.md"));
+    // Read signed then guard: a negative top_n would wrap to a huge size_t.
+    const long tn = ou.value("top_n", 20L);
+    if (tn < 0) throw std::runtime_error("output.top_n must be >= 0");
+    c.out.top_n = static_cast<std::size_t>(tn);
+    c.out.dump_tree = ou.value("dump_tree", false);
+  }
+
+  for (const auto& r : d.at("equations")) {
+    EquationSpec rs;
+    rs.name = r.at("name").get<std::string>();
+    rs.equation_file = r.at("equation_file").get<std::string>();
+    c.equations.push_back(std::move(rs));
+  }
+  return c;
+}
+
+std::shared_ptr<IndexSpaceRegistry> make_registry(const std::string& conv) {
+  if (conv == "min_sr") return make_min_sr_spaces();
+  if (conv == "sr") return make_sr_spaces();
+  if (conv == "mr") return make_mr_spaces();
+  if (conv == "f12") return make_F12_sr_spaces();
+  throw std::runtime_error("unknown convention (want min_sr|sr|mr|f12): " +
+                           conv);
+}
+
+// Start from one of SeQuant's standard registries then override each named
+// space's approximate size and field from the config.
+void setup_context(const Config& cfg) {
+  auto isr = make_registry(cfg.convention);
+  // Register auxiliary factorization spaces (Κ for df, L for thc) so
+  // pre-factorized DF/THC equations parse; validated to df|thc in load_config.
+  for (const auto& a : cfg.aux) {
+    if (a == "df")
+      add_df_spaces(isr);
+    else
+      add_thc_spaces(isr);
+  }
+
+  // Field is a single global choice for the computation, so apply it to every
+  // registry space (not just the sized ones) the way MPQC's
+  // scoped_sequant_field does: bra<->ket symmetry must be consistent for all
+  // used spaces, else a used-but-unsized space keeps the registry's default
+  // (complex) and its tensors canonicalize wrongly over a real field. Collect
+  // keys first (iteration is const) then set via mutable lookup.
+  const Field field = cfg.real_field ? Field::Real : Field::Complex;
+  std::vector<std::wstring> keys;
+  for (const auto& sp : *isr) keys.push_back(sp.base_key());
+  for (const auto& key : keys)
+    if (IndexSpace* sp = isr->retrieve_ptr(key)) sp->field(field);
+
+  for (const auto& [label, size] : cfg.sizes) {
+    IndexSpace* sp = isr->retrieve_ptr(toUtf16(label));
+    if (!sp) throw std::runtime_error("unknown index space: " + label);
+    sp->approximate_size(size);
+  }
+
+  auto ctx = get_default_context();
+  auto copts = CanonicalizeOptions::default_options().copy_and_set(
+      CanonicalizationMethod::Complete);
+  ctx.set(Vacuum::SingleProduct)
+      .set(cfg.spinor ? SPBasis::Spinor : SPBasis::Spinfree)
+      .set(isr)
+      .set(BraKetSlotTypesetting::Naive)
+      .set(copts);
+  set_default_context(ctx);
+
+  TensorCanonicalizer::register_instance(
+      std::make_shared<DefaultTensorCanonicalizer>());
+}
+
+ResultExpr parse_equation(const std::string& text) {
+  return deserialize<ResultExpr>(std::string_view{text});
+}
+
+/// Cost Analysis
+
+ObjectiveFunction objective_of(const std::string& s) {
+  if (s == "dense_size") return ObjectiveFunction::DenseSize;
+  if (s == "dense_peak_size") return ObjectiveFunction::DensePeakSize;
+  if (s == "dense_flops") return ObjectiveFunction::DenseFLOPs;
+  throw std::runtime_error(
+      "unknown objective (want dense_size|dense_peak_size|dense_flops): " + s);
+}
+
+// Classify each index as Occupied / Virtual / other (aux) by matching its space
+// Type against the registry's designated hole/particle types.
+std::string space_signature(const EvalExpr& ev) {
+  if (!ev.is_tensor()) return "scalar";
+  const auto& t = ev.as_tensor();
+  const auto& reg = *get_default_context().index_space_registry();
+  const auto hole = reg.hole_space(/*nulltype_ok=*/true);
+  const auto particle = reg.particle_space(/*nulltype_ok=*/true);
+  std::string result;
+  for (const auto& idx : t.const_braketaux_indices()) {
+    const auto ty = idx.space().type();
+    if (ty == hole)
+      result += 'O';
+    else if (ty == particle)
+      result += 'V';
+    else
+      result += 'X';  // neither hole nor particle (e.g. df/thc aux)
+  }
+  return result;
+}
+
+AsyCost result_memory(const EvalExpr& ev) {
+  if (!ev.is_tensor()) return AsyCost::zero();
+  const auto& t = ev.as_tensor();
+  AsyCost::ExponentMap exponents;
+  // Include aux (braketaux, not braket) so auxiliary dimensions are counted in
+  // the intermediate's memory/cost, consistent with space_signature.
+  for (const auto& idx : t.const_braketaux_indices()) ++exponents[idx.space()];
+  return AsyCost{std::move(exponents)};
+}
+
+OptimizeOptions make_opts(const Config& cfg) {
+  OptimizeOptions o;
+  o.objective_function = objective_of(cfg.optimize.objective);
+  o.reorder =
+      cfg.optimize.reorder ? ReorderSum::Reorder : ReorderSum::NoReorder;
+  o.CSE.subnet = cfg.optimize.cse_subnet;
+  // volatile leaf (label match), MPQC-style volatile weighting.
+  if (!cfg.optimize.volatile_leaf.empty()) {
+    const std::wstring vl = toUtf16(cfg.optimize.volatile_leaf);
+    o.batch_policy.is_volatile_leaf = [vl](const Tensor& t) {
+      return t.label() == vl;
+    };
+    // MPQC: volatile_weight pinned to 1.0 when caching is off.
+    o.volatile_weight = cfg.cache.enabled ? cfg.optimize.volatile_weight : 1.0;
+    if (cfg.cache.enabled) {
+      o.roofline.machine_balance = cfg.optimize.machine_balance;
+      o.roofline.fast_mem_elems = cfg.optimize.fast_mem_elems;
+    }
+  }
+  return o;
+}
+
+CellResult process(const Config& cfg, const ExprPtr& rhs, const Tensor& head) {
+  CellResult res;
+
+  ExprPtr input = rhs->clone();
+  auto popped = pop_tensor(input, reserved::symm_label());
+  if (!popped.has_value()) pop_tensor(input, reserved::antisymm_label());
+
+  flatten(input);
+
+  const auto opts = make_opts(cfg);
+
+  auto scan = [&](const TreeNode& tree) {
+    tree.visit_internal([&](const TreeNode& node) {
+      if (node.root())
+        return;  // per-term root is the output, not an intermediate
+      auto it = res.catalog.find(node);
+      if (it == res.catalog.end()) {
+        Record rec;
+        rec.uses = 1;
+        rec.memory = result_memory(*node);
+        rec.flops = Flops{}(node);
+        rec.label = node->label();
+        rec.spaces = space_signature(*node);
+        res.catalog.emplace(node, std::move(rec));
+      } else {
+        ++it->second.uses;
+      }
+    });
+  };
+
+  auto handle = [&](const ExprPtr& summand) {
+    ExprPtr os = optimize(summand, opts);
+    auto tree = binarize<EvalExpr>(ResultExpr{head, os});
+    scan(tree);
+    // Peak working set = heaviest single-contraction over the schedule; report
+    // the max across the term-by-term summand trees. Total FLOPs = their sum.
+    const AsyCost ps = min_storage(tree);
+    if (res.peak_storage < ps) res.peak_storage = ps;
+    res.total_flops = res.total_flops + asy_cost(tree);
+    res.trees.push_back(std::move(tree));
+  };
+
+  if (input->is<Sum>()) {
+    const auto& sum = input->as<Sum>();
+    res.n_terms = sum.size();
+    for (std::size_t i = 0; i < sum.size(); ++i) handle(sum.summand(i));
+  } else {
+    res.n_terms = 1;
+    handle(input);
+  }
+
+  for (const auto& [_, rec] : res.catalog) {
+    if (rec.memory == AsyCost::zero()) continue;
+    ++res.n_distinct;
+    if (res.largest_mem < rec.memory) res.largest_mem = rec.memory;
+    if (rec.uses >= 2) ++res.n_reused;
+  }
+  return res;
+}
+
+/// Cache Simulation
+
+bool is_leaf_labeled(const TreeNode& n, const std::wstring& label) {
+  return n.leaf() && n->is_tensor() && n->as_tensor().label() == label;
+}
+
+SimResult simulate_cache(
+    const Config& cfg,
+    const std::vector<std::pair<std::string, CellResult>>& results) {
+  std::vector<TreeNode> forest;
+  for (const auto& [name, cell] : results)
+    for (const auto& tree : cell.trees) forest.push_back(tree);
+
+  const std::wstring vl = toUtf16(cfg.optimize.volatile_leaf);
+  auto is_volatile = [vl](const TreeNode& n) { return is_leaf_labeled(n, vl); };
+  auto cm = cache_manager(
+      forest, is_volatile, cfg.cache.min_repeats,
+      [](const TreeNode& n) { return result_memory(*n).ops(); },
+      cfg.cache.max_footprint);
+
+  SimResult r;
+  cm.for_each_key([&](const TreeNode& n) {
+    const AsyCost mem = result_memory(*n);
+    r.cached_footprint += mem;
+    ++r.n_cached;
+    if (cm.persistent(n)) {
+      r.persistent_footprint += mem;
+      ++r.n_persistent;
+    }
+  });
+  return r;
+}
+
+std::string read_file(const std::filesystem::path& p) {
+  std::ifstream in(p);
+  if (!in) throw std::runtime_error("cannot open equation file: " + p.string());
+  return std::string(std::istreambuf_iterator<char>(in), {});
+}
+
+}  // namespace
+
+int main(int argc, char** argv) {
+  set_locale();
+  sequant::detail::OpIdRegistrar op_id_registrar;
+
+  CLI::App app("SeQuant expression cost/factorization analysis");
+  argv = app.ensure_utf8(argv);
+  std::filesystem::path driver;
+  app.add_option("--driver", driver, "Path to the JSON driver file")
+      ->required();
+  bool omit_revision = false;
+  app.add_flag("--omit-revision", omit_revision,
+               "Omit the git-revision line (for reproducible test diffs)");
+  CLI11_PARSE(app, argc, argv);
+
+  try {
+    std::ifstream din(driver);
+    if (!din) {
+      std::cerr << "cannot open driver: " << driver << "\n";
+      return 1;
+    }
+    json d = json::parse(din, nullptr, true, /*skip_comments=*/true);
+    std::filesystem::current_path(
+        std::filesystem::absolute(driver).parent_path());
+
+    Config cfg = load_config(d);
+    cfg.out.omit_revision = omit_revision;
+    setup_context(cfg);  // registry baked before any parse
+
+    std::vector<std::pair<std::string, CellResult>> results;
+    for (const auto& r : cfg.equations) {
+      const ResultExpr res = parse_equation(read_file(r.equation_file));
+      CellResult cell = process(cfg, res.expression(), res.result_as_tensor());
+
+      if (cfg.out.dump_tree) {
+        const std::string dump_path = r.name + ".tree.txt";
+        std::ofstream to(dump_path);
+        if (!to)
+          throw std::runtime_error("cannot open dump file: " + dump_path);
+        for (const auto& tree : cell.trees) to << full_expr(tree) << "\n";
+      }
+      results.emplace_back(r.name, std::move(cell));
+    }
+
+    const SimResult sim =
+        cfg.cache.enabled ? simulate_cache(cfg, results) : SimResult{};
+
+    std::ofstream out(cfg.out.path);
+    if (!out)
+      throw std::runtime_error("cannot open output file: " + cfg.out.path);
+    write_report(cfg, results, sim, out);
+    // Print the absolute path: the working directory was switched to the
+    // driver's parent above, so a relative out.path lands there, not in the
+    // user's invocation directory.
+    std::cout << "Report written to " << std::filesystem::absolute(cfg.out.path)
+              << "\n";
+    return 0;
+
+  } catch (const std::exception& e) {
+    std::cerr << "error: " << e.what() << "\n";
+    return 1;
+  }
+}
diff --git a/utilities/cost_analysis/cost_analysis.hpp b/utilities/cost_analysis/cost_analysis.hpp
new file mode 100644
index 0000000000..91d5596f5a
--- /dev/null
+++ b/utilities/cost_analysis/cost_analysis.hpp
@@ -0,0 +1,99 @@
+#ifndef SEQUANT_UTILITIES_COST_ANALYSIS_COST_ANALYSIS_HPP
+#define SEQUANT_UTILITIES_COST_ANALYSIS_COST_ANALYSIS_HPP
+
+/// \file
+/// Data types shared between the pipeline (`cost_analysis.cpp`) and the
+/// reporter (`report.cpp`): the parsed JSON `Config` and the per-equation
+/// analysis results. Logic lives in the `.cpp` files.
+
+#include <SeQuant/core/asy_cost.hpp>
+#include <SeQuant/core/eval/eval_expr.hpp>
+#include <SeQuant/core/eval/eval_node.hpp>
+#include <SeQuant/core/eval/eval_node_compare.hpp>
+#include <SeQuant/core/expr.hpp>
+#include <SeQuant/core/expressions/result_expr.hpp>
+
+#include <cstddef>
+#include <map>
+#include <string>
+#include <unordered_map>
+#include <vector>
+
+/// `optimize` driver block; fields map onto `sequant::OptimizeOptions` (see
+/// there for their semantics). `volatile_leaf` is the amplitude label solved
+/// for (empty => no volatile weighting).
+struct OptSpec {
+  std::string objective = "dense_flops";
+  bool reorder = true, cse_subnet = false;
+  std::string volatile_leaf = "R";
+  double volatile_weight = 10.0, machine_balance = 0.0, fast_mem_elems = 0.0;
+};
+
+/// `cache` driver block: parameters of the gated cache-reuse simulation.
+struct CacheSpec {
+  bool enabled = true;
+  std::size_t min_repeats = 2;
+  double max_footprint = 0.0;
+};
+
+/// `output` driver block: report path and verbosity.
+struct OutSpec {
+  std::string path = "cost_analysis.md";
+  std::size_t top_n = 20;  ///< rows in the largest/expensive tables
+  bool dump_tree = false;  ///< also write each result's binarized tree
+  bool omit_revision =
+      false;  ///< drop the git-revision line (reproducible diffs)
+};
+
+/// One `equations[]` entry: a named equation file to analyze.
+struct EquationSpec {
+  std::string name, equation_file;
+};
+
+/// The parsed JSON driver.
+struct Config {
+  bool spinor = true;
+  bool real_field = true;  // true -> Field::Real (the utility's default); the
+                           // field is a single global choice for a computation
+  std::string convention = "sr";  // standard registry: min_sr|sr|mr|f12
+  std::vector<std::string>
+      aux;  // extra factorization spaces to register so
+            // DF/THC equations parse: df (Κ) and/or thc (L)
+  std::map<std::string, std::size_t> sizes;  // space label -> approximate size,
+                                             // overriding the standard registry
+  OptSpec optimize;
+  CacheSpec cache;
+  OutSpec out;
+  std::vector<EquationSpec> equations;
+};
+
+/// A binarized evaluation-tree node, keyed in the catalog by structural
+/// hash/equality so equal intermediates across terms collapse to one entry.
+using TreeNode = sequant::EvalNode<sequant::EvalExpr>;
+using Hasher = sequant::TreeNodeHasher<TreeNode>;
+using Comparator = sequant::TreeNodeEqualityComparator<TreeNode>;
+
+/// Catalog entry for one distinct intermediate (eval-tree node).
+struct Record {
+  std::size_t uses = 0;            ///< how many terms build this node
+  sequant::AsyCost memory, flops;  ///< result footprint and local op count
+  std::string label, spaces;       ///< tensor label and O/V/X signature
+};
+
+/// Cost analysis of one result equation.
+struct CellResult {
+  std::size_t n_terms = 0, n_distinct = 0, n_reused = 0;
+  sequant::AsyCost largest_mem, peak_storage, total_flops;
+  std::unordered_map<TreeNode, Record, Hasher, Comparator> catalog;
+  /// Binarized tree per summand, built once and reused by the catalog scan,
+  /// the cache simulation, and --dump_tree.
+  std::vector<TreeNode> trees;
+};
+
+/// Cache-simulation totals across all analyzed results.
+struct SimResult {
+  std::size_t n_cached = 0, n_persistent = 0;
+  sequant::AsyCost cached_footprint, persistent_footprint;
+};
+
+#endif  // SEQUANT_UTILITIES_COST_ANALYSIS_COST_ANALYSIS_HPP
diff --git a/utilities/cost_analysis/examples/.gitignore b/utilities/cost_analysis/examples/.gitignore
new file mode 100644
index 0000000000..f40c9f9ab3
--- /dev/null
+++ b/utilities/cost_analysis/examples/.gitignore
@@ -0,0 +1,5 @@
+# Generated by running cost_analysis / the reference CTests against the examples
+# in this directory; the frozen reference reports are the *.md.expected files.
+*.md
+!*.md.expected
+*.tree.txt
diff --git a/utilities/cost_analysis/examples/R2.tree.txt.expected b/utilities/cost_analysis/examples/R2.tree.txt.expected
new file mode 100644
index 0000000000..672d15649d
--- /dev/null
+++ b/utilities/cost_analysis/examples/R2.tree.txt.expected
@@ -0,0 +1,3 @@
+(g(a_1,a_2,i_1,i_2) * 1)
+((g(a_3,a_4,i_1,i_2) * t(i_3,i_4,a_3,a_4)) * R(a_1,a_2,i_3,i_4))
+((g(a_3,a_4,i_1,i_2) * t(i_3,i_4,a_3,a_4)) * t(a_1,a_2,i_3,i_4))
diff --git a/utilities/cost_analysis/examples/ccsd_r2.inp b/utilities/cost_analysis/examples/ccsd_r2.inp
new file mode 100644
index 0000000000..408e4c7528
--- /dev/null
+++ b/utilities/cost_analysis/examples/ccsd_r2.inp
@@ -0,0 +1 @@
+R{a1,a2;i1,i2} = g{a1,a2;i1,i2} + g{a3,a4;i1,i2} t{i3,i4;a3,a4} R{a1,a2;i3,i4} + g{a3,a4;i1,i2} t{i3,i4;a3,a4} t{a1,a2;i3,i4}
diff --git a/utilities/cost_analysis/examples/ccsd_r2.json b/utilities/cost_analysis/examples/ccsd_r2.json
new file mode 100644
index 0000000000..e6a77f52b8
--- /dev/null
+++ b/utilities/cost_analysis/examples/ccsd_r2.json
@@ -0,0 +1,30 @@
+{
+  "context": {
+    "spbasis": "spinor",
+    "field": "real"
+  },
+  "sizes": {
+    "i": 10,
+    "a": 38
+  },
+  "optimize": {
+    "objective": "dense_flops",
+    "reorder": true,
+    "volatile_leaf": "R"
+  },
+  "cache": {
+    "enabled": true,
+    "min_repeats": 2
+  },
+  "output": {
+    "path": "ccsd_r2.md",
+    "top_n": 20,
+    "dump_tree": true
+  },
+  "equations": [
+    {
+      "name": "R2",
+      "equation_file": "ccsd_r2.inp"
+    }
+  ]
+}
diff --git a/utilities/cost_analysis/examples/ccsd_r2.md.expected b/utilities/cost_analysis/examples/ccsd_r2.md.expected
new file mode 100644
index 0000000000..b76065f994
--- /dev/null
+++ b/utilities/cost_analysis/examples/ccsd_r2.md.expected
@@ -0,0 +1,36 @@
+# SeQuant cost analysis
+
+Sizes in MB (8 bytes/element).
+
+## R2
+
+Terms: 3; distinct intermediates: 1; reused: 1; largest: 0.0763 MB; peak storage: 2.2797 MB.
+
+Total operations (symbolic): 8*i^4a^2 + i^2a^2
+
+### Largest intermediates
+
+| Rank | Result | Spaces | Memory (MB) | Uses | Construction | Local ops |
+|---|---|---|---|---|---|---|
+| 1 | I(i_4,i_3,i_2,i_1) | OOOO | 0.0763 | 2 | (g(a_3,a_4,i_1,i_2) * t(i_3,i_4,a_3,a_4)) | 2*i^4a^2 |
+
+### Most expensive contractions
+
+| Rank | Result | Spaces | Uses | Construction | Operations |
+|---|---|---|---|---|---|
+| 1 | I(i_4,i_3,i_2,i_1) | OOOO | 2 | (g(a_3,a_4,i_1,i_2) * t(i_3,i_4,a_3,a_4)) | 2*i^4a^2 |
+
+### Shape census
+
+| Shape | Size (MB) | # Distinct | Total uses |
+|---|---|---|---|
+| OOOO | 0.0763 | 1 | 2 |
+
+## Cache (gated, volatile leaf = "R", min_repeats = 2)
+
+| Metric | Value |
+|---|---|
+| Cached | 1 |
+| Persistent | 1 |
+| Persistent footprint (MB) | 0.0763 |
+| Total cached footprint (MB) | 0.0763 |
diff --git a/utilities/cost_analysis/examples/df_r1.inp b/utilities/cost_analysis/examples/df_r1.inp
new file mode 100644
index 0000000000..4f480e24ef
--- /dev/null
+++ b/utilities/cost_analysis/examples/df_r1.inp
@@ -0,0 +1 @@
+R{a1,a2;i1,i2} = B{a1,i1;Κ1} B{a2,i2;Κ2} c{Κ1,Κ2}
diff --git a/utilities/cost_analysis/examples/df_r1.json b/utilities/cost_analysis/examples/df_r1.json
new file mode 100644
index 0000000000..6ad0d96440
--- /dev/null
+++ b/utilities/cost_analysis/examples/df_r1.json
@@ -0,0 +1,8 @@
+{
+  "context": { "spbasis": "spinor", "field": "complex", "convention": "sr", "aux": ["df"] },
+  "sizes": { "i": 10, "a": 38, "Κ": 400 },
+  "optimize": { "objective": "dense_flops", "reorder": true, "volatile_leaf": "R" },
+  "cache": { "enabled": true, "min_repeats": 2 },
+  "output": { "path": "df_r1.md", "top_n": 20 },
+  "equations": [ { "name": "R1", "equation_file": "df_r1.inp" } ]
+}
diff --git a/utilities/cost_analysis/examples/df_r1.md.expected b/utilities/cost_analysis/examples/df_r1.md.expected
new file mode 100644
index 0000000000..0741c67c49
--- /dev/null
+++ b/utilities/cost_analysis/examples/df_r1.md.expected
@@ -0,0 +1,36 @@
+# SeQuant cost analysis
+
+Sizes in MB (16 bytes/element).
+
+## R1
+
+Terms: 1; distinct intermediates: 1; reused: 0; largest: 2.3193 MB; peak storage: 6.8420 MB.
+
+Total operations (symbolic): 2*i^2a^2Κ + 2*iaΚ^2
+
+### Largest intermediates
+
+| Rank | Result | Spaces | Memory (MB) | Uses | Construction | Local ops |
+|---|---|---|---|---|---|---|
+| 1 | I(i_1,a_1,Κ_2) | OVX | 2.3193 | 1 | (B(a_1,i_1,Κ_1) * c(Κ_1,Κ_2)) | 2*iaΚ^2 |
+
+### Most expensive contractions
+
+| Rank | Result | Spaces | Uses | Construction | Operations |
+|---|---|---|---|---|---|
+| 1 | I(i_1,a_1,Κ_2) | OVX | 1 | (B(a_1,i_1,Κ_1) * c(Κ_1,Κ_2)) | 2*iaΚ^2 |
+
+### Shape census
+
+| Shape | Size (MB) | # Distinct | Total uses |
+|---|---|---|---|
+| OVX | 2.3193 | 1 | 1 |
+
+## Cache (gated, volatile leaf = "R", min_repeats = 2)
+
+| Metric | Value |
+|---|---|
+| Cached | 0 |
+| Persistent | 0 |
+| Persistent footprint (MB) | 0.0000 |
+| Total cached footprint (MB) | 0.0000 |
diff --git a/utilities/cost_analysis/report.cpp b/utilities/cost_analysis/report.cpp
new file mode 100644
index 0000000000..736476f87f
--- /dev/null
+++ b/utilities/cost_analysis/report.cpp
@@ -0,0 +1,145 @@
+#include "report.hpp"
+
+#include <SeQuant/core/eval/eval_node.hpp>
+#include <SeQuant/version.hpp>
+
+#include <algorithm>
+#include <format>
+#include <map>
+
+using namespace sequant;
+
+namespace {
+// Bytes per stored scalar (real=8, complex=16), threaded from the field.
+double mb(const AsyCost& c, double bytes_per_elem) {
+  return c.ops() * bytes_per_elem / (1024.0 * 1024.0);
+}
+
+void report_largest(const CellResult& cell, std::size_t top_n,
+                    double bytes_per_elem, std::ostream& out) {
+  std::vector<std::pair<const TreeNode*, const Record*>> rows;
+  for (const auto& [node, rec] : cell.catalog)
+    if (rec.memory != AsyCost::zero()) rows.emplace_back(&node, &rec);
+  // Sort by size (largest first); tie-break on the construction string so the
+  // row order is deterministic regardless of unordered_map iteration order.
+  std::ranges::sort(rows, [bytes_per_elem](const auto& a, const auto& b) {
+    const double sa = mb(a.second->memory, bytes_per_elem),
+                 sb = mb(b.second->memory, bytes_per_elem);
+    if (sa != sb) return sb < sa;
+    return full_expr(*a.first) < full_expr(*b.first);
+  });
+  out << "| Rank | Result | Spaces | Memory (MB) | Uses | Construction | "
+         "Local ops |\n|---|---|---|---|---|---|---|\n";
+  for (std::size_t i = 0; i < rows.size() && i < top_n; ++i) {
+    const auto& [node, rec] = rows[i];
+    out << std::format("| {} | {} | {} | {:.4f} | {} | {} | {} |\n", i + 1,
+                       rec->label, rec->spaces, mb(rec->memory, bytes_per_elem),
+                       rec->uses, full_expr(*node), rec->flops.text());
+  }
+}
+
+void report_expensive(const CellResult& cell, std::size_t top_n,
+                      std::ostream& out) {
+  std::vector<std::pair<const TreeNode*, const Record*>> rows;
+  for (const auto& [node, rec] : cell.catalog)
+    if (rec.flops != AsyCost::zero()) rows.emplace_back(&node, &rec);
+  // Most operations first, tie-broken on the construction string for
+  // determinism. The op count only orders rows
+  std::ranges::sort(rows, [](const auto& a, const auto& b) {
+    const double fa = a.second->flops.ops(), fb = b.second->flops.ops();
+    if (fa != fb) return fb < fa;
+    return full_expr(*a.first) < full_expr(*b.first);
+  });
+  out << "| Rank | Result | Spaces | Uses | Construction | "
+         "Operations |\n|---|---|---|---|---|---|\n";
+  for (std::size_t i = 0; i < rows.size() && i < top_n; ++i) {
+    const auto& [node, rec] = rows[i];
+    out << std::format("| {} | {} | {} | {} | {} | {} |\n", i + 1, rec->label,
+                       rec->spaces, rec->uses, full_expr(*node),
+                       rec->flops.text());
+  }
+}
+
+void report_shape_histogram(const CellResult& cell, double bytes_per_elem,
+                            std::ostream& out) {
+  struct Agg {
+    std::size_t count = 0, uses = 0;
+    AsyCost mem;
+  };
+  std::map<std::string, Agg> by_shape;
+  for (const auto& [_, rec] : cell.catalog) {
+    if (rec.memory == AsyCost::zero()) continue;
+    auto& g = by_shape[rec.spaces];
+    ++g.count;
+    g.uses += rec.uses;
+    // Representative size of the shape = its largest instance. A single O/V/X
+    // signature can group differing sizes only when multiple aux spaces (Κ, L)
+    // are registered; taking the max keeps this deterministic either way.
+    if (mb(g.mem, bytes_per_elem) < mb(rec.memory, bytes_per_elem))
+      g.mem = rec.memory;
+  }
+  std::vector<std::pair<std::string, Agg>> rows(by_shape.begin(),
+                                                by_shape.end());
+  std::ranges::sort(rows, [bytes_per_elem](const auto& a, const auto& b) {
+    const double sa = mb(a.second.mem, bytes_per_elem),
+                 sb = mb(b.second.mem, bytes_per_elem);
+    if (sa != sb) return sb < sa;
+    return a.first < b.first;  // tie-break on shape string
+  });
+  out << "| Shape | Size (MB) | # Distinct | Total uses |\n|---|---|---|---|\n";
+  for (const auto& [shape, g] : rows)
+    out << std::format("| {} | {:.4f} | {} | {} |\n", shape,
+                       mb(g.mem, bytes_per_elem), g.count, g.uses);
+}
+
+}  // namespace
+
+std::string full_expr(const TreeNode& n) {
+  if (n.leaf()) return n->label();
+  if (n->op_type() == EvalOp::Adjoint) return full_expr(n.left()) + "+";
+  const char* op = n->op_type() == EvalOp::Product ? " * " : " + ";
+  return "(" + full_expr(n.left()) + op + full_expr(n.right()) + ")";
+}
+
+void write_report(
+    const Config& cfg,
+    const std::vector<std::pair<std::string, CellResult>>& results,
+    const SimResult& sim, std::ostream& out) {
+  const double bytes_per_elem = cfg.real_field ? 8.0 : 16.0;
+  // Sections are separated by a leading blank line rather than a trailing one,
+  // so the report ends with a single newline (matches the committed fixtures).
+  out << "# SeQuant cost analysis\n";
+  if (!cfg.out.omit_revision)
+    out << std::format("\n_SeQuant: {}_\n", git_revision());
+  out << std::format("\nSizes in MB ({} bytes/element).\n", bytes_per_elem);
+  for (const auto& [name, cell] : results) {
+    out << std::format("\n## {}\n", name);
+    out << std::format(
+        "\nTerms: {}; distinct intermediates: {}; reused: {}; largest: {:.4f} "
+        "MB; peak storage: {:.4f} MB.\n",
+        cell.n_terms, cell.n_distinct, cell.n_reused,
+        mb(cell.largest_mem, bytes_per_elem),
+        mb(cell.peak_storage, bytes_per_elem));
+    // Op-count only; no concrete FLOP number.
+    out << std::format("\nTotal operations (symbolic): {}\n",
+                       cell.total_flops.text());
+    out << "\n### Largest intermediates\n\n";
+    report_largest(cell, cfg.out.top_n, bytes_per_elem, out);
+    out << "\n### Most expensive contractions\n\n";
+    report_expensive(cell, cfg.out.top_n, out);
+    out << "\n### Shape census\n\n";
+    report_shape_histogram(cell, bytes_per_elem, out);
+  }
+  if (cfg.cache.enabled) {
+    out << std::format(
+        "\n## Cache (gated, volatile leaf = \"{}\", min_repeats = {})\n\n",
+        cfg.optimize.volatile_leaf, cfg.cache.min_repeats);
+    out << "| Metric | Value |\n|---|---|\n";
+    out << std::format("| Cached | {} |\n", sim.n_cached);
+    out << std::format("| Persistent | {} |\n", sim.n_persistent);
+    out << std::format("| Persistent footprint (MB) | {:.4f} |\n",
+                       mb(sim.persistent_footprint, bytes_per_elem));
+    out << std::format("| Total cached footprint (MB) | {:.4f} |\n",
+                       mb(sim.cached_footprint, bytes_per_elem));
+  }
+}
diff --git a/utilities/cost_analysis/report.hpp b/utilities/cost_analysis/report.hpp
new file mode 100644
index 0000000000..251a79dffd
--- /dev/null
+++ b/utilities/cost_analysis/report.hpp
@@ -0,0 +1,22 @@
+#ifndef SEQUANT_UTILITIES_COST_ANALYSIS_REPORT_HPP
+#define SEQUANT_UTILITIES_COST_ANALYSIS_REPORT_HPP
+
+#include "cost_analysis.hpp"
+
+#include <ostream>
+#include <string>
+#include <utility>
+#include <vector>
+
+/// Serializes a binarized tree back to a parenthesized string (also used by the
+/// pipeline's --dump_tree output).
+std::string full_expr(const TreeNode& n);
+
+/// Writes the full Markdown report (per-result tables + cache section) to \p
+/// out.
+void write_report(
+    const Config& cfg,
+    const std::vector<std::pair<std::string, CellResult>>& results,
+    const SimResult& sim, std::ostream& out);
+
+#endif  // SEQUANT_UTILITIES_COST_ANALYSIS_REPORT_HPP