diff --git a/docs/plasticity.md b/docs/plasticity.md
new file mode 100644
index 0000000..482ca08
--- /dev/null
+++ b/docs/plasticity.md
@@ -0,0 +1,359 @@
+# Plasticity materials
+
+How the plasticity models are put together after the split of
+`small_strain_plasticity`, what each one requires, and which assumptions are
+load-bearing.
+
+Current as of `refactor/scalar-newton-split` (PRs #39 and #40).
+
+---
+
+## 1. The layout
+
+```
+                    ┌──────────────────────┐
+                    │   newton_scalar<T>   │   plain algorithm
+                    │  no graph, no props  │   returns {x, converged, iterations}
+                    └──────────┬───────────┘
+                     used by   │
+              ┌────────────────┴────────────────┐
+              │                                 │
+   ┌──────────▼──────────┐          ┌───────────▼───────────┐
+   │   backward_euler    │          │     local_newton      │
+   │  MATERIAL, graph-   │          │  MATERIAL, driven by  │
+   │  driven             │          │  another material     │
+   │  "function" REQUIRED│          │  exposes solve(eval)  │
+   └──────────┬──────────┘          └───────────┬───────────┘
+              │                                 │  material_ref
+   ┌──────────▼──────────┐          ┌───────────▼───────────┐
+   │ autocatalytic_      │          │  j2_plasticity        │
+   │ reaction, curing    │          │  drucker_prager_      │
+   │                     │          │  plasticity           │
+   └─────────────────────┘          └───────────────────────┘
+
+   j2_rk_plasticity iterates its own Butcher tableau and uses neither.
+```
+
+Three plasticity materials, none of them templated on a yield function any
+more:
+
+| material | integrator | yield surface | flow |
+|---|---|---|---|
+| `j2_plasticity` | backward Euler (radial return) | von Mises cylinder | associative |
+| `drucker_prager_plasticity` | backward Euler + apex branch | DP cone | **non**-associative (β ≠ η) |
+| `j2_rk_plasticity` | Runge–Kutta, any Butcher tableau | von Mises cylinder | associative |
+
+---
+
+## 2. Interfaces
+
+### `j2_plasticity`
+
+| parameter | meaning |
+|---|---|
+| `hardening_source` | material publishing `hardening_stress` and `hardening_modulus` |
+| `strain_source` | material publishing `strain` |
+| `solver_source` | a `local_newton` |
+| `K`, `G` | bulk and shear moduli — the elastic stiffness is built from these |
+| `sigma_0` | initial yield stress |
+
+Outputs: `stress`, `tangent`, and the history pair `plastic_strain`,
+`equivalent_plastic_strain`.
+
+### `drucker_prager_plasticity`
+
+Same, minus `K`, plus the cone:
+
+| parameter | meaning |
+|---|---|
+| `eta` | friction (pressure coefficient) |
+| `beta` | dilatancy — `beta != eta` is what makes the flow non-associative |
+| `K_bulk` | bulk modulus, used both for the cone apex and for the elastic stiffness |
+
+These used to arrive inside a `yield_function` **C++ object**, which the JSON
+reader cannot convert — so the material could not be configured from a document
+at all. As plain scalars it can (issue #33).
+
+### `j2_rk_plasticity`
+
+Takes `K`, `G`, `sigma_0`, the two sources, plus `tolerance`, `max_iter` and a
+`tableau` pointer. Available tableaus: `forward_euler`, `explicit_midpoint`,
+`rk4`, `implicit_euler`, `implicit_midpoint`, `crank_nicolson`, `sdirk3`,
+`gauss_legendre_4`.
+
+### Hardening
+
+Both publish `hardening_stress` (H) and `hardening_modulus` (dH/dκ), and read
+`equivalent_plastic_strain` from the plasticity material — a `Local` edge,
+because H depends on κ, which is the unknown being solved for.
+
+| material | law | parameters |
+|---|---|---|
+| `linear_isotropic_hardening` | `H = K κ` | `source`, `K` |
+| `exponential_isotropic_hardening` | `H = K_inf (1 − e^{−δκ})` | `source`, `K_inf`, `delta` |
+
+---
+
+## 3. What is assumed, and where it bites
+
+### Isotropic elasticity is required, not preferred
+
+Every closed form below rests on `C_e` being isotropic:
+
+```
+C_e : N        = 2G dev(N) + K tr(N) I
+N : C_e : N    = 3G                        (J2, N deviatoric)
+C_e : X : C_e  = 4G² X                     (X deviatoric in both index pairs)
+```
+
+This is why the plasticity materials **build their own** `C_e` from `K` and `G`
+rather than reading a rank-4 tangent from an elastic material. Accepting an
+arbitrary `C_e` advertised a generality none of them can honour — a caller
+wiring an anisotropic tangent would have got silently wrong answers.
+
+It also keeps `linear_elasticity` out of plasticity graphs, which matters for a
+second reason: its `stress` output is `C : ε` with `ε_p` ignored, so inside a
+plasticity graph it is not the stress of anything and diverges as plastic strain
+accumulates.
+
+| step | κ | `stress` (real) | `elastic::stress` | error |
+|---|---|---|---|---|
+| 40 | 0.0094 | 106.25 | 107.69 | +1.4% |
+| 120 | 0.1094 | 306.25 | 323.08 | +5.5% |
+| 200 | 0.2094 | 506.25 | 538.46 | **+6.4%** |
+
+`linear_elasticity` remains a valid material — for a genuinely elastic model
+`C : ε` *is* the answer, and `isotropic_damage` consumes both its stress and its
+tangent legitimately. It is simply the wrong dependency for plasticity.
+
+### Clamps belong to the caller
+
+`newton_scalar` does not clamp its result. `max(x, 0)` is a statement about a
+plastic multiplier and `abs(x)` about a curing degree; neither is about Newton's
+method. Each lives in the material that owns the assumption (issue #13).
+
+### The return map solves conditionally; a graph property is evaluated unconditionally
+
+This is why `local_newton` exists rather than everything using
+`backward_euler`.
+
+A property's update callback runs whenever the graph updates. **Plasticity only
+solves when the trial state exceeds yield** — and in a real analysis most
+integration points are elastic most of the time. Both return maps short-circuit
+before touching the solver:
+
+```cpp
+if (sig_eq - sigma_0 - H <= 0) { /* elastic: stress = trial, tangent = C_e */ return; }
+```
+
+Making the solve a property callback would run a Newton at every elastic point:
+
+| | measured |
+|---|---|
+| elastic step | ~92 ns |
+| plastic step | ~168 ns |
+
+Roughly 80 % more work, paid exactly where a real analysis spends most of its
+time.
+
+It is worse than a cost. At `Δλ = 0` on an elastic step the residual is
+negative, so Newton drives `Δλ` negative; holding it at zero needs a clamp —
+which is precisely the `max(x, 0)` that issue #13 objects to. **The clamp and
+the graph-driving are the same problem**: the graph mode has no way to express
+*do not solve*.
+
+Drucker-Prager cannot be graph-driven for a second, independent reason: it
+solves up to **twice** per update — an apex pre-check, a smooth-cone solve, then
+an apex fallback if that fails to converge — and chooses the branch on the first
+solve's convergence. A property edge carries one number, not "and it failed, so
+take the other branch".
+
+J2 solves once and has no such branch, so *only* the conditional argument
+applies to it. It could be graph-driven; it would simply cost more than it
+saves, and would leave two patterns where there is now one.
+
+Convergence therefore travels *with* the result (`{x, converged, iterations}`)
+rather than being queried from the solver afterwards, where it went stale
+between Drucker-Prager's two solves.
+
+---
+
+## 4. The consistent tangent
+
+For the smooth return, with `M` the yield normal and `N` the flow normal
+(`M == N` when associative):
+
+```
+A            = C_e − Δλ (C_e : dN/dσ : C_e)
+dλ/dε        = (M : C_e) / (M : C_e : N + H′)
+C_consistent = A − (C_e : N) ⊗ dλ/dε
+```
+
+`M ≠ N` under non-associative flow makes the tangent **major-asymmetric**,
+measured at 0.6–0.8 % of peak magnitude on multiaxial DP paths.
+
+J2 collapses this to the standard closed form:
+
+```
+C = C_e − (6G²Δλ/σ_eq) IIdev + (4G²Δλ/σ_eq − 4G²/(3G+H′)) N ⊗ N
+```
+
+which is algebraically identical, not an approximation — confirmed by
+bit-identical stress against the general path.
+
+### Apex return (Drucker-Prager only)
+
+When the deviatoric correction would overshoot the cone tip
+(`G Δλ ≥ √J₂`), the return goes to the apex instead: deviatoric stress to zero,
+pressure pinned at `(k + H)/η`, volumetric plastic flow from `β`. Its tangent is
+a **branch** tangent — valid only for perturbations staying on the apex — because
+the return map is non-smooth there.
+
+---
+
+## 5. Verification
+
+`tangent_checker` compares the analytical tangent against a central difference
+through the graph. Coverage:
+
+| path | purpose |
+|---|---|
+| uniaxial | the historical case |
+| pure shear, biaxial, triaxial, mixed dev+shear | exercise genuine non-associativity |
+| hydrostatic | the **only** path that reaches the apex branch |
+
+All agree to ~1e-10, the apex to 8e-9.
+
+Two coverage lessons worth keeping:
+
+- **A single load path proves one path.** The apex return was executed by *no*
+  test for its entire existence — instrumenting the predicate gave
+  `APEX_HITS=0` across every binary — because every path was uniaxial and the
+  apex sits on the hydrostatic axis. It was unreachable by construction, not by
+  oversight.
+- **A tolerance set by the worst step licenses errors in all the others.**
+  `J2TangentTest` bounded a whole run at `0.1` because the elastic→plastic
+  transition step is genuinely inexact. A 0.5 % error injected into the tangent
+  landed at 2.4e-4 and passed, while real plastic steps sit at 4.4e-10. The two
+  regimes are now bounded separately.
+
+---
+
+## 6. Performance
+
+Two binaries built from the actual commits, run interleaved so each pair sees
+the same machine state; 15 pairs; speedup computed per pair.
+
+| | old (median) | new (median) | paired speedup | range |
+|---|---|---|---|---|
+| J2 | 1466.2 ns/step | **211.5** | **7.00×** | 5.54–8.84× |
+| Drucker-Prager | 1176.3 ns/step | **652.4** | **1.79×** | 1.58–1.90× |
+
+The range is the honest figure: the *same* binary measured 1171–2253 ns across
+runs on this machine, so any single-run comparison is worth about one
+significant digit.
+
+J2 gains more because it also sheds machinery it never used. DP keeps the
+non-associative structure and the apex branch and gains only the arithmetic —
+which is the correct outcome. **The generality DP pays for is generality DP
+uses.**
+
+---
+
+## 7. Building a model
+
+```cpp
+// solver — one instance can serve several materials
+p.insert<std::string>("name", "solver");
+ctx.create<local_newton<policy>>(p);
+
+// hardening — reads back from the plasticity material
+p.clear();
+p.insert<std::string>("name", "hardening");
+p.insert<std::string>("source", "j2");
+p.insert<T>("K", 1000.0);
+ctx.create<linear_isotropic_hardening<policy>>(p);
+
+// plasticity — no elastic material anywhere
+p.clear();
+p.insert<std::string>("name", "j2");
+p.insert<std::string>("hardening_source", "hardening");
+p.insert<std::string>("strain_source", "stepper");
+p.insert<std::string>("solver_source", "solver");
+p.insert<T>("K", 166.67);
+p.insert<T>("G", 76.92);
+p.insert<T>("sigma_0", 50.0);
+ctx.create<j2_plasticity<policy>>(p);
+```
+
+Note `hardening.source = "j2"` while `j2.hardening_source = "hardening"`. The
+cycle is deliberate and is why those edges are `Local`: H depends on κ, which is
+what the return map solves for, so the hardening material is re-evaluated inside
+the Newton loop through `update_source()`.
+
+The same model in JSON, which is possible for Drucker-Prager only since its cone
+parameters became plain scalars:
+
+```json
+{"type": "drucker_prager_plasticity", "name": "dp",
+ "hardening_source": "hardening", "strain_source": "stepper",
+ "solver_source": "solver",
+ "G": 76.92, "sigma_0": 20.0,
+ "eta": 0.1, "beta": 0.05, "K_bulk": 166.67}
+```
+
+---
+
+## 8. What is still generic, and why
+
+`plasticity_utils`' free functions — `compute_trial`, `evaluate_at_state`,
+`compute_tangent` — remain templated on a yield function, and that generality is
+real: `drucker_prager_plasticity` instantiates them with the DP cone and
+`j2_rk_plasticity` with the von Mises cylinder. Two callers, two yield
+functions, shared return-mapping algebra.
+
+That is the distinction worth holding onto. A template parameter with **one**
+argument is indirection — it was removed from `small_strain_plasticity` and from
+`rk_plasticity`. A template parameter with two genuinely different arguments is
+what templates are for.
+
+## 9. Settled: the solver is reached by `material_ref`, not by the graph
+
+The return maps hold a `material_ref<local_newton>` and call `solve()`. The
+alternative — publishing a `yielding` flag so a graph-driven `backward_euler`
+knows when to iterate — was considered and **rejected**.
+
+The pattern would work: `strain_threshold_yield` already publishes
+`is_yielding` and `isotropic_damage` consumes it, so a gating flag is native to
+this codebase. It would retire `local_newton` and, with it, `material_ref` —
+about 99 lines of core machinery whose only two call sites are these.
+
+It was rejected on what it would cost to express:
+
+- **The material splits into phases.** `compute()` currently does trial →
+  solve → update in one callback. Graph-driven it becomes publish
+  `yielding`/`residual`/`jacobian`, solver runs, read `delta`, update. Two or
+  three properties where there is one, and graph dispatch measures ~30 ns even
+  for a trivial graph.
+- **A flag cannot carry Drucker-Prager's branch.** The apex fallback fires when
+  the *smooth solve fails*, which is not known until after it runs. Expressible
+  only as trial → smooth solver → `smooth_converged` → apex solver gated on
+  that flag → update: two solver instances, three flags, four phases, replacing
+  `if (!smooth.converged) do_apex_return(...)`.
+- **The elastic saving is partial** anyway — the solver's callback still fires
+  and returns early on the flag.
+
+The cost of the machinery is real but bounded; the cost of removing it is a
+graph harder to read than the code it replaces.
+
+**One constraint this decision carries.** `material_ref` bypasses the
+topological sort, so the plasticity↔solver ordering is not an edge the engine
+knows about. That is safe only because `local_newton` holds **no per-solve
+state** — its `solve()` is `const` and returns everything it computes. Give it
+mutable state and the ordering becomes real and unenforced.
+
+## 10. Known gaps
+
+- **The apex tangent is a branch tangent.** Verified consistent *on* the branch;
+  a perturbation that leaves the apex back onto the smooth cone is not covered,
+  and cannot be by a central difference.
diff --git a/include/numsim-materials/materials/drucker_prager_plasticity.h b/include/numsim-materials/materials/drucker_prager_plasticity.h
index 7cf5885..ebfffad 100644
--- a/include/numsim-materials/materials/drucker_prager_plasticity.h
+++ b/include/numsim-materials/materials/drucker_prager_plasticity.h
@@ -3,6 +3,7 @@
 
 #include <cmath>
 #include <concepts>
+#include <algorithm>
 #include <stdexcept>
 #include <utility>
 #include <tmech/tmech.h>
@@ -10,7 +11,7 @@
 #include "numsim-materials/core/material_ref.h"
 #include "numsim-materials/materials/drucker_prager_yield_function.h"
 #include "numsim-materials/materials/plasticity_utils.h"
-#include "numsim-materials/solvers/backward_euler.h"
+#include "numsim-materials/solvers/local_newton.h"
 
 namespace numsim::materials {
 
@@ -29,7 +30,7 @@ class drucker_prager_plasticity final
   using tensor2 = tmech::tensor<value_type, Dim, 2>;
   using tensor4 = tmech::tensor<value_type, Dim, 4>;
   using yield_fn = drucker_prager_yield_function<value_type, base::Dim>;
-  using solver_type = backward_euler<Traits>;
+  using solver_type = local_newton<Traits>;
 
   template <typename... Args>
   explicit drucker_prager_plasticity(Args&&... args)
@@ -106,18 +107,20 @@ class drucker_prager_plasticity final
       return;
     }
 
-    const auto dlambda = solve_smooth_newton(ts.eval.modified_sig_eq, kappa_n);
+    const auto smooth = solve_smooth_newton(ts.eval.modified_sig_eq, kappa_n);
 
-    // If smooth Newton fails and apex is available, try apex as fallback.
-    if (!m_solver.get().converged()) {
-      do_apex_return(ts.eval.sig, C_e, kappa_n);
-      if (!m_solver.get().converged())
+    // If the smooth return fails, the apex is the remaining branch.
+    if (!smooth.converged) {
+      if (!do_apex_return(ts.eval.sig, C_e, kappa_n))
         throw std::runtime_error(
             "drucker_prager_plasticity: both smooth and apex Newton failed");
       return;
     }
 
-    do_smooth_return(ts.eval, C_e, kappa_n, dlambda);
+    // dlambda >= 0 is a statement about the plastic multiplier, enforced here
+    // rather than inside a general scalar solver (see #13).
+    do_smooth_return(ts.eval, C_e, kappa_n,
+                     std::max(smooth.x, value_type{0}));
   }
 
 private:
@@ -131,8 +134,12 @@ class drucker_prager_plasticity final
 
   /// Scalar Newton solve: r(Δλ) = phi - G_eff·Δλ - Y0 - H(κ_n + Δλ) = 0.
   /// Used for both the smooth and apex returns with different (phi, G_eff).
-  value_type solve_scalar_return(value_type phi_trial, value_type G_eff,
-                                  value_type kappa_n) {
+  /// Returns the solver's result, not a bare number: convergence travels WITH
+  /// the value instead of being queried from the solver afterwards, where it
+  /// went stale between the smooth and apex solves.
+  typename solver_type::result solve_scalar_return(value_type phi_trial,
+                                                   value_type G_eff,
+                                                   value_type kappa_n) {
     auto eval = [&](value_type dl) -> std::pair<value_type, value_type> {
       m_kappa.new_value() = kappa_n + dl;
       m_H.update_source();
@@ -143,7 +150,8 @@ class drucker_prager_plasticity final
   }
 
   /// Smooth-cone return Newton: phi = modified_sig_eq, G_eff from yield function.
-  value_type solve_smooth_newton(value_type phi_trial, value_type kappa_n) {
+  typename solver_type::result solve_smooth_newton(value_type phi_trial,
+                                                   value_type kappa_n) {
     return solve_scalar_return(phi_trial, m_yf.effective_modulus(m_G), kappa_n);
   }
 
@@ -172,11 +180,14 @@ class drucker_prager_plasticity final
 
   /// Apex return: deviatoric stress vanishes, only volumetric Newton.
   /// dev(ε_p) = dev(ε), tr(ε_p) += β·Δκ. Tangent is rank-1 volumetric.
-  void do_apex_return(const tensor2& sig_trial, const tensor4& C_e,
+  /// @return whether the apex Newton converged.
+  bool do_apex_return(const tensor2& sig_trial, const tensor4& C_e,
                       value_type kappa_n) {
     const auto phi_apex = m_yf.apex_modified_sig_eq(sig_trial);
     const auto G_eff_apex = m_yf.apex_effective_modulus();
-    const auto dkappa = solve_scalar_return(phi_apex, G_eff_apex, kappa_n);
+    const auto sol = solve_scalar_return(phi_apex, G_eff_apex, kappa_n);
+    if (!sol.converged) return false;
+    const auto dkappa = std::max(sol.x, value_type{0});
 
     m_eps_p.new_value() = m_yf.apex_plastic_strain(
         m_strain.get(), m_eps_p.old_value(), dkappa);
@@ -187,6 +198,7 @@ class drucker_prager_plasticity final
     // remain on the active apex branch (q stays at 0).
     m_H.update_source();
     m_tangent = m_yf.apex_tangent(m_dH.get());
+    return true;
   }
 
 private:
diff --git a/include/numsim-materials/materials/j2_plasticity.h b/include/numsim-materials/materials/j2_plasticity.h
index ac3959d..265604f 100644
--- a/include/numsim-materials/materials/j2_plasticity.h
+++ b/include/numsim-materials/materials/j2_plasticity.h
@@ -2,13 +2,14 @@
 #define J2_PLASTICITY_H
 
 #include <cmath>
+#include <algorithm>
 #include <stdexcept>
 #include <string>
 #include <tmech/tmech.h>
 #include "numsim-materials/core/material_base.h"
 #include "numsim-materials/core/material_ref.h"
 #include "numsim-materials/materials/plasticity_utils.h"
-#include "numsim-materials/solvers/backward_euler.h"
+#include "numsim-materials/solvers/local_newton.h"
 
 namespace numsim::materials {
 
@@ -59,7 +60,7 @@ class j2_plasticity final
   using base::Dim;
   using tensor2 = tmech::tensor<value_type, Dim, 2>;
   using tensor4 = tmech::tensor<value_type, Dim, 4>;
-  using solver_type = backward_euler<Traits>;
+  using solver_type = local_newton<Traits>;
 
   template <typename... Args>
   explicit j2_plasticity(Args&&... args)
@@ -150,10 +151,13 @@ class j2_plasticity final
       return {sig_eq - G_eff * dl - m_sigma_0 - m_H.get(),
               -G_eff - m_dH.get()};
     };
-    const auto dlambda = m_solver.get().solve(eval);
-    if (!m_solver.get().converged())
+    const auto sol = m_solver.get().solve(eval);
+    if (!sol.converged)
       throw std::runtime_error(
           "j2_plasticity: return-mapping Newton failed to converge");
+    // dlambda >= 0 is a statement about the plastic multiplier, so it is
+    // enforced here rather than inside a general scalar solver (see #13).
+    const auto dlambda = std::max(sol.x, value_type{0});
 
     m_eps_p.new_value() = m_eps_p.old_value() + dlambda * N;
     m_kappa.new_value() = kappa_n + dlambda;
diff --git a/include/numsim-materials/materials/rk_plasticity.h b/include/numsim-materials/materials/j2_rk_plasticity.h
similarity index 93%
rename from include/numsim-materials/materials/rk_plasticity.h
rename to include/numsim-materials/materials/j2_rk_plasticity.h
index a072e5a..e796536 100644
--- a/include/numsim-materials/materials/rk_plasticity.h
+++ b/include/numsim-materials/materials/j2_rk_plasticity.h
@@ -1,5 +1,5 @@
-#ifndef RK_PLASTICITY_H
-#define RK_PLASTICITY_H
+#ifndef J2_RK_PLASTICITY_H
+#define J2_RK_PLASTICITY_H
 
 #include <cmath>
 #include <vector>
@@ -17,23 +17,23 @@ namespace numsim::materials {
 /// Each implicit stage solves F = 0 for Δλ_i. Explicit stages use
 /// the consistency condition. Shares trial/tangent code with
 /// small_strain_plasticity via plasticity_utils.h.
-template<typename Traits, typename YieldFunction>
-class rk_plasticity final
-    : public material_base<rk_plasticity<Traits, YieldFunction>, Traits> {
+template<typename Traits>
+class j2_rk_plasticity final
+    : public material_base<j2_rk_plasticity<Traits>, Traits> {
 public:
-  using base = material_base<rk_plasticity<Traits, YieldFunction>, Traits>;
+  using base = material_base<j2_rk_plasticity<Traits>, Traits>;
   using value_type = typename base::value_type;
   using input_parameter_controller = typename base::input_parameter_controller;
   static constexpr auto Dim = base::Dim;
   using tensor2 = tmech::tensor<value_type, Dim, 2>;
   using tensor4 = tmech::tensor<value_type, Dim, 4>;
-  using yield_fn = YieldFunction;
+  using yield_fn = j2_yield_function<value_type, base::Dim>;
 
   template <typename... Args>
-  explicit rk_plasticity(Args&&... args)
+  explicit j2_rk_plasticity(Args&&... args)
       : base(std::forward<Args>(args)...),
         m_stress(base::template add_output<tensor2>(
-            "stress", &rk_plasticity::compute)),
+            "stress", &j2_rk_plasticity::compute)),
         m_tangent(base::template add_output<tensor4>("tangent")),
         m_eps_p(base::template add_history_output<tensor2>("plastic_strain")),
         m_kappa(base::template add_history_output<value_type>("equivalent_plastic_strain")),
@@ -211,10 +211,7 @@ class rk_plasticity final
   std::vector<double> m_diag;
 };
 
-template<typename Traits>
-using j2_rk_plasticity = rk_plasticity<Traits,
-    j2_yield_function<typename Traits::value_type, Traits::Dim>>;
 
 } // namespace numsim::materials
 
-#endif // RK_PLASTICITY_H
+#endif // J2_RK_PLASTICITY_H
diff --git a/include/numsim-materials/solvers/backward_euler.h b/include/numsim-materials/solvers/backward_euler.h
index e6707aa..0b378b3 100644
--- a/include/numsim-materials/solvers/backward_euler.h
+++ b/include/numsim-materials/solvers/backward_euler.h
@@ -2,20 +2,32 @@
 #define BACKWARD_EULER_H
 
 #include <cmath>
+#include <string>
 #include "numsim-materials/core/material_base.h"
+#include "numsim-materials/solvers/newton_scalar.h"
 
 namespace numsim::materials {
 
-/// Backward Euler solver as a material.
+/// Backward-Euler update solved by a scalar Newton the PROPERTY GRAPH drives.
 ///
-/// Consumes: function_name::residual, function_name::jacobian (Local edges)
-/// Produces: "delta" (converged increment)
+/// Reads "residual" and "jacobian" from the material named by "function", and
+/// publishes the increment as "delta". The function material re-evaluates its
+/// residual against the current delta through update_source(), which is the
+/// circularity that makes this work at all.
 ///
-/// The Newton iteration calls update_source() on the residual/jacobian inputs
-/// to re-evaluate them at each trial increment.
+/// "function" is REQUIRED. It used to default to an empty string, which
+/// silently selected a second, callback-driven mode inside this same class: no
+/// inputs were created, update() was never bound, and "delta" stayed at zero.
+/// A consumer reading it -- autocatalytic_reaction does -- then froze at its
+/// start value for the whole analysis with no error anywhere. Measured: a cure
+/// that should reach 1.0 sat at 0.01 for 30 steps. That mode is now
+/// local_newton, chosen by naming a different type rather than by omitting a
+/// parameter.
+///
+/// Parameters:
+///   "name", "function", "tolerance", "max_iter"
 template<typename Traits>
-class backward_euler final
-    : public material_base<backward_euler<Traits>, Traits> {
+class backward_euler final : public material_base<backward_euler<Traits>, Traits> {
 public:
   using base = material_base<backward_euler<Traits>, Traits>;
   using value_type = typename base::value_type;
@@ -26,89 +38,59 @@ class backward_euler final
       : base(std::forward<Args>(args)...),
         m_delta(base::template add_output<value_type>("delta")),
         m_func_name(base::template get_parameter<std::string>("function")),
-        m_tol(base::template get_parameter<value_type>("tolerance")),
-        m_max_iter(base::template get_parameter<int>("max_iter"))
+        m_solver(base::template get_parameter<value_type>("tolerance"),
+                 base::template get_parameter<int>("max_iter")),
+        m_residual(base::template add_input<value_type>(
+            m_func_name, "residual", EdgeKind::Local)),
+        m_jacobian(base::template add_input<value_type>(
+            m_func_name, "jacobian", EdgeKind::Local))
   {
-    // If a function name is provided, set up graph-driven iteration
-    if (!m_func_name.empty()) {
-      m_residual = &base::template add_input<value_type>(
-          m_func_name, "residual", EdgeKind::Local);
-      m_jacobian = &base::template add_input<value_type>(
-          m_func_name, "jacobian", EdgeKind::Local);
-      // Bind update callback for graph-driven mode
-      if (auto p = base::m_property_handler.find(base::m_name, "delta"))
-        (*p)->traits().update = [this]() { this->update(); };
-    }
+    if (auto p = base::m_property_handler.find(base::m_name, "delta"))
+      (*p)->traits().update = [this]() { this->update(); };
   }
 
   static input_parameter_controller parameters() {
     input_parameter_controller para{base::parameters()};
-    para.template insert<std::string>("function")
-        .template add<set_default>(std::string{});
+    para.template insert<std::string>("function").template add<is_required>();
     para.template insert<value_type>("tolerance")
-        .template add<set_default>(value_type{5e-12});
-    para.template insert<int>("max_iter")
-        .template add<set_default>(100);
+        .template add<set_default>(value_type{1e-10});
+    para.template insert<int>("max_iter").template add<set_default>(50);
     return para;
   }
 
-  /// Property-graph driven iteration: reads residual/jacobian via update_source.
-  /// Used when the solver is in the graph (e.g., curing reaction).
+  /// Whether the last update() converged.
+  ///
+  /// Previously never set on this path at all: the loop broke on tolerance, on
+  /// a singular jacobian and on exhausting its budget, and all three looked
+  /// identical from outside.
+  [[nodiscard]] bool converged() const noexcept { return m_converged; }
+
   void update() override {
-    if (!m_residual || !m_jacobian) return;
-    m_delta = value_type{5e-12};
-    for (int i = 0; i < m_max_iter; ++i) {
-      m_residual->update_source();
-      const auto& r = m_residual->get();
-      if (std::abs(r) <= m_tol) break;
-      m_jacobian->update_source();
-      const auto& j = m_jacobian->get();
-      if (std::abs(j) < value_type{1e-30}) break;
-      auto step = r / j;
-      // Damped Newton: halve step if it produces NaN
-      for (int k = 0; k < 5; ++k) {
-        auto candidate = m_delta - step;
-        m_delta = candidate;
-        m_residual->update_source();
-        auto r_new = m_residual->get();
-        if (!std::isnan(r_new) && !std::isinf(r_new)) break;
-        m_delta = candidate + step; // restore
-        step *= value_type{0.5};
-      }
-    }
-    // Ensure positive increment (curing degree can only increase)
-    m_delta = std::abs(m_delta);
-  }
+    // The seed is nonzero because a jacobian evaluated at exactly zero is
+    // singular for the rate laws this drives.
+    const auto r = m_solver.solve(
+        [this](value_type x) -> std::pair<value_type, value_type> {
+          m_delta = x;
+          m_residual.update_source();
+          const auto res = m_residual.get();
+          m_jacobian.update_source();
+          return {res, m_jacobian.get()};
+        },
+        value_type{5e-12});
 
-  /// Direct call: another material provides eval(x) → {residual, jacobian}.
-  /// Used when the caller drives the iteration (e.g., plasticity return mapping).
-  /// Sets m_converged to indicate whether the iteration converged.
-  /// The returned value is clamped to be non-negative — for plasticity, a
-  /// negative plastic-multiplier increment is unphysical (backward plastic flow).
-  template<typename Eval>
-  value_type solve(Eval&& eval, value_type x0 = value_type{0}) {
-    auto x = x0;
-    for (int i = 0; i < m_max_iter; ++i) {
-      auto [r, dr] = eval(x);
-      if (std::abs(r) < m_tol) { m_converged = true; return std::max(x, value_type{0}); }
-      if (std::abs(dr) < value_type{1e-30}) { m_converged = false; return std::max(x, value_type{0}); }
-      x -= r / dr;
-    }
-    m_converged = false;
-    return std::max(x, value_type{0});
+    m_converged = r.converged;
+    // Sign convention owned here rather than by the solver: the quantities this
+    // integrates -- degree of cure, and similar -- only increase.
+    m_delta = std::abs(r.x);
   }
 
-  /// Whether the last solve() call converged.
-  bool converged() const { return m_converged; }
-
 private:
   value_type& m_delta;
   const std::string& m_func_name;
-  const value_type& m_tol;
-  const int& m_max_iter;
-  const input_property<value_type, property_traits>* m_residual{nullptr};
-  const input_property<value_type, property_traits>* m_jacobian{nullptr};
-  bool m_converged{true};
+  newton_scalar<value_type> m_solver;
+  const input_property<value_type, property_traits>& m_residual;
+  const input_property<value_type, property_traits>& m_jacobian;
+  bool m_converged{false};
 };
 
 } // namespace numsim::materials
diff --git a/include/numsim-materials/solvers/local_newton.h b/include/numsim-materials/solvers/local_newton.h
new file mode 100644
index 0000000..0ae13d8
--- /dev/null
+++ b/include/numsim-materials/solvers/local_newton.h
@@ -0,0 +1,69 @@
+#ifndef LOCAL_NEWTON_H
+#define LOCAL_NEWTON_H
+
+#include <string>
+#include "numsim-materials/core/material_base.h"
+#include "numsim-materials/solvers/newton_scalar.h"
+
+namespace numsim::materials {
+
+/// A scalar Newton another MATERIAL drives, by calling solve() with its own
+/// residual.
+///
+/// The counterpart to backward_euler, which the property graph drives. The two
+/// used to be one class selected by whether a "function" parameter happened to
+/// be set, which meant an empty "function" silently chose this behaviour --
+/// leaving a graph-driven consumer reading a delta of 0 forever, with no error
+/// anywhere. Splitting them makes the choice a type rather than a defaulted
+/// string.
+///
+/// Publishes no properties and takes no inputs. It is a material only so its
+/// tolerance and iteration budget are configurable from the same document as
+/// everything else, and so one instance can serve several materials.
+///
+/// Used by the return maps, which cannot be graph-driven: they solve twice per
+/// update with different residuals (smooth cone, then apex) and choose the
+/// branch on the first solve's convergence. A graph edge carries one number,
+/// not that.
+///
+/// Parameters:
+///   "name", "tolerance", "max_iter"
+template <typename Traits>
+class local_newton final : public material_base<local_newton<Traits>, Traits> {
+public:
+  using base = material_base<local_newton<Traits>, Traits>;
+  using value_type = typename base::value_type;
+  using input_parameter_controller = typename base::input_parameter_controller;
+  using result = typename newton_scalar<value_type>::result;
+
+  template <typename... Args>
+  explicit local_newton(Args&&... args)
+      : base(std::forward<Args>(args)...),
+        m_solver(base::template get_parameter<value_type>("tolerance"),
+                 base::template get_parameter<int>("max_iter")) {}
+
+  static input_parameter_controller parameters() {
+    input_parameter_controller para{base::parameters()};
+    para.template insert<value_type>("tolerance")
+        .template add<set_default>(value_type{1e-10});
+    para.template insert<int>("max_iter").template add<set_default>(50);
+    return para;
+  }
+
+  /// Solve with the caller's residual. @p eval maps x -> {residual, jacobian}.
+  ///
+  /// The result carries its own convergence flag, so it cannot be read without
+  /// being available -- unlike a converged() queried separately, which a caller
+  /// can forget and which goes stale between solves.
+  template <typename Eval>
+  result solve(Eval&& eval, value_type x0 = value_type{}) const {
+    return m_solver.solve(std::forward<Eval>(eval), x0);
+  }
+
+private:
+  newton_scalar<value_type> m_solver;
+};
+
+}  // namespace numsim::materials
+
+#endif  // LOCAL_NEWTON_H
diff --git a/include/numsim-materials/solvers/newton_scalar.h b/include/numsim-materials/solvers/newton_scalar.h
new file mode 100644
index 0000000..5cd4aa4
--- /dev/null
+++ b/include/numsim-materials/solvers/newton_scalar.h
@@ -0,0 +1,73 @@
+#ifndef NEWTON_SCALAR_H
+#define NEWTON_SCALAR_H
+
+#include <cmath>
+#include <utility>
+
+namespace numsim::materials {
+
+/// Scalar Newton iteration. Plain algorithm: no properties, no graph presence,
+/// no material_base.
+///
+/// Split out of backward_euler, which had grown two Newton loops -- one driven
+/// by graph properties, one by a caller's lambda -- with separate damping and
+/// separate clamping. The algorithm is the same in both cases; only where the
+/// residual comes from differs, and that belongs to the caller.
+///
+/// Deliberately does NOT clamp its result. backward_euler forced x >= 0 for
+/// plasticity's dlambda and |x| for a curing degree; both are statements about
+/// a particular unknown, not about Newton's method, and a general solver that
+/// silently enforces one is wrong for every other caller (see #13). Callers
+/// clamp what they own.
+template <typename T>
+class newton_scalar {
+public:
+  struct result {
+    T x{};                  ///< the iterate, converged or not
+    bool converged{false};
+    int iterations{0};
+  };
+
+  newton_scalar(T tolerance, int max_iterations) noexcept
+      : m_tol(tolerance), m_max_iter(max_iterations) {}
+
+  /// @param eval  x -> {residual, jacobian}
+  /// @param x0    initial iterate
+  ///
+  /// Reports convergence rather than returning a bare number: a caller that
+  /// does not check gets a non-converged iterate, which is the same failure
+  /// backward_euler's graph path had -- there it was unreportable, here it is
+  /// merely unchecked.
+  template <typename Eval>
+  result solve(Eval&& eval, T x0 = T{}) const {
+    result r{x0, false, 0};
+    for (int i = 0; i < m_max_iter; ++i) {
+      const auto [residual, jacobian] = eval(r.x);
+      ++r.iterations;
+      if (std::abs(residual) < m_tol) {
+        r.converged = true;
+        return r;
+      }
+      // A vanishing jacobian is a stall, not a convergence.
+      if (std::abs(jacobian) < T{1e-30}) return r;
+      r.x -= residual / jacobian;
+    }
+    // Exhausted the budget: the last step may still have landed on the root,
+    // so the residual is checked once more rather than assumed bad.
+    const auto [residual, jacobian] = eval(r.x);
+    (void)jacobian;
+    r.converged = std::abs(residual) < m_tol;
+    return r;
+  }
+
+  [[nodiscard]] T tolerance() const noexcept { return m_tol; }
+  [[nodiscard]] int max_iterations() const noexcept { return m_max_iter; }
+
+private:
+  T m_tol;
+  int m_max_iter;
+};
+
+}  // namespace numsim::materials
+
+#endif  // NEWTON_SCALAR_H
diff --git a/tests/plot_data.cpp b/tests/plot_data.cpp
index b32db4b..9b9984e 100644
--- a/tests/plot_data.cpp
+++ b/tests/plot_data.cpp
@@ -8,7 +8,7 @@
 #include "numsim-materials/materials/drucker_prager_yield_function.h"
 #include "numsim-materials/materials/drucker_prager_plasticity.h"
 #include "numsim-materials/materials/j2_plasticity.h"
-#include "numsim-materials/solvers/backward_euler.h"
+#include "numsim-materials/solvers/local_newton.h"
 #include "numsim-materials/postprocessing/numerical_diff_checker.h"
 
 using policy = numsim::materials::material_policy_default;
@@ -68,7 +68,7 @@ run_result run_j2(T increment, int steps,
 
   p.clear();
   p.insert<std::string>("name", "solver");
-  ctx.create<numsim::materials::backward_euler<policy>>(p);
+  ctx.create<numsim::materials::local_newton<policy>>(p);
 
   p.clear();
   p.insert<std::string>("name", "hardening");
@@ -133,7 +133,7 @@ run_result run_dp(T increment, int steps,
 
   p.clear();
   p.insert<std::string>("name", "solver");
-  ctx.create<numsim::materials::backward_euler<policy>>(p);
+  ctx.create<numsim::materials::local_newton<policy>>(p);
 
   p.clear();
   p.insert<std::string>("name", "hardening");
diff --git a/tests/test_drucker_prager.cpp b/tests/test_drucker_prager.cpp
index c495c32..2f03c7d 100644
--- a/tests/test_drucker_prager.cpp
+++ b/tests/test_drucker_prager.cpp
@@ -9,8 +9,8 @@
 #include "numsim-materials/materials/linear_isotropic_hardening.h"
 #include "numsim-materials/materials/drucker_prager_yield_function.h"
 #include "numsim-materials/materials/drucker_prager_plasticity.h"
-#include "numsim-materials/materials/rk_plasticity.h"
-#include "numsim-materials/solvers/backward_euler.h"
+#include "numsim-materials/materials/j2_rk_plasticity.h"
+#include "numsim-materials/solvers/local_newton.h"
 #include "numsim-materials/solvers/butcher_tableau.h"
 #include "numsim-materials/postprocessing/numerical_diff_checker.h"
 
@@ -104,7 +104,7 @@ class DruckerPragerTest : public ::testing::Test {
 
     p.clear();
     p.insert<std::string>("name", "solver");
-    ctx.create<numsim::materials::backward_euler<policy>>(p);
+    ctx.create<numsim::materials::local_newton<policy>>(p);
 
     p.clear();
     p.insert<std::string>("name", "hardening");
@@ -209,7 +209,7 @@ class DPTangentTest : public ::testing::Test {
 
     p.clear();
     p.insert<std::string>("name", "solver");
-    ctx.create<numsim::materials::backward_euler<policy>>(p);
+    ctx.create<numsim::materials::local_newton<policy>>(p);
 
     p.clear();
     p.insert<std::string>("name", "hardening");
@@ -282,7 +282,7 @@ T run_dp_max_tangent_error(T increment, int steps) {
 
   p.clear();
   p.insert<std::string>("name", "solver");
-  ctx.create<numsim::materials::backward_euler<policy>>(p);
+  ctx.create<numsim::materials::local_newton<policy>>(p);
 
   p.clear();
   p.insert<std::string>("name", "hardening");
@@ -370,7 +370,7 @@ T max_tangent_error(std::vector<double> direction, T increment, int steps) {
 
   p.clear();
   p.insert<std::string>("name", "solver");
-  ctx.create<numsim::materials::backward_euler<policy>>(p);
+  ctx.create<numsim::materials::local_newton<policy>>(p);
 
   p.clear();
   p.insert<std::string>("name", "hardening");
@@ -467,7 +467,7 @@ TEST(DruckerPragerApex, HydrostaticTensionReachesTheApex) {
 
   p.clear();
   p.insert<std::string>("name", "solver");
-  ctx.create<numsim::materials::backward_euler<policy>>(p);
+  ctx.create<numsim::materials::local_newton<policy>>(p);
 
   p.clear();
   p.insert<std::string>("name", "hardening");
@@ -527,7 +527,7 @@ TEST(DruckerPragerApex, ApexStateIsAdmissible) {
   ctx.create<numsim::materials::linear_elasticity<policy>>(p);
   p.clear();
   p.insert<std::string>("name", "solver");
-  ctx.create<numsim::materials::backward_euler<policy>>(p);
+  ctx.create<numsim::materials::local_newton<policy>>(p);
   p.clear();
   p.insert<std::string>("name", "hardening");
   p.insert<std::string>("source", "dp");
diff --git a/tests/test_j2_plasticity.cpp b/tests/test_j2_plasticity.cpp
index a7b3429..b24f045 100644
--- a/tests/test_j2_plasticity.cpp
+++ b/tests/test_j2_plasticity.cpp
@@ -8,8 +8,8 @@
 #include "numsim-materials/materials/linear_isotropic_hardening.h"
 #include "numsim-materials/materials/drucker_prager_plasticity.h"
 #include "numsim-materials/materials/j2_plasticity.h"
-#include "numsim-materials/materials/rk_plasticity.h"
-#include "numsim-materials/solvers/backward_euler.h"
+#include "numsim-materials/materials/j2_rk_plasticity.h"
+#include "numsim-materials/solvers/local_newton.h"
 #include "numsim-materials/solvers/butcher_tableau.h"
 #include "numsim-materials/postprocessing/numerical_diff_checker.h"
 
@@ -44,7 +44,7 @@ class J2PlasticityTest : public ::testing::Test {
     // Newton-Raphson solver
     p.clear();
     p.insert<std::string>("name", "solver");
-    ctx.create<numsim::materials::backward_euler<policy>>(p);
+    ctx.create<numsim::materials::local_newton<policy>>(p);
 
     // Linear isotropic hardening (Local edge — called in inner loop)
     p.clear();
@@ -158,7 +158,7 @@ class J2TangentTest : public ::testing::Test {
 
     p.clear();
     p.insert<std::string>("name", "solver");
-    ctx.create<numsim::materials::backward_euler<policy>>(p);
+    ctx.create<numsim::materials::local_newton<policy>>(p);
 
     p.clear();
     p.insert<std::string>("name", "hardening");
diff --git a/tests/test_materials.cpp b/tests/test_materials.cpp
index 44e09cb..74f1c2c 100644
--- a/tests/test_materials.cpp
+++ b/tests/test_materials.cpp
@@ -4,6 +4,7 @@
 #include "numsim-materials/materials/tensor_component_stepper.h"
 #include "numsim-materials/materials/linear_elasticity.h"
 #include "numsim-materials/materials/scalar_stepper.h"
+#include "numsim-materials/solvers/local_newton.h"
 #include "numsim-materials/materials/scalar_identity_weight.h"
 #include "numsim-materials/materials/autocatalytic_reaction.h"
 #include "numsim-materials/solvers/backward_euler.h"
@@ -159,3 +160,67 @@ TEST(CuringSimulation, ConvergesToFullCure) {
 }
 
 } // namespace
+
+namespace {
+namespace nm_be = numsim::materials;
+
+/// backward_euler is the GRAPH-driven solver, so "function" is required.
+///
+/// It used to default to empty, which silently selected a second,
+/// callback-driven mode inside the same class: no inputs were created, update()
+/// was never bound, and "delta" stayed at zero. A consumer reading it froze at
+/// its start value for the whole analysis with no error -- measured, a cure
+/// that should reach 1.0 sat at 0.01 for 30 steps. That mode is local_newton
+/// now, chosen by naming a type rather than by omitting a parameter.
+TEST(BackwardEulerSetup, RequiresAFunctionToSolve) {
+  using policy = nm_be::material_policy_default;
+  nm_be::material_context<policy> ctx;
+  policy::ParameterHandler p;
+  p.insert<std::string>("name", "solver");
+  EXPECT_THROW(ctx.create<nm_be::backward_euler<policy>>(p),
+               std::invalid_argument)
+      << "omitting \"function\" must fail at setup, not produce an inert "
+         "solver whose consumers silently never advance";
+}
+
+/// A function material that does not publish residual/jacobian is caught at
+/// finalize, by name.
+TEST(BackwardEulerSetup, RejectsAFunctionWithoutResidualOrJacobian) {
+  using policy = nm_be::material_policy_default;
+  using T2 = policy::value_type;
+  nm_be::material_context<policy> ctx;
+  policy::ParameterHandler p;
+  p.insert<std::string>("name", "drv");
+  p.insert<T2>("increment", T2{0.1});
+  ctx.create<nm_be::scalar_stepper<policy>>(p);
+  p.clear();
+  p.insert<std::string>("name", "solver");
+  p.insert<std::string>("function", std::string("drv"));
+  ctx.create<nm_be::backward_euler<policy>>(p);
+  EXPECT_THROW(ctx.finalize(), std::runtime_error);
+}
+
+/// local_newton is the material-driven one: no function, no graph inputs, and
+/// its result carries convergence so a caller cannot read the number without
+/// the flag being at hand.
+TEST(LocalNewtonSolver, SolvesAndReportsConvergence) {
+  using policy = nm_be::material_policy_default;
+  using T2 = policy::value_type;
+  nm_be::material_context<policy> ctx;
+  policy::ParameterHandler p;
+  p.insert<std::string>("name", "solver");
+  auto& s = ctx.create<nm_be::local_newton<policy>>(p);
+  ctx.finalize();
+
+  // x^2 - 4 = 0 from x0 = 3
+  const auto ok = s.solve(
+      [](T2 x) { return std::pair<T2, T2>{x * x - T2{4}, T2{2} * x}; }, T2{3});
+  EXPECT_TRUE(ok.converged);
+  EXPECT_NEAR(ok.x, 2.0, 1e-10);
+
+  // No root: x^2 + 1 = 0. Must report failure rather than a plausible number.
+  const auto bad = s.solve(
+      [](T2 x) { return std::pair<T2, T2>{x * x + T2{1}, T2{2} * x}; }, T2{1});
+  EXPECT_FALSE(bad.converged);
+}
+}  // namespace