diff --git a/.claude/skills/benchmark/SKILL.md b/.claude/skills/benchmark/SKILL.md
index 981aabe..1d4ff4c 100644
--- a/.claude/skills/benchmark/SKILL.md
+++ b/.claude/skills/benchmark/SKILL.md
@@ -1,7 +1,7 @@
 ---
 name: benchmark
-description: 'Benchmark GIR compiled evaluator vs interpreter. Args: example number (default: 2), `--fit` for full fit, `-n N` for repetitions, `--nfev` for residual-evaluation counts, `--plan-time` for planning vs fit wall time.'
-argument-hint: '[example_num] [--fit] [-n N] [--nfev] [--plan-time]'
+description: 'Benchmark GIR compiled evaluator vs interpreter. Args: example number (default: 2), `--fit` for full fit, `-n N` for repetitions, `--nfev` for residual-evaluation counts, `--plan-time` for planning vs fit wall time, `--par-variability` for fit robustness vs perturbed initial guesses (`--starts N`).'
+argument-hint: '[example_num] [--fit] [-n N] [--nfev] [--plan-time] [--par-variability] [--starts N]'
 disable-model-invocation: true
 allowed-tools: Bash(.venv/bin/python .claude/skills/benchmark/benchmark_gir.py *)
 ---
@@ -11,6 +11,6 @@ source of truth for this skill.
 
 Pass through the user's optional benchmark arguments unchanged:
 
-- `[example_num] [--fit] [-n N] [--nfev] [--plan-time]`
+- `[example_num] [--fit] [-n N] [--nfev] [--plan-time] [--par-variability] [--starts N]`
 
 When this wrapper and the shared doc differ, follow the shared doc.
diff --git a/.claude/skills/benchmark/benchmark_gir.py b/.claude/skills/benchmark/benchmark_gir.py
index a23f127..9239341 100644
--- a/.claude/skills/benchmark/benchmark_gir.py
+++ b/.claude/skills/benchmark/benchmark_gir.py
@@ -427,6 +427,125 @@ def timed_sched(*args, **kwargs):
     return plan_total, fit_time
 
 
+# ------------------------------------------------------------------
+# Parameter variability (fit robustness vs. perturbed initial guesses)
+# ------------------------------------------------------------------
+
+
+#
+def capture_par_variability(example_num, *, n_starts=4):
+    """Run the standard pipeline from perturbed initial guesses.
+
+    Runs the baseline + fit_2d pipeline once from the example's own
+    initial values (reference) and ``n_starts`` more times with every
+    free parameter's init scaled by a fixed factor ladder, clipped
+    into bounds. Deterministic by construction -- no RNG -- so
+    repeated invocations give identical numbers.
+
+    This is a diagnostic (robustness report), not a pass/fail test:
+    a parameter whose fitted value depends on its start indicates a
+    shallow or multi-modal objective, or init-dependent machinery
+    (e.g. conv-kernel support sized from the initial value).
+    """
+
+    ladder = [0.6, 0.8, 1.25, 1.6]
+    fitted_runs = []
+    labels = []
+    redchis = []
+
+    for run in range(n_starts + 1):
+        file, dynamics_calls = load_example(example_num, add_dynamics=False)
+        file.define_baseline(
+            time_start=0, time_stop=10, time_type="ind", show_plot=False
+        )
+        file.fit_baseline(model_name="2D", stages=2, try_ci=0)
+        for call in dynamics_calls:
+            file.add_time_dependence(**call)
+
+        model = file.model_active
+        assert model is not None
+        free_names = [
+            name for name, par in model.lmfit_pars.items() if par.vary and not par.expr
+        ]
+
+        label = "reference" if run == 0 else f"start {run}"
+        if run > 0:
+            for i, name in enumerate(free_names):
+                par = model.lmfit_pars[name]
+                factor = ladder[(i + run - 1) % len(ladder)]
+                value = par.value * factor
+                if abs(par.value) < 1e-12:
+                    # multiplicative perturbation is a no-op at 0; nudge
+                    # by a fraction of the bound range when available
+                    if np.isfinite(par.min) and np.isfinite(par.max):
+                        value = par.value + 0.1 * (par.max - par.min)
+                    else:
+                        continue
+                if np.isfinite(par.min) and np.isfinite(par.max):
+                    span = par.max - par.min
+                    value = min(
+                        max(value, par.min + 0.01 * span),
+                        par.max - 0.01 * span,
+                    )
+                elif value <= par.min or value >= par.max:
+                    continue  # would leave bounds; keep original init
+                par.value = value
+
+        t0 = time.perf_counter()
+        file.fit_2d(model_name="2D", stages=2, try_ci=0)
+        wall = time.perf_counter() - t0
+
+        df = file.get_fit_results(fit_type="2d")
+        fitted = dict(zip(df["name"], df["value"], strict=True))
+        fitted_runs.append({name: fitted[name] for name in free_names})
+
+        try:
+            redchi = model.result[1].redchi
+        except (AttributeError, IndexError, TypeError):
+            redchi = float("nan")
+        labels.append(label)
+        redchis.append(redchi)
+        print(f"  {label:12s} redchi={redchi:10.4g}  wall={wall:6.2f} s")
+
+    # Separate the two failure signals: a run that converged to a worse
+    # optimum (secondary minimum) vs. spread among runs that reached the
+    # best optimum (identifiability / init-dependent machinery).
+    redchi_arr = np.array(redchis)
+    best_redchi = np.nanmin(redchi_arr)
+    ok_idx = [
+        i
+        for i in range(len(redchis))
+        if np.isfinite(redchi_arr[i]) and redchi_arr[i] <= best_redchi * 1.01
+    ]
+    off_idx = [i for i in range(len(redchis)) if i not in ok_idx]
+
+    print()
+    if off_idx:
+        off_str = ", ".join(f"{labels[i]} (redchi {redchis[i]:.4g})" for i in off_idx)
+        print(f"  off-optimum runs (excluded from spread): {off_str}")
+    print(f"  {'parameter':32s}{'reference':>12s}{'spread':>12s}{'rel':>8s}")
+    worst_rel = 0.0
+    worst_name = "-"
+    for name in fitted_runs[0]:
+        values = np.array([fitted_runs[i][name] for i in ok_idx])
+        spread = float(values.max() - values.min())
+        scale = max(abs(float(values.mean())), 1e-300)
+        rel = spread / scale
+        if rel > worst_rel:
+            worst_rel, worst_name = rel, name
+        flag = "  <-- start-sensitive" if rel > 0.01 else ""
+        print(
+            f"  {name:32s}{fitted_runs[0][name]:12.6g}{spread:12.3g}"
+            f"{rel * 100:7.2f}%{flag}"
+        )
+    print()
+    print(
+        f"  {len(ok_idx)}/{len(redchis)} runs at best optimum; "
+        f"worst spread among them: {worst_name} ({worst_rel * 100:.2f}%)"
+    )
+    return worst_rel, worst_name, len(off_idx)
+
+
 # ------------------------------------------------------------------
 # Full-fit benchmark
 # ------------------------------------------------------------------
@@ -536,11 +655,29 @@ def bench_fit(example_num, dynamics_calls, *, n_reps=3):
             "wall time. If --example is 0, runs all examples."
         ),
     )
+    parser.add_argument(
+        "--par-variability",
+        action="store_true",
+        help=(
+            "Fit from deterministically perturbed initial guesses and "
+            "report the spread of fitted values across starts (no RNG, "
+            "reproducible). --starts controls the number of perturbed "
+            "starts. If --example is 0, runs all examples."
+        ),
+    )
+    parser.add_argument(
+        "--starts",
+        type=int,
+        default=4,
+        help="Perturbed starts for --par-variability (default: 4)",
+    )
     args = parser.parse_args()
 
-    # --nfev / --plan-time with example 0 means "all examples"; skip the
+    # Multi-example modes with example 0 mean "all examples"; skip the
     # single-example preamble so the capture fn can load each one itself.
-    skip_preamble = (args.nfev or args.plan_time) and args.example == 0
+    skip_preamble = (
+        args.nfev or args.plan_time or args.par_variability
+    ) and args.example == 0
     if not skip_preamble:
         folder = _find_example_folder(args.example)
         print(f"Example: {folder.name}")
@@ -615,6 +752,36 @@ def bench_fit(example_num, dynamics_calls, *, n_reps=3):
             print(f"PLAN-TIME CAPTURE -- example {args.example:02d}")
             print("=" * 60)
             capture_plan_time(args.example)
+    elif args.par_variability:
+        if args.example == 0:
+            summaries: dict[int, str] = {}
+            for n in range(1, 5):
+                print()
+                print("=" * 60)
+                print(f"PARAMETER VARIABILITY -- example {n:02d}")
+                print("=" * 60)
+                try:
+                    worst_rel, worst_name, n_off = capture_par_variability(
+                        n, n_starts=args.starts
+                    )
+                    summaries[n] = (
+                        f"worst {worst_rel * 100:6.2f}%  ({worst_name})"
+                        f"{f'  [{n_off} off-optimum]' if n_off else ''}"
+                    )
+                except FileNotFoundError as e:
+                    print(f"  skipped: {e}")
+                    summaries[n] = "skipped"
+            print()
+            print("=" * 60)
+            print("PARAMETER VARIABILITY SUMMARY")
+            for n, summary_str in summaries.items():
+                print(f"  example {n:02d}:  {summary_str}")
+            print("=" * 60)
+        else:
+            print("=" * 60)
+            print(f"PARAMETER VARIABILITY -- example {args.example:02d}")
+            print("=" * 60)
+            capture_par_variability(args.example, n_starts=args.starts)
     else:
         bench_per_call(file, n_calls=args.calls)
 
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 48a85c5..47792b7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,24 @@ and this project adheres to [Semantic Versioning](https://semver.org/).
 This file is maintained using the shared changelog workflow in
 [`docs/ai/changelog.md`](docs/ai/changelog.md).
 
+## [0.10.0] - 2026-07-08
+
+### Added
+
+- **`stepFun(t, A, t0)` dynamics primitive**: a causal step (0 before `t0`, `A` after). Since dynamics components combine by addition, `stepFun` is the dedicated way to model baselines and plateaus (e.g. `expFun` + `stepFun` sharing `t0` gives a decay to a nonzero plateau); `erfFun` is its Gaussian-broadened counterpart. Registered on both the mcp and compiled (GIR) evaluation paths.
+
+### Changed
+
+- **Breaking: `y0` removed from all dynamics functions** (`linFun`, `expFun`, `sinFun`, `sinDivX`, `erfFun`, `sqrtFun`). The per-function offset broke the causality convention `f(t < t0) = 0` inconsistently (`sqrtFun`/`erfFun` leaked `y0` before `t0`; the others clamped to 0) and duplicated what an additive component expresses directly. Migration for model YAMLs: delete `y0:` lines that were `0` (the common case); replace a nonzero `y0` with a `stepFun` component sharing the function's `t0`. Stale `y0:` entries now fail model validation with a parameter-count error. The energy-domain `Offset(x, y0)` background is unaffected.
+- **Breaking: `Voigt` amplitude is now grid-independent.** The profile was normalized by its maximum over the sampled energy window, so the value at a fixed energy depended on the window/grid, and a peak center pushed outside the fit window rescaled the in-window tail up to amplitude `A`. Normalization now uses the analytic peak value `wofz(i·W/(2·SD·√2))` at `dx = 0`; `A` is the true peak height regardless of the sampled grid. Fitted `A` values may shift slightly relative to earlier releases (the two normalizations agree only when the sampled grid contains the exact peak). `voigtCONV` is unchanged (its kernel axis is symmetric around the peak, where max-over-grid is exact).
+
+### Fixed
+
+- **`expRiseCONV` docstring claimed the kernel is causal.** It is deliberately anti-causal — the mirror of `expDecayCONV`, nonzero only for `x ≤ 0`, so the convolved response rises before the excitation and saturates at `t0`. The docstring now says so.
+- **`boxCONV` docstring claimed "smooth edges"** — the kernel is a hard `|x| ≤ width/2` threshold; docstring corrected.
+- **`GLS`/`GLP` docstrings** now state the valid mixing range `m ∈ [0, 1]` (for `GLP`, `m < 0` makes the denominator `1 + 4·m·u²` cross zero, producing NaN/inf).
+- **`DS` docstring** now documents that the asymmetric tail extends toward `x < x0` (correct for a kinetic-energy axis; mirrored on a binding-energy axis).
+
 ## [0.9.3] - 2026-06-29
 
 ### Added
diff --git a/CLAUDE.md b/CLAUDE.md
index eac8197..2aa7aa6 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -5,6 +5,7 @@
 - **Token Efficiency:** Be concise. Reference file paths and line numbers rather than quoting large code blocks.
 - **Subagent Protocol:** Use subagents for repo-wide scans, parallel research, or scanning large directories. Instruct them to return only concise summaries to keep the main context window lean.
 - **Guardrails:** Use this file for coding guardrails; heavier review checklists live in `docs/ai/code-review.md`.
+- **Task recipes & skills:** Common repo tasks have step-by-step recipes in `docs/ai/*.md` (add-function, changelog, check-docs, check-example, benchmark, bump-versions, code-review), each mirrored by a thin `/`-invokable wrapper in `.claude/skills/` whose source of truth is the `docs/ai` file. Most set `disable-model-invocation: true`, so they will NOT appear in an agent's auto-invoke list—before hand-rolling one of these tasks, consult the matching `docs/ai` recipe directly.
 - **Renaming / API changes:** grep the entire repo—notebooks, YAML, tests, and docs all reference the public API.
 
 
diff --git a/TODO.md b/TODO.md
index 806b3c4..e6c3c22 100644
--- a/TODO.md
+++ b/TODO.md
@@ -2,15 +2,15 @@
 
 ## Fitting
 
-- [ ] **Mismatched initial guesses**: round-trip tests — one each for basic, profile, profile+dynamics.
-- [ ] **Conv kernel support is sized from the initial parameter value**: `Component.create_t_kernel` builds the kernel time axis once at model construction (`t_range = par_init * kernel_width`, e.g. ±4·SD for `gaussCONV`) and never rebuilds it. When the fitted width grows past its init, the kernel is silently truncated, biasing the recovered width (observed while building `04_parameter_profiles`, 2026-06-11: truth SD=10, init 5 → fitted SD≈10.8 at any SNR; 03's SD 0.148 vs truth 0.15 with init 0.1 is likely the same effect, mild). The GIR path snapshots the same static axis (`graph_ir.py`, `kernel_time`). Fix candidates: rebuild the kernel axis when the kernel parameter value changes; size the support from the parameter's max bound; or at minimum warn when `fitted_value * kernel_width` exceeds the kernel range. Workaround used in the examples: initialize conv widths generously above the expected value (commented in the YAMLs).
+- [ ] **Enable 1D time-trace fitting (post-SbS kinetics)**: wire standalone `TIME_1D` dynamics models (already evaluable on the mcp path per `docs/design/supported_models.md`, but connected to no fit method) to a fit entry point, so parameter-vs-time traces extracted from SbS results can be fit to `functions/time.py` dynamics (`expFun` sums, IRF convolution) inside the package instead of in external scripts. Frame it as the diagnostic/initialization rung before `fit_2d`, not a statistical equivalent (per-slice correlations and unpropagated SbS uncertainties make two-step inferior). Design direction: promote SbS results into a time-axis `File` — the promoted object inherits limits/CI/MCMC/save/export machinery, and trace fits land in a `SavedFitSlot` like every other fit type (no parallel fitting pipeline). Requires `File` to support time as the primary axis, which today it assumes is energy. Weighting traces by per-slice stderr ties into the `sigma_type` expansion item below.
+- [ ] **Conv kernel support is sized from the initial parameter value**: `Component.create_t_kernel` builds the kernel time axis once at model construction (`t_range = par_init * kernel_width`, e.g. ±4·SD for `gaussCONV`) and never rebuilds it. When the fitted width grows past its init, the kernel is silently truncated, biasing the recovered width (observed while building `04_parameter_profiles`, 2026-06-11: truth SD=10, init 5 → fitted SD≈10.8 at any SNR; 03's SD 0.148 vs truth 0.15 with init 0.1 is likely the same effect, mild). The GIR path snapshots the same static axis (`graph_ir.py`, `kernel_time`). Fix candidates: rebuild the kernel axis when the kernel parameter value changes; size the support from the parameter's max bound; or at minimum warn when `fitted_value * kernel_width` exceeds the kernel range. Workaround used in the examples: initialize conv widths generously above the expected value (commented in the YAMLs). When fixing, add a deterministic regression test that inits the conv width well below truth (the failure mode above); broad robustness-vs-start checking deliberately lives in the benchmark skill (`--par-variability`, added 2026-07-06), not the test suite, to avoid weak-or-flaky convergence asserts.
 
 
 ## Noise and simulation
 
 - [ ] **Simulator noise-language cleanup**: align simulator docs/metadata with the fit-results noise schema. Keep simulator `noise_type` meaning "noise distribution / random generator" (`gaussian`, `poisson`, `none`), not sigma shape. Fix the stale `Simulator.set_noise_type()` docstring that mentions `uniform`; clarify `detection` vs. `noise_type` vs. `noise_level`; and, for analog Gaussian simulations, consider saving the derived `sigma_data = noise_level * max(abs(clean_data))` alongside existing metadata. For parameter sweeps, store derived `sigma_data` per configuration when it depends on each clean dataset.
 - [ ] **Warn on subcycle-boundary time samples**: in `Dynamics.normalize_time`, emit a warning when a time sample lands within epsilon of a subcycle boundary (`|t * f_eff - round(t * f_eff)| < eps`). Subcycle masks switch discretely there, so the sample's assignment — and hence the model prediction for that whole row — is sensitive to floating-point representation of the time axis. This silently biased a multi-cycle fit when synthetic data was generated on `np.arange` axes but fit against their `%.6e`-rounded CSV reload (see `03_multi_cycle_dynamics/data/generate_data.ipynb`, 2026-06-11). A warning turns that silent bias into a one-line diagnosis; it also flags the physically ambiguous case of measuring exactly at the switching instant.
-- [ ] **Future `sigma_type` expansion in FitResults**: after the first constant, user-supplied sigma schema lands, extend uncertainty handling beyond scalar `sigma_data`. Keep `noise_type` for the statistical assumption/distribution and use `sigma_type` for sigma shape: initially `constant`, later `per_spectrum` and `per_point`. Add HDF5 storage, validation, baseline/SBS/2D alignment, `compare_models()` behavior, and tests for vector/matrix sigma. Defer automatic Poisson-derived sigma until residual-space variance propagation is explicit.
+- [ ] **Future `sigma_type` expansion in FitResults**: the constant, user-supplied sigma schema has landed (`SIGMA_TYPE_CONSTANT` in [fit_io.py](src/trspecfit/utils/fit_io.py) ~L54; `validate_noise_metadata` hard-locks `sigma_type` to `"constant"`), so this is now unblocked: extend uncertainty handling beyond scalar `sigma_data`. Keep `noise_type` for the statistical assumption/distribution and use `sigma_type` for sigma shape: initially `constant`, later `per_spectrum` and `per_point`. Add HDF5 storage, validation, baseline/SBS/2D alignment, `compare_models()` behavior, and tests for vector/matrix sigma. Defer automatic Poisson-derived sigma until residual-space variance propagation is explicit.
 
 ## Performance & architecture
 
@@ -21,28 +21,31 @@
   - **Persist `correl` and `acceptance_fraction` into the slots**: 2026-06 added live-only accessors (`get_correlations`, `get_conf_intervals`, `get_mcmc`) reading `model.result` as a stopgap for notebook 12, so these are NOT yet saved. Add per-parameter correlations to the slot `params` payload and `acceptance_fraction` to the slot `mcmc` payload, with `.fit.h5` read/write support and save/load round-trip tests, so they survive persistence like the rest of the slot.
   - **Relocate the live accessors to `FitResults`**: 2026-06 added `File.get_correlations`, `File.get_conf_intervals`, `File.get_mcmc` (and the private `File._result_model` resolver) reading `model.result[...]` directly. These conceptually belong on `FitResults` (like `compare_models`, which already lives there with `File.compare_models` as thin sugar). The existing `File.get_fit_results` is in the same boat. Decide whether all of these should move into `FitResults` (with thin `File.*` sugar that delegates), and whether they read live `model.result` or persisted slots — then move them and update callers (notebook 12 reads them).
   - The raw list-index access (`result[1..4]`) and the deeper unified-results-object question are deferred to this item.
-- [ ] **Disentangle plotting from saving/conversion in the fit pipeline**: figure rendering is currently entangled with data conversion and file IO. `fitlib.results_to_df` (results → DataFrame) and `File._save_2d_fit_legacy` / `_save_sbs_fit_legacy` (CSV writers) also render figures, and `fit_baseline` / `fit_slice_by_slice` / `fit_2d` reach plotting only by calling that save-legacy path. Split into (a) compute/convert, (b) explicit save/export, (c) an explicit plotting API — cf. the `_save_img_flag` helper and `FitResults.plot_residuals`. Until then, the fit methods gate display on `show_output` and saving on `auto_export` via a `save_files` flag through the legacy methods. Scope: beyond the examples-upgrade branch.
+- [ ] **Disentangle plotting from saving/conversion in the fit pipeline**: figure rendering is currently entangled with data conversion and file IO. `fitlib.results_to_df` (results → DataFrame) and `File._save_2d_fit_legacy` / `_save_sbs_fit_legacy` (CSV writers) also render figures, and `fit_slice_by_slice` / `fit_2d` reach plotting only by calling that save-legacy path (`fit_baseline` is already disentangled — it calls `fitlib.plt_fit_res_1d` directly, gated by `_save_img_flag`; use it as the template). Blocks the legacy-shim removal item in the v1.0.0 checklist, since the `_save_*_legacy` impls are the live SbS/2D plotting path. Split into (a) compute/convert, (b) explicit save/export, (c) an explicit plotting API — cf. the `_save_img_flag` helper and `FitResults.plot_residuals`. Until then, the fit methods gate display on `show_output` and saving on `auto_export` via a `save_files` flag through the legacy methods. Scope: beyond the examples-upgrade branch.
 - [ ] **Decouple directory creation from path computation (`create_model_path`)**: `Project.create_model_path` ([trspecfit.py](src/trspecfit/trspecfit.py) ~L2167) `mkdir`s the `{path_results}/{file}/{fit_type}/{model}/` tree as a side effect of *computing* the path. Every fit method calls it to build the `save_path` handed to `fit_wrapper` (`fit_baseline` ~L2572, `fit_2d` ~L3983, spectrum ~L2804, sbs ~L3120), and the project plot-grid display builds image paths through it too (`fit_baselines` ~L1070, `fit_2d` ~L1432) — so empty `*_fits/` dir trees appear even when `auto_export: False` and nothing is written (observed in 21, 2026-06-27; every fitting example already does this). `fit_wrapper` only writes when `save_output=1 if auto_export else 0`. Fix: make `create_model_path` return the path *without* `mkdir`, and `mkdir`-on-write at the actual write sites (`fit_wrapper`'s save branch + the explicit `save_*`/`export_*` functions); the plot-grid readers then just build a path and `.exists()`-check. Do **not** gate on `auto_export` directly — explicit `save_fit()`/`export_fit()` legitimately need the dir and aren't `auto_export`-driven. Cosmetic only (the dirs are gitignored, no functional bug); blast radius spans baseline/spectrum/sbs/2d + savers, needs no-regression tests. Related to the plotting/saving disentanglement item above; scope: beyond the examples-upgrade branch.
 
 ## Testing
 
-- [ ] **Pylance/Pyright `| None` noise in tests**: ~280 Pyright errors across tests, all from accessing `File`/`Model` attributes typed as `ndarray | None`. Current `# type guard` asserts are inconsistent and don't propagate through helper methods. Find a cleaner pattern (e.g. `TypeGuard`, narrowing wrapper, or Pyright config) and apply consistently.
+- [ ] **Pylance/Pyright `| None` noise in tests**: several hundred Pyright errors across tests, mostly from accessing `File`/`Model` attributes typed as `ndarray | None` (`[tool.pyright]` only includes `src`, so these surface in editors, not CI). Current `# type guard` asserts (~110) are inconsistent and don't propagate through helper methods. Find a cleaner pattern (e.g. `TypeGuard`, narrowing wrapper, or Pyright config) and apply consistently.
 
 ## User and AI ergonomics
 
-- [ ] **Curate the public API surface before v1.0.0**: audit user-facing classes (`File`, `Project`, `Simulator`, `Model`, etc.) and decide which methods/attributes should be discoverable in notebooks and docs. Add curated `__dir__()` output for autocomplete, keep `__all__`/API docs aligned, and gradually rename or deprecate internal helper methods that should not look like primary user workflows. This should improve both human notebook ergonomics and AI/LLM efficiency by making the intended workflow surface smaller, clearer, and easier to infer.
-- [ ] **Document API tiers**: add a short guide that separates stable user API (`Project`, `File`, `Simulator`, `PlotConfig`), advanced public API (`mcp.Model`, `Component`, `Par`, `ParameterSweep`, `MC`), and internal implementation modules (`graph_ir`, `eval_1d`, `eval_2d`, low-level parsing/HDF5 helpers). Use this as the source of truth for docs, tests, examples, and AI-agent guidance.
 - [ ] **Add tool-neutral agent orientation**: add `AGENTS.md` or `docs/ai/agent-orientation.md` pointing agents to `CLAUDE.md`, `TODO.md`, `PLAN.md`, `docs/design/repo_architecture.md`, supported-model docs, common commands, and API-change guardrails. Keep it concise so any LLM can quickly find the intended workflow and repo boundaries.
-- [ ] **Add more AI-friendly task recipes**: extend `docs/ai/` with checklists for common repo changes, such as adding YAML syntax, adding plotting options, changing fitting workflows, modifying GIR/evaluator behavior, extending save/load fields, and preparing a release.
-- [ ] **`trspec.fit` project website**: stand up the site as the home for fully rendered example notebooks (plots, fit tables) so users can browse without installing — preferred over Read the Docs, which has build timeouts/execution limits for notebooks. Room to grow into interactive tooling: a model-builder UI (GUI that outputs the YAML model files) and other future work. Keep git notebook sources stripped (`nbstripout`); render/execute at publish time so outputs never drift from the code.
+- [ ] **Add more AI-friendly task recipes**: `docs/ai/` already covers add-function, check-example, check-docs, code-review, changelog, and benchmark. Extend with checklists for the remaining common changes: adding YAML syntax, adding plotting options, changing fitting workflows, modifying GIR/evaluator behavior, extending save/load fields, and preparing a release.
 - [ ] **Revisit default `Project.name` for saving examples**: `Project.name` defaults to `"test"` ([trspecfit.py](src/trspecfit/trspecfit.py) ~L202), so a bare `save_fit()` / `save_fits()` with no path writes `./fit_results/test.fit.h5`. The example `project.yaml` files set no `name:`. When working through the saving notebooks (`11_save_load_export`, `20_multi_file_independent_fit`), decide on a convention — set a meaningful `name:` per example so the default archive path is self-describing, and/or keep using explicit content-named paths (`comparison.fit.h5`, `batch.fit.h5`). Pick one and apply consistently.
 - [ ] **Add minimal runnable workflow examples**: supplement notebooks with small script-like examples or docs snippets for the canonical public workflows: load data, load a model, set limits, fit baseline, fit 2D, inspect results, simulate data, and run a parameter sweep.
-- [ ] **Improve public validation errors**: make user-facing errors state what failed, where it failed (file/model/component/parameter when applicable), and what the user or agent should change next. Prioritize YAML parsing, model loading, fit setup, and unsupported-model fallback paths.
-- [ ] **Tighten public type hints and aliases**: reduce ambiguous `Any` on public APIs, document key aliases such as `ModelRef`, and keep return types crisp for IDEs, Pyright, generated docs, and LLM code navigation.
+- [ ] **Improve public validation errors**: YAML parsing and model loading already meet the bar — `ModelValidationError` messages carry model/component/parameter context plus hints (verified 2026-07-06). Remaining scope: audit fit setup and unsupported-model fallback paths and bring them to the same standard (state what failed, where, and what to change next).
 - [ ] **Document fast verification slices**: add focused pytest commands for common edits (public workflow/API, YAML parser, functions, GIR/evaluator, plotting) so contributors and agents can validate changes quickly before running the full suite.
 
-## Build & release
+## Road to v1.0.0 & first adopters
 
-- [ ] **Automate tagging and pushing**: automate `git tag v1.2.3` + `git push v1.2.3` as part of the release workflow.
-- [ ] **Remove legacy/backwards-compat code**: before v1.0.0 release, audit codebase for legacy fallbacks and backwards compatibility shims and consider removing. Known shims slated for removal:
-  - `File.save_sbs_fit` / `File.save_2d_fit` (deprecated wrappers — replace internal callers and drop the public methods plus their `_save_*_fit_legacy` impls; users should migrate to `File.export_fit(fit_type=...)`).
+Rationale (2026-07-06): lock-in comes from schemas and public names, not from users — so stabilize the contract surfaces (1–5) before the website (6) and outreach (tracked outside the repo: `~/Desktop/trspecfit_outreach.md`). Early TR-XPS/TR-XAS adopters are then a feature, not a risk: they supply the usage signal that feature decisions currently lack, and their needs are a strict subset of the general "parameterized lineshapes vs. control axis" problem, so serving them does not distort the architecture.
+
+1. [ ] **Version the model YAML schema**: model YAML files are the artifact user groups accumulate, and currently the only unversioned contract — the fit archive already carries `SCHEMA_VERSION` ([fit_io.py](src/trspecfit/utils/fit_io.py) ~L46) and refuses mismatched reads/appends. Add an optional top-level format-version key in `utils/parsing.py` (absent = current version) so future syntax changes can warn or migrate instead of silently misparsing old model files.
+2. [ ] **Curate the public API surface**: audit user-facing classes (`File`, `Project`, `Simulator`, `Model`, etc.) and decide which methods/attributes should be discoverable in notebooks and docs. Add curated `__dir__()` output for autocomplete, keep `__all__`/API docs aligned, and gradually rename or deprecate internal helper methods that should not look like primary user workflows. This should improve both human notebook ergonomics and AI/LLM efficiency by making the intended workflow surface smaller, clearer, and easier to infer. Prerequisite for outreach: once external groups have notebooks, every public-looking name is frozen in practice.
+3. [ ] **Document API tiers**: add a short guide that separates stable user API (`Project`, `File`, `Simulator`, `PlotConfig`), advanced public API (`mcp.Model`, `Component`, `Par`, `ParameterSweep`, `MC`), and internal implementation modules (`graph_ir`, `eval_1d`, `eval_2d`, low-level parsing/HDF5 helpers). Use this as the source of truth for docs, tests, examples, and AI-agent guidance.
+4. [ ] **Remove legacy/backwards-compat code**: audit codebase for legacy fallbacks and backwards-compatibility shims and consider removing before v1.0.0. Full removal depends on the plotting/saving disentanglement item (Performance & architecture) — the `_save_*_fit_legacy` impls are still the live SbS/2D plotting path. Known shims slated for removal:
+   - `File.save_sbs_fit` / `File.save_2d_fit` (deprecated wrappers — replace internal callers and drop the public methods plus their `_save_*_fit_legacy` impls; users should migrate to `File.export_fit(fit_type=...)`).
+   - Legacy pre-rename format branch in [sweep.py](src/trspecfit/utils/sweep.py) ~L624 and the legacy `save_img` int-mapping helper in [plot.py](src/trspecfit/utils/plot.py) ~L813.
+5. [ ] **State a pre-1.0 stability & deprecation policy**: short docs note establishing the social contract for early adopters — pre-1.0 public API may change, but only with a `DeprecationWarning` and a changelog entry; versioned schemas (`.fit.h5` archive, model YAML) refuse or migrate rather than misread. Cheap to write, and it is what makes early adoption safe for both sides.
+6. [ ] **`trspec.fit` project website**: stand up the site as the home for fully rendered example notebooks (plots, fit tables) so users can browse without installing — preferred over Read the Docs, which has build timeouts/execution limits for notebooks. Room to grow into interactive tooling: a model-builder UI (GUI that outputs the YAML model files) and other future work. Keep git notebook sources stripped (`nbstripout`); render/execute at publish time so outputs never drift from the code.
diff --git a/docs/ai/add-function.md b/docs/ai/add-function.md
index a058f81..30520f6 100644
--- a/docs/ai/add-function.md
+++ b/docs/ai/add-function.md
@@ -132,7 +132,7 @@ Time dynamics functions:
       dispatch.
 - [ ] Verify `can_lower_2d()` accepts the new function where intended.
 - [ ] If lowering fails, recommend a signature that matches the compiled
-      dynamics pattern (typically `func(t, par1, ..., t0, y0)` with vectorized
+      dynamics pattern (typically `func(t, par1, ..., t0)` with vectorized
       NumPy math and no extra runtime context).
 
 Convolution kernels:
@@ -159,6 +159,28 @@ Profile functions:
       current profile model (`func(x, par1, ...)` returning broadcast-friendly
       NumPy arrays).
 
+### GIR parity test (required when GIR support lands)
+
+The unit tests in §3 call the function directly, so they do NOT exercise the
+compiled enum/dispatch/param-count wiring. A wrong `DYNAMICS_DISPATCH` /
+`CONV_KERNEL_DISPATCH` param count, a missing enum entry, or a dispatch typo
+passes every unit and signature-discovery test and only breaks when the
+function is actually lowered. So when the function works on GIR:
+
+- [ ] Add an end-to-end parity case that drives the function through the
+      compiled path and compares against the MCP interpreter, via the shared
+      `_compare_evaluator_vs_interpreter` + `_perturb_theta` helpers. Use
+      `tests/test_evaluate_1d.py` and/or `tests/test_evaluate_2d.py` as
+      applicable, matching whichever of `can_lower_1d()` / `can_lower_2d()` the
+      function targets above. Dynamics drivers and convolution kernels have
+      parametrized lists to extend in the 2D file (`_DYN_DRIVER_PARITY_CASES`,
+      `_KERNEL_PARITY_CASES`); energy/profile functions have per-function parity
+      tests to mirror.
+- [ ] Where the function has several parameters, perturb more than one so a
+      parameter-order swap can't slip through. Watch for init values that mask
+      an ordering bug: if two args are both 0 at init, set one nonzero so a
+      swap actually changes the output.
+
 Before finishing, explicitly state one of:
 
 - `Function is available on MCP only; GIR falls back intentionally.`
diff --git a/docs/ai/benchmark.md b/docs/ai/benchmark.md
index 1378bef..1d2163c 100644
--- a/docs/ai/benchmark.md
+++ b/docs/ai/benchmark.md
@@ -41,6 +41,8 @@ Parse the arguments:
 - First positional integer -> `--example N` (default: `2`)
 - `--fit` -> include full-fit benchmark
 - `-n N` -> fit repetitions (default: `3`)
+- `--par-variability` -> fit-robustness report (see below); `--starts N` sets
+  the number of perturbed starts (default: `4`)
 
 Run:
 
@@ -72,6 +74,35 @@ table at the end.
 .venv/bin/python .claude/skills/benchmark/benchmark_gir.py --example <N> --plan-time
 ```
 
+## Parameter variability (fit robustness vs. initial guesses)
+
+`--par-variability` runs the baseline + `fit_2d` pipeline once from the
+example's own initial values (reference) and `--starts` more times (default 4)
+with every free parameter's init scaled by a fixed factor ladder
+(x0.6 / x0.8 / x1.25 / x1.6, clipped into bounds). No RNG — repeated
+invocations give identical numbers.
+
+The report separates the two failure signals:
+
+- **off-optimum runs** — a start that converged to a worse redchi (secondary
+  minimum). Excluded from the spread statistics and listed separately.
+- **spread among best-optimum runs** — fitted-value spread across starts that
+  reached the same optimum. Nonzero spread here indicates a flat objective
+  direction (the parameter is not identifiable from the data) or
+  init-dependent machinery: state derived once from initial parameter values
+  and never rebuilt during the fit, such as a convolution-kernel support.
+  Parameters above 1% relative spread are flagged `start-sensitive`.
+
+This is a diagnostic, not a pass/fail test — it deliberately lives here rather
+than in the pytest suite, where perturbed-start convergence asserts would be
+either weak (tiny perturbations) or flaky (aggressive ones). Accepts
+`--example 0` for an all-examples summary.
+
+```bash
+.venv/bin/python .claude/skills/benchmark/benchmark_gir.py --example <N> --par-variability
+.venv/bin/python .claude/skills/benchmark/benchmark_gir.py --example 0 --par-variability --starts 6
+```
+
 ## Profiling (GIR path only)
 
 For flamegraphs of the GIR hot path, use `--profile` to run a GIR-only loop
diff --git a/docs/design/lowered_evaluator.md b/docs/design/lowered_evaluator.md
index aac661e..db49784 100644
--- a/docs/design/lowered_evaluator.md
+++ b/docs/design/lowered_evaluator.md
@@ -324,25 +324,23 @@ Only Shirley has a true data dependency on the peak sum.
 #### 1.5 Time-dependent model as a graph
 
 When `GLP_01_A` has a dynamics model (`expFun` with params `A`, `tau`,
-`t0`, `y0`):
+`t0`):
 
 ```
 Additional nodes:
  20: OPT_PARAM     "GLP_01_A_expFun_01_A"     value=5.0
  21: OPT_PARAM     "GLP_01_A_expFun_01_tau"   value=100.0
  22: OPT_PARAM     "GLP_01_A_expFun_01_t0"    value=0.0
- 23: STATIC_PARAM   "GLP_01_A_expFun_01_y0"   value=0.0
- 24: DYNAMICS_TRACE "GLP_01_A_dynamics"         function_name="expFun"  package="time"
- 25: PARAM_PLUS_TRACE "GLP_01_A_resolved"
+ 23: DYNAMICS_TRACE "GLP_01_A_dynamics"         function_name="expFun"  package="time"
+ 24: PARAM_PLUS_TRACE "GLP_01_A_resolved"
 
 Additional edges:
- 20 -> 24  PARAM_INPUT(pos=0)    # A -> dynamics
- 21 -> 24  PARAM_INPUT(pos=1)    # tau -> dynamics
- 22 -> 24  PARAM_INPUT(pos=2)    # t0 -> dynamics
- 23 -> 24  PARAM_INPUT(pos=3)    # y0 -> dynamics
-  2 -> 25  BASE_INPUT            # GLP_01_A base value
- 24 -> 25  TRACE_INPUT           # dynamics trace
- 25 -> 10  PARAM_INPUT(pos=0)    # resolved A -> GLP_01 (replaces edge 2->10)
+ 20 -> 23  PARAM_INPUT(pos=0)    # A -> dynamics
+ 21 -> 23  PARAM_INPUT(pos=1)    # tau -> dynamics
+ 22 -> 23  PARAM_INPUT(pos=2)    # t0 -> dynamics
+  2 -> 24  BASE_INPUT            # GLP_01_A base value
+ 23 -> 24  TRACE_INPUT           # dynamics trace
+ 24 -> 10  PARAM_INPUT(pos=0)    # resolved A -> GLP_01 (replaces edge 2->10)
 ```
 
 #### 1.6 Expression model as a graph
@@ -870,9 +868,9 @@ from `functions/time.py`:
 
 ```python
 DYNAMICS_DISPATCH = {
-    0: fcts_time.expFun,      # (t, A, tau, t0, y0) -> (n_time,)
-    1: fcts_time.sinFun,      # (t, A, f, phi, t0, y0) -> (n_time,)
-    2: fcts_time.linFun,      # (t, m, t0, y0) -> (n_time,)
+    0: fcts_time.expFun,      # (t, A, tau, t0) -> (n_time,)
+    1: fcts_time.sinFun,      # (t, A, f, phi, t0) -> (n_time,)
+    2: fcts_time.linFun,      # (t, m, t0) -> (n_time,)
     ...
 }
 ```
@@ -1030,8 +1028,9 @@ This replaces `2 * n_free_params + 1` evaluator calls per iteration with
 1 evaluator call + 1 Jacobian call. For 4 free params, that's 9 -> 2
 calls, ~4.5x fewer evaluations per iteration.
 
-**Variable projection (VARPRO):** Linear parameters (amplitudes `A`,
-offset `y0`, slope `m`) can be solved in closed form given the nonlinear
+**Variable projection (VARPRO):** Linear parameters (amplitudes `A`, the
+`Offset` background `y0`, slope `m`) can be solved in closed form given
+the nonlinear
 parameters. Reduces optimizer dimensionality. The graph makes identifying
 linear params straightforward: any param that appears as a linear factor
 in its component's function.
diff --git a/docs/design/repo_architecture.md b/docs/design/repo_architecture.md
index 8a12509..b7e34d1 100644
--- a/docs/design/repo_architecture.md
+++ b/docs/design/repo_architecture.md
@@ -212,8 +212,10 @@ current peak sum (e.g. Shirley). Add new peak or background shapes here.
 ### `functions/time.py`
 
 Dynamics and convolution kernels. Dynamics functions (e.g. `expFun`,
-`sinFun`, `linFun`, `erfFun`, `sqrtFun`) share signature
-`func(t, par1, ..., t0, y0)` with the invariant `f(t < t0) = 0`.
+`sinFun`, `linFun`, `erfFun`, `sqrtFun`, `stepFun`) share signature
+`func(t, par1, ..., t0)` with the invariant `f(t < t0) = 0`; constant
+offsets are their own additive component (`stepFun`, or `erfFun` for a
+broadened onset) rather than a `y0` parameter on every function.
 Convolution kernels are named `funcCONV` (e.g. `gaussCONV`) with a
 companion `funcCONV_kernel_width(...)` returning the kernel-width multiplier.
 Add new time-domain behavior or IRF kernels here.
diff --git a/examples/fitting_workflows/01_basic_fitting/data/generate_data.ipynb b/examples/fitting_workflows/01_basic_fitting/data/generate_data.ipynb
index 401e81c..5491338 100644
--- a/examples/fitting_workflows/01_basic_fitting/data/generate_data.ipynb
+++ b/examples/fitting_workflows/01_basic_fitting/data/generate_data.ipynb
@@ -169,7 +169,7 @@
     "| `gaussCONV SD` | 3 | IRF width |\n",
     "| `expFun A` | 5 | Peak-shift amplitude |\n",
     "| `expFun tau` | 50 | Shift decay time constant |\n",
-    "| `expFun t0` / `y0` | 0 / 0 | Time zero / offset |"
+    "| `expFun t0` | 0 | Time zero |"
    ]
   }
  ],
diff --git a/examples/fitting_workflows/01_basic_fitting/data/models_time_truth.yaml b/examples/fitting_workflows/01_basic_fitting/data/models_time_truth.yaml
index d3dcaa1..360be80 100644
--- a/examples/fitting_workflows/01_basic_fitting/data/models_time_truth.yaml
+++ b/examples/fitting_workflows/01_basic_fitting/data/models_time_truth.yaml
@@ -8,4 +8,3 @@ MonoExpIRF:
         A: [5, True]
         tau: [50, True]
         t0: [0, True]
-        y0: [0, True]
diff --git a/examples/fitting_workflows/01_basic_fitting/models_time.yaml b/examples/fitting_workflows/01_basic_fitting/models_time.yaml
index 8fa610d..640f464 100644
--- a/examples/fitting_workflows/01_basic_fitting/models_time.yaml
+++ b/examples/fitting_workflows/01_basic_fitting/models_time.yaml
@@ -8,4 +8,3 @@ MonoExpIRF:
         A: [5, True, 0, 15]
         tau: [10, True, 1, 100]
         t0: [0, False]
-        y0: [0, False]
diff --git a/examples/fitting_workflows/02_dependent_parameters/models_time.yaml b/examples/fitting_workflows/02_dependent_parameters/models_time.yaml
index 9475358..580cfb7 100644
--- a/examples/fitting_workflows/02_dependent_parameters/models_time.yaml
+++ b/examples/fitting_workflows/02_dependent_parameters/models_time.yaml
@@ -7,4 +7,3 @@ sine_wave:
         f: [0.25, True, 1E-3, 10]
         phi: [0, True, -0.75, 0.75]
         t0: [0, False]
-        y0: [0, False]
diff --git a/examples/fitting_workflows/03_multi_cycle_dynamics/data/models_time_truth.yaml b/examples/fitting_workflows/03_multi_cycle_dynamics/data/models_time_truth.yaml
index eca07ad..c6a84cf 100644
--- a/examples/fitting_workflows/03_multi_cycle_dynamics/data/models_time_truth.yaml
+++ b/examples/fitting_workflows/03_multi_cycle_dynamics/data/models_time_truth.yaml
@@ -13,7 +13,6 @@ MonoExpNeg:
         A: [-2, True]
         tau: [0.4, True]
         t0: [0, True]
-        y0: [0, True]
 
 # subcycle 2: amplitude mirrors subcycle 1 via expression, tau independent
 MonoExpPos:
@@ -21,4 +20,3 @@ MonoExpPos:
         A: ["-expFun_01_A"]
         tau: [0.8, True]
         t0: [0, True]
-        y0: [0, True]
diff --git a/examples/fitting_workflows/03_multi_cycle_dynamics/models_time.yaml b/examples/fitting_workflows/03_multi_cycle_dynamics/models_time.yaml
index 0fccaea..5076e97 100644
--- a/examples/fitting_workflows/03_multi_cycle_dynamics/models_time.yaml
+++ b/examples/fitting_workflows/03_multi_cycle_dynamics/models_time.yaml
@@ -12,13 +12,12 @@ IRF:
         SD: [0.1, True, 0, 1]
 
 # Subcycle 1: negative exponential shift, relaxing back to 0.
-# t0/y0 are on the subcycle's local clock, which resets at each subcycle start.
+# t0 is on the subcycle's local clock, which resets at each subcycle start.
 MonoExpNeg:
     expFun:
         A: [-1.5, True, -5, -1E-3]
         tau: [1, True, 0.05, 5]
         t0: [0, False]
-        y0: [0, False]
 
 # Subcycle 2: amplitude mirrors subcycle 1 via an expression (no extra free
 # parameter), while tau is fit independently — expressions work across
@@ -28,4 +27,3 @@ MonoExpPos:
         A: ["-expFun_01_A"]
         tau: [1, True, 0.05, 5]
         t0: [0, False]
-        y0: [0, False]
diff --git a/examples/fitting_workflows/04_parameter_profiles/data/models_time_truth.yaml b/examples/fitting_workflows/04_parameter_profiles/data/models_time_truth.yaml
index 0af788e..6435757 100644
--- a/examples/fitting_workflows/04_parameter_profiles/data/models_time_truth.yaml
+++ b/examples/fitting_workflows/04_parameter_profiles/data/models_time_truth.yaml
@@ -14,4 +14,3 @@ BandBendingRecovery:
         A: [0.5, True]
         tau: [100, True]
         t0: [0, True]
-        y0: [0, True]
diff --git a/examples/fitting_workflows/04_parameter_profiles/models_time.yaml b/examples/fitting_workflows/04_parameter_profiles/models_time.yaml
index 7129d22..70c9b11 100644
--- a/examples/fitting_workflows/04_parameter_profiles/models_time.yaml
+++ b/examples/fitting_workflows/04_parameter_profiles/models_time.yaml
@@ -26,4 +26,3 @@ BandBendingRecovery:
         A:   [0.3, True, 0, 1]      # collapse amplitude (eV/nm)
         tau: [50, True, 5, 500]     # recovery time constant (ps)
         t0:  [0, False]
-        y0:  [0, False]
diff --git a/examples/fitting_workflows/10_model_comparison/models_time.yaml b/examples/fitting_workflows/10_model_comparison/models_time.yaml
index b30d354..48469f4 100644
--- a/examples/fitting_workflows/10_model_comparison/models_time.yaml
+++ b/examples/fitting_workflows/10_model_comparison/models_time.yaml
@@ -10,7 +10,6 @@ MonoExpPosIRF:
         A: [2, True, 0, 10]
         tau: [30, True, 1, 100]
         t0: [0, False]
-        y0: [0, False]
 
 # MonoExpPos — same exponential rise, no IRF convolution. Sharp turn-on
 # at t = t0; will undershoot the smoothed rise around t = 0.
@@ -19,4 +18,3 @@ MonoExpPos:
         A: [2, True, 0, 10]
         tau: [30, True, 1, 100]
         t0: [0, False]
-        y0: [0, False]
diff --git a/examples/fitting_workflows/20_multi_file_independent_fit/models_time.yaml b/examples/fitting_workflows/20_multi_file_independent_fit/models_time.yaml
index ceac58f..1bed166 100644
--- a/examples/fitting_workflows/20_multi_file_independent_fit/models_time.yaml
+++ b/examples/fitting_workflows/20_multi_file_independent_fit/models_time.yaml
@@ -8,4 +8,3 @@ MonoExpPosIRF:
         A: [5, file, 0, 15]
         tau: [10, project, 1, 100]
         t0: [0, False]
-        y0: [0, False]
diff --git a/examples/fitting_workflows/21_multi_file_shared_fit/data/generate_data.ipynb b/examples/fitting_workflows/21_multi_file_shared_fit/data/generate_data.ipynb
index ac8f331..ec2134a 100644
--- a/examples/fitting_workflows/21_multi_file_shared_fit/data/generate_data.ipynb
+++ b/examples/fitting_workflows/21_multi_file_shared_fit/data/generate_data.ipynb
@@ -177,8 +177,7 @@
     "| `gaussCONV SD` | 3 | IRF width |\n",
     "| `expFun A` | varies | Shift amplitude |\n",
     "| `expFun tau` | 50 | Decay time constant |\n",
-    "| `expFun t0` | 0 | Time zero |\n",
-    "| `expFun y0` | 0 | Offset |"
+    "| `expFun t0` | 0 | Time zero |"
    ]
   }
  ],
diff --git a/examples/fitting_workflows/21_multi_file_shared_fit/data/models_time_truth.yaml b/examples/fitting_workflows/21_multi_file_shared_fit/data/models_time_truth.yaml
index 6c30ffe..9a48413 100644
--- a/examples/fitting_workflows/21_multi_file_shared_fit/data/models_time_truth.yaml
+++ b/examples/fitting_workflows/21_multi_file_shared_fit/data/models_time_truth.yaml
@@ -7,4 +7,3 @@ MonoExpPosIRF:
         A: [5, True, 0, 15]
         tau: [50, True, 1, 100]
         t0: [0, False, 0, 1]
-        y0: [0, False, 0, 1]
diff --git a/examples/fitting_workflows/21_multi_file_shared_fit/models_time.yaml b/examples/fitting_workflows/21_multi_file_shared_fit/models_time.yaml
index ceac58f..1bed166 100644
--- a/examples/fitting_workflows/21_multi_file_shared_fit/models_time.yaml
+++ b/examples/fitting_workflows/21_multi_file_shared_fit/models_time.yaml
@@ -8,4 +8,3 @@ MonoExpPosIRF:
         A: [5, file, 0, 15]
         tau: [10, project, 1, 100]
         t0: [0, False]
-        y0: [0, False]
diff --git a/examples/synthetic_data/01_simulator/models_time.yaml b/examples/synthetic_data/01_simulator/models_time.yaml
index 2efeb06..fbad753 100644
--- a/examples/synthetic_data/01_simulator/models_time.yaml
+++ b/examples/synthetic_data/01_simulator/models_time.yaml
@@ -8,20 +8,18 @@ IRF:
         SD: [5.0E-2, True, 0, 1]
              
 MonoExpPos:
-    # expFun(t, A, tau, t0, y0)
+    # expFun(t, A, tau, t0)
     expFun:
         A: [1, True, 0, 5]
         tau: [2.5, True, 1, 10]
         t0: [0, False, 0, 1]
-        y0: [0, False, 0, 1]
 
 MonoExpNeg:
-    # expFun(t, A, tau, t0, y0)
+    # expFun(t, A, tau, t0)
     expFun:
         A: [-1, True, -5, 0]
         tau: [2.5, True, 1, 10]
         t0: [0, False, 0, 1]
-        y0: [0, False, 0, 1]
 
 MonoExpPosIRF:
     gaussCONV:
@@ -30,4 +28,3 @@ MonoExpPosIRF:
         A: [5, True, 0, 15]
         tau: [50, True, 1, 100]
         t0: [0, False, 0, 1]
-        y0: [0, False, 0, 1]
diff --git a/examples/synthetic_data/02_ml_training_data/models_time.yaml b/examples/synthetic_data/02_ml_training_data/models_time.yaml
index 4fd2839..9af37f8 100644
--- a/examples/synthetic_data/02_ml_training_data/models_time.yaml
+++ b/examples/synthetic_data/02_ml_training_data/models_time.yaml
@@ -8,20 +8,18 @@ IRF:
         SD: [5.0E-2, True, 0, 1]
              
 MonoExpPos:
-    # expFun(t, A, tau, t0, y0)
+    # expFun(t, A, tau, t0)
     expFun:
         A: [1, True, 0, 5]
         tau: [2.5, True, 1, 10]
         t0: [0, False, 0, 1]
-        y0: [0, False, 0, 1]
 
 MonoExpNeg:
-    # expFun(t, A, tau, t0, y0)
+    # expFun(t, A, tau, t0)
     expFun:
         A: [-1, True, -5, 0]
         tau: [2.5, True, 1, 10]
         t0: [0, False, 0, 1]
-        y0: [0, False, 0, 1]
 
 MonoExpPosIRF:
     gaussCONV:
@@ -30,4 +28,3 @@ MonoExpPosIRF:
         A: [5, True, 0, 15]
         tau: [50, True, 1, 100]
         t0: [0, False, 0, 1]
-        y0: [0, False, 0, 1]
diff --git a/pyproject.toml b/pyproject.toml
index c27ad97..b06690d 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,9 +4,9 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "trspecfit"
-version = "0.9.3"
+version = "0.10.0"
 authors = [
-    {name = "Johannes Mahl", email = "johannes.a.mahl@gmail.com"},
+    {name = "Johannes Mahl", email = "infinitymonkeyatwork@gmail.com"},
 ]
 description = "Fit 2D time- and energy-resolved spectroscopy data"
 keywords = [
diff --git a/src/trspecfit/eval_2d.py b/src/trspecfit/eval_2d.py
index b92f496..88310bf 100644
--- a/src/trspecfit/eval_2d.py
+++ b/src/trspecfit/eval_2d.py
@@ -29,12 +29,13 @@
 # ---------------------------------------------------------------------------
 
 DYNAMICS_DISPATCH: dict[int, tuple] = {
-    DynFuncKind.EXPFUN: (fcts_time.expFun, 4),
-    DynFuncKind.SINFUN: (fcts_time.sinFun, 5),
-    DynFuncKind.LINFUN: (fcts_time.linFun, 3),
-    DynFuncKind.SINDIVX: (fcts_time.sinDivX, 4),
-    DynFuncKind.ERFFUN: (fcts_time.erfFun, 4),
-    DynFuncKind.SQRTFUN: (fcts_time.sqrtFun, 3),
+    DynFuncKind.EXPFUN: (fcts_time.expFun, 3),
+    DynFuncKind.SINFUN: (fcts_time.sinFun, 4),
+    DynFuncKind.LINFUN: (fcts_time.linFun, 2),
+    DynFuncKind.SINDIVX: (fcts_time.sinDivX, 3),
+    DynFuncKind.ERFFUN: (fcts_time.erfFun, 3),
+    DynFuncKind.SQRTFUN: (fcts_time.sqrtFun, 2),
+    DynFuncKind.STEPFUN: (fcts_time.stepFun, 2),
 }
 
 # Convolution kernel dispatch: kernel function evaluated on the frozen
diff --git a/src/trspecfit/functions/energy.py b/src/trspecfit/functions/energy.py
index 94f7554..e8646c1 100644
--- a/src/trspecfit/functions/energy.py
+++ b/src/trspecfit/functions/energy.py
@@ -255,7 +255,7 @@ def Voigt(x: np.ndarray, A: float, x0: float, SD: float, W: float) -> np.ndarray
     x : ndarray
         Energy axis
     A : float
-        Peak amplitude (maximum value, approximately, for narrow peaks)
+        Peak amplitude (maximum value, at x = x0)
     x0 : float
         Peak center position
     SD : float
@@ -270,8 +270,9 @@ def Voigt(x: np.ndarray, A: float, x0: float, SD: float, W: float) -> np.ndarray
     """
 
     voigt = np.real(wofz(((x - x0) + 1j * (W / 2)) / SD / np.sqrt(2)))
-    max_voigt = np.max(voigt, axis=-1, keepdims=True)
-    return np.asarray(A * voigt / max_voigt)
+    # analytic peak value (profile at x = x0), keeps A grid-independent
+    peak_voigt = np.real(wofz(1j * (W / 2) / SD / np.sqrt(2)))
+    return np.asarray(A * voigt / peak_voigt)
 
 
 #
@@ -295,6 +296,8 @@ def GLS(x: np.ndarray, A: float, x0: float, F: float, m: float) -> np.ndarray:
         - m = 1: Pure Lorentzian
         - 0 < m < 1: Weighted mixture
         Typical value: m ≈ 0.3
+        Keep m in [0, 1] (set parameter bounds accordingly); values
+        outside give negative mixture weights and unphysical shapes.
 
     Returns
     -------
@@ -327,6 +330,8 @@ def GLP(x: np.ndarray, A: float, x0: float, F: float, m: float) -> np.ndarray:
         - m = 1: Pure Lorentzian
         - 0 < m < 1: Hybrid shape
         Typical value: m ≈ 0.3
+        Keep m in [0, 1] (set parameter bounds accordingly); m < 0 makes
+        the denominator 1 + 4*m*u² cross zero, producing NaN/inf.
 
     Returns
     -------
@@ -359,6 +364,8 @@ def DS(x: np.ndarray, A: float, x0: float, F: float, alpha: float) -> np.ndarray
         - 0 < alpha < 0.3: Typical for metals (e.g., Al: 0.10-0.15)
         - Larger alpha: Stronger asymmetry, more pronounced tail
         - Range: typically 0-0.5 for physical systems
+        The asymmetric tail extends toward x < x0, correct for a
+        kinetic-energy axis; on a binding-energy axis it appears mirrored.
 
     Returns
     -------
diff --git a/src/trspecfit/functions/time.py b/src/trspecfit/functions/time.py
index 4b36a22..001bfc8 100644
--- a/src/trspecfit/functions/time.py
+++ b/src/trspecfit/functions/time.py
@@ -6,11 +6,10 @@
 Use CamelCase naming (UpperCamelCase or lowerCamelCase) for function names.
 
 **Dynamics Functions:**
-Signature: func(t, par1, par2, ..., t0, y0)
+Signature: func(t, par1, par2, ..., t0)
 - t: Time axis (numpy array)
 - par1, par2, ...: Function-specific parameters
 - t0: Time zero (function starts at this time)
-- y0: Offset value (baseline)
 - Returns: f(t) = 0 for t < t0, dynamics for t >= t0
 
 **Convolution Kernels:**
@@ -23,12 +22,15 @@
 **Time Zero Convention:**
 All dynamics functions are zero before t0 and activate at t >= t0.
 This reflects physical causality: response begins after excitation.
+Exception: erfFun is a smoothed step, so it is nonzero (approaching 0)
+shortly before t0; it crosses A/2 at t0.
 
-**Offset Convention:**
-Parameter y0 sets the asymptotic value or baseline.
-- Decays: approach y0 as t → ∞
-- Rises: start from 0, reach y0 + A
-- Oscillations: oscillate around y0
+**Offsets and Baselines:**
+Dynamics components combine by addition, so constant offsets are their
+own component rather than a parameter of every function:
+- stepFun: sharp onset to a constant value A at t0
+- erfFun: Gaussian-broadened onset, rises from 0 to A around t0
+For example, a decay to a nonzero plateau is expFun + stepFun sharing t0.
 
 **Time Resolution:**
 Functions inherit time axis from Dynamics model. Consider:
@@ -42,7 +44,6 @@
 - A: Amplitude (change in signal)
 - tau: Time constant (decay/rise time, 1/e point)
 - t0: Time zero (start of dynamics)
-- y0: Offset/baseline value
 - f: Frequency (for oscillations)
 - phi: Phase (for oscillations)
 - SD: Standard deviation (for Gaussian kernels)
@@ -93,7 +94,35 @@ def none(t: np.ndarray) -> np.ndarray:
 
 
 #
-def linFun(t: np.ndarray, m: float, t0: float, y0: float) -> np.ndarray:
+def stepFun(t: np.ndarray, A: float, t0: float) -> np.ndarray:
+    """
+    Step function (constant offset switching on at t0).
+
+    The causal offset primitive: dynamics components combine by addition,
+    so add stepFun to model baselines or plateaus (e.g. expFun + stepFun
+    sharing t0 gives a decay to a nonzero plateau). For a Gaussian-broadened
+    onset use erfFun instead.
+
+    Parameters
+    ----------
+    t : ndarray
+        Time axis
+    A : float
+        Step height (constant value for t >= t0)
+    t0 : float
+        Time zero (onset of the step)
+
+    Returns
+    -------
+    ndarray
+        Step function: 0 for t<t0, A for t>=t0
+    """
+
+    return np.where(t < t0, 0.0, A)
+
+
+#
+def linFun(t: np.ndarray, m: float, t0: float) -> np.ndarray:
     """
     Linear dynamics (constant rate of change).
 
@@ -107,20 +136,18 @@ def linFun(t: np.ndarray, m: float, t0: float, y0: float) -> np.ndarray:
         - m < 0: Linear decrease
     t0 : float
         Time zero (start of linear change)
-    y0 : float
-        Offset value at t0 (initial value)
 
     Returns
     -------
     ndarray
-        Linear function: 0 for t<t0, m*(t-t0)+y0 for t>=t0
+        Linear function: 0 for t<t0, m*(t-t0) for t>=t0
     """
 
-    return np.where(t < t0, 0.0, m * (t - t0) + y0)
+    return np.where(t < t0, 0.0, m * (t - t0))
 
 
 #
-def expFun(t: np.ndarray, A: float, tau: float, t0: float, y0: float) -> np.ndarray:
+def expFun(t: np.ndarray, A: float, tau: float, t0: float) -> np.ndarray:
     """
     Exponential decay or rise dynamics.
 
@@ -130,29 +157,25 @@ def expFun(t: np.ndarray, A: float, tau: float, t0: float, y0: float) -> np.ndar
         Time axis
     A : float
         Amplitude (initial change at t0).
-        - A > 0: Decay from y0+A to y0
-        - A < 0: Rise from y0 to y0+|A|
+        - A > 0: Jumps to A at t0, decays toward 0
+        - A < 0: Jumps to -|A| at t0, rises toward 0
     tau : float
         Time constant (1/e time). Units: [time units]
         At t = t0 + tau, signal changes by factor of e (≈2.718)
     t0 : float
         Time zero (start of exponential)
-    y0 : float
-        Asymptotic value (baseline as t → ∞)
 
     Returns
     -------
     ndarray
-        Exponential: 0 for t<t0, A*exp(-(t-t0)/tau)+y0 for t>=t0
+        Exponential: 0 for t<t0, A*exp(-(t-t0)/tau) for t>=t0
     """
 
-    return np.where(t < t0, 0.0, A * np.exp(-1 / tau * (t - t0)) + y0)
+    return np.where(t < t0, 0.0, A * np.exp(-1 / tau * (t - t0)))
 
 
 #
-def sinFun(
-    t: np.ndarray, A: float, f: float, phi: float, t0: float, y0: float
-) -> np.ndarray:
+def sinFun(t: np.ndarray, A: float, f: float, phi: float, t0: float) -> np.ndarray:
     """
     Sinusoidal oscillations (coherent dynamics).
 
@@ -172,20 +195,19 @@ def sinFun(
         - phi = π: Starts at zero (negative slope)
     t0 : float
         Time zero (start of oscillation)
-    y0 : float
-        Offset (center line of oscillation)
 
     Returns
     -------
     ndarray
-        Sinusoid: 0 for t<t0, A*sin(2πf(t-t0)+phi)+y0 for t>=t0
+        Sinusoid: 0 for t<t0, A*sin(2πf(t-t0)+phi) for t>=t0
+        Oscillates around 0; add stepFun to shift the center line.
     """
 
-    return np.where(t < t0, 0.0, A * np.sin(2 * np.pi * f * (t - t0) + phi) + y0)
+    return np.where(t < t0, 0.0, A * np.sin(2 * np.pi * f * (t - t0) + phi))
 
 
 #
-def sinDivX(t: np.ndarray, A: float, f: float, t0: float, y0: float) -> np.ndarray:
+def sinDivX(t: np.ndarray, A: float, f: float, t0: float) -> np.ndarray:
     """
     Damped sinc function: sin(x)/x oscillation.
 
@@ -199,50 +221,49 @@ def sinDivX(t: np.ndarray, A: float, f: float, t0: float, y0: float) -> np.ndarr
         Frequency in [1/time units]
     t0 : float
         Time zero (start of oscillation)
-    y0 : float
-        Offset value
 
     Returns
     -------
     ndarray
-        Sinc oscillation: 0 for t<t0, A*sin(2πf(t-t0))/(2πf(t-t0))+y0 for t>=t0
+        Sinc oscillation: 0 for t<t0, A*sin(2πf(t-t0))/(2πf(t-t0)) for t>=t0
     """
 
     # np.sinc(u) = sin(pi*u)/(pi*u), so u=2*f*(t-t0) gives sin(2*pi*f*dt)/(2*pi*f*dt)
-    return np.where(t < t0, 0.0, A * np.sinc(2 * f * (t - t0)) + y0)
+    return np.where(t < t0, 0.0, A * np.sinc(2 * f * (t - t0)))
 
 
 #
-def erfFun(t: np.ndarray, A: float, SD: float, t0: float, y0: float) -> np.ndarray:
+def erfFun(t: np.ndarray, A: float, SD: float, t0: float) -> np.ndarray:
     """
     Error function rise (step with Gaussian broadening).
-    erfFun ≈ step ⊗ Gaussian(SD)
+    erfFun ≈ stepFun ⊗ Gaussian(SD)
+
+    As a smoothed step this is the one dynamics function that is nonzero
+    (approaching 0) shortly before t0; it crosses A/2 at t0 and rises to A.
 
     Parameters
     ----------
     t : ndarray
         Time axis
     A : float
-        Amplitude (total change from initial to final value)
+        Amplitude (final value, asymptote as t → ∞)
     SD : float
         Standard deviation of Gaussian broadening (rise time ~2.355*SD)
         Smaller SD → sharper rise
     t0 : float
         Center of rise (50% point)
-    y0 : float
-        Final value (asymptote as t → ∞)
 
     Returns
     -------
     ndarray
-        Error function: A/2 * (1 + erf((t-t0)/(SD*√2))) + y0
+        Error function: A/2 * (1 + erf((t-t0)/(SD*√2)))
     """
 
-    return np.asarray(A / 2 * (1 + erf((t - t0) / (SD * np.sqrt(2)))) + y0)
+    return np.asarray(A / 2 * (1 + erf((t - t0) / (SD * np.sqrt(2)))))
 
 
 #
-def sqrtFun(t: np.ndarray, A: float, t0: float, y0: float) -> np.ndarray:
+def sqrtFun(t: np.ndarray, A: float, t0: float) -> np.ndarray:
     """
     Square root rise (diffusion dynamics).
 
@@ -254,17 +275,15 @@ def sqrtFun(t: np.ndarray, A: float, t0: float, y0: float) -> np.ndarray:
         Amplitude scaling factor
     t0 : float
         Time zero (start of diffusion)
-    y0 : float
-        Offset value
 
     Returns
     -------
     ndarray
-        Square root rise: 0 for t<t0, A*√(t-t0)+y0 for t>=t0
+        Square root rise: 0 for t<t0, A*√(t-t0) for t>=t0
     """
 
     # numpy array .clip sets all t<t0 to zero
-    return np.asarray(A * np.sqrt((t - t0).clip(0)) + y0)
+    return np.asarray(A * np.sqrt((t - t0).clip(0)))
 
 
 #
@@ -434,7 +453,11 @@ def expDecayCONV_kernel_width(tau: float | None = None) -> int:
 #
 def expRiseCONV(x: np.ndarray, tau: float) -> np.ndarray:
     """
-    Causal exponential rise kernel.
+    Anti-causal exponential rise kernel (mirror of expDecayCONV).
+
+    The kernel is nonzero only for x <= 0, so the convolved response at
+    time t draws from the signal at later times: it rises before the
+    excitation and saturates at t0.
 
     Parameters
     ----------
@@ -474,7 +497,7 @@ def boxCONV(x: np.ndarray, width: float) -> np.ndarray:
     Returns
     -------
     ndarray
-        Rectangular function: 1 inside width, 0 outside (with smooth edges)
+        Rectangular function: 1 inside width, 0 outside (hard edges)
     """
 
     return np.where(np.abs(x) <= width / 2, 1.0, 0.0)
diff --git a/src/trspecfit/graph_ir.py b/src/trspecfit/graph_ir.py
index 804f590..64f5866 100644
--- a/src/trspecfit/graph_ir.py
+++ b/src/trspecfit/graph_ir.py
@@ -128,6 +128,7 @@ class DynFuncKind(IntEnum):
     SINDIVX = 3
     ERFFUN = 4
     SQRTFUN = 5
+    STEPFUN = 6
 
 
 #
@@ -157,6 +158,7 @@ class ParamSourceKind(IntEnum):
     "sinDivX": DynFuncKind.SINDIVX,
     "erfFun": DynFuncKind.ERFFUN,
     "sqrtFun": DynFuncKind.SQRTFUN,
+    "stepFun": DynFuncKind.STEPFUN,
 }
 
 _FUNCTION_NAME_TO_PROFILE_FUNC: dict[str, ProfileFuncKind] = {
@@ -2894,6 +2896,7 @@ def schedule_2d(graph: GraphIR) -> ScheduledPlan2D:
         int(DynFuncKind.SINDIVX): fcts_time.sinDivX,
         int(DynFuncKind.ERFFUN): fcts_time.erfFun,
         int(DynFuncKind.SQRTFUN): fcts_time.sqrtFun,
+        int(DynFuncKind.STEPFUN): fcts_time.stepFun,
     }
 
     _CONV_KERNEL_DISPATCH: dict[int, Callable[..., Any]] = {
diff --git a/tests/models/file_time.yaml b/tests/models/file_time.yaml
index 74e1cdc..034defe 100644
--- a/tests/models/file_time.yaml
+++ b/tests/models/file_time.yaml
@@ -19,12 +19,11 @@ IRF:
         SD: [5.0E-2, True, 0, 1]
              
 MonoExpPos:
-    # expFun(t, A, tau, t0, y0)
+    # expFun(t, A, tau, t0)
     expFun:
         A: [1, True, 0, 5]
         tau: [2.5, True, 1, 10]
         t0: [0, False, 0, 1]
-        y0: [0, False, 0, 1]
 
 # stronger amplitude for profile round-trip tests (~20% modulation)
 MonoExpPosStrong:
@@ -32,15 +31,55 @@ MonoExpPosStrong:
         A: [2, True, 0, 10]
         tau: [2.5, True, 1, 10]
         t0: [0, False, 0, 1]
-        y0: [0, False, 0, 1]
 
 MonoExpNeg:
-    # expFun(t, A, tau, t0, y0)
+    # expFun(t, A, tau, t0)
     expFun:
         A: [-1, True, -5, 0]
         tau: [2.5, True, 1, 10]
         t0: [0, False, 0, 1]
-        y0: [0, False, 0, 1]
+
+MonoStep:
+    # stepFun(t, A, t0)
+    stepFun:
+        A: [1.5, True, 0, 5]
+        t0: [0, False, 0, 1]
+
+MonoSin:
+    # sinFun(t, A, f, phi, t0)
+    # phi initialized nonzero (with t0=0) so a phi/t0 ordering bug in the
+    # lowered path diverges from the interpreter even at initial theta.
+    sinFun:
+        A: [1, True, 0, 5]
+        f: [0.05, True, 0.01, 0.2]
+        phi: [0.3, True, -3.15, 3.15]
+        t0: [0, False, 0, 1]
+
+MonoLin:
+    # linFun(t, m, t0)
+    linFun:
+        m: [0.02, True, -1, 1]
+        t0: [0, False, 0, 1]
+
+MonoSinDivX:
+    # sinDivX(t, A, f, t0)
+    sinDivX:
+        A: [1, True, 0, 5]
+        f: [0.05, True, 0.01, 0.2]
+        t0: [0, False, 0, 1]
+
+MonoErf:
+    # erfFun(t, A, SD, t0)
+    erfFun:
+        A: [1, True, 0, 5]
+        SD: [5, True, 1, 20]
+        t0: [20, False, 0, 50]
+
+MonoSqrt:
+    # sqrtFun(t, A, t0)
+    sqrtFun:
+        A: [0.2, True, 0, 5]
+        t0: [0, False, 0, 1]
 
 MonoExpPosIRF:
     gaussCONV:
@@ -49,7 +88,6 @@ MonoExpPosIRF:
         A: [1, True, 0, 5]
         tau: [2.5, True, 1, 10]
         t0: [0, False, 0, 1]
-        y0: [0, False, 0, 1]
 
 MonoExpPosLorentzIRF:
     lorentzCONV:
@@ -58,7 +96,6 @@ MonoExpPosLorentzIRF:
         A: [1, True, 0, 5]
         tau: [2.5, True, 1, 10]
         t0: [0, False, 0, 1]
-        y0: [0, False, 0, 1]
 
 MonoExpPosVoigtIRF:
     voigtCONV:
@@ -68,7 +105,6 @@ MonoExpPosVoigtIRF:
         A: [1, True, 0, 5]
         tau: [2.5, True, 1, 10]
         t0: [0, False, 0, 1]
-        y0: [0, False, 0, 1]
 
 MonoExpPosExpSymIRF:
     expSymCONV:
@@ -77,7 +113,6 @@ MonoExpPosExpSymIRF:
         A: [1, True, 0, 5]
         tau: [2.5, True, 1, 10]
         t0: [0, False, 0, 1]
-        y0: [0, False, 0, 1]
 
 MonoExpPosExpDecayIRF:
     expDecayCONV:
@@ -86,7 +121,6 @@ MonoExpPosExpDecayIRF:
         A: [1, True, 0, 5]
         tau: [2.5, True, 1, 10]
         t0: [0, False, 0, 1]
-        y0: [0, False, 0, 1]
 
 MonoExpPosExpRiseIRF:
     expRiseCONV:
@@ -95,7 +129,6 @@ MonoExpPosExpRiseIRF:
         A: [1, True, 0, 5]
         tau: [2.5, True, 1, 10]
         t0: [0, False, 0, 1]
-        y0: [0, False, 0, 1]
 
 MonoExpPosBoxIRF:
     boxCONV:
@@ -104,7 +137,6 @@ MonoExpPosBoxIRF:
         A: [1, True, 0, 5]
         tau: [2.5, True, 1, 10]
         t0: [0, False, 0, 1]
-        y0: [0, False, 0, 1]
 
 # Empty subcycle placeholder (applies to all times in multi-cycle)
 ModelNone:
@@ -118,7 +150,6 @@ MonoExpPosExpr:
         A: ["-expFun_01_A"]
         tau: ["expFun_01_tau"]
         t0: [0, False, 0, 1]
-        y0: [0, False, 0, 1]
 
 # Bi-exponential: two expFun, second shares t0 with first via expression.
 # Single-cycle, no convolution -- fully lowerable.
@@ -127,12 +158,10 @@ BiExpSharedT0:
         A: [3, True, 0, 10]
         tau: [5, True, 1, 20]
         t0: [0, False, 0, 1]
-        y0: [0, False, 0, 1]
     expFun:
         A: [-1, True, -5, 0]
         tau: [20, True, 5, 50]
         t0: ["expFun_01_t0"]
-        y0: [0, False, 0, 1]
 
 # convolution as last component (should fail ordering validation)
 conv_last:
@@ -140,6 +169,5 @@ conv_last:
         A: [1, True, 0, 5]
         tau: [2.5, True, 1, 10]
         t0: [0, False, 0, 1]
-        y0: [0, False, 0, 1]
     gaussCONV:
         SD: [5.0E-2, True, 0, 1]
\ No newline at end of file
diff --git a/tests/models/project_time.yaml b/tests/models/project_time.yaml
index e0326b8..084d814 100644
--- a/tests/models/project_time.yaml
+++ b/tests/models/project_time.yaml
@@ -6,7 +6,6 @@ MonoExpProject:
         A: [5, "file", 0, 15]
         tau: [5, "project", 1, 100]
         t0: [0, "static", 0, 1]
-        y0: [0, "static", 0, 1]
 
 # Bi-exponential decay with Gaussian IRF — t0 of second component is
 # an expression referencing t0 of first (project-vary). Tests that
@@ -20,9 +19,7 @@ BiExpProject:
         A: [2, "file", -20, 20]
         tau: [2, "project", 0.1, 50]
         t0: [3, "project", -5, 10]
-        y0: [0, "static", 0, 1]
     expFun:
         A: [1, "file", -20, 20]
         tau: [20, "project", 1, 200]
         t0: ["expFun_01_t0"]
-        y0: [0, "static", 0, 1]
diff --git a/tests/test_config_functions.py b/tests/test_config_functions.py
index 2117fef..864b8c6 100644
--- a/tests/test_config_functions.py
+++ b/tests/test_config_functions.py
@@ -51,10 +51,13 @@ def test_LinBack_strips_spectrum(self):
     # time: dynamics functions — first arg (t) stripped
     #
     def test_expFun(self):
-        assert get_function_parameters("expFun") == ["A", "tau", "t0", "y0"]
+        assert get_function_parameters("expFun") == ["A", "tau", "t0"]
 
     def test_sinFun(self):
-        assert get_function_parameters("sinFun") == ["A", "f", "phi", "t0", "y0"]
+        assert get_function_parameters("sinFun") == ["A", "f", "phi", "t0"]
+
+    def test_stepFun(self):
+        assert get_function_parameters("stepFun") == ["A", "t0"]
 
     def test_none_has_no_params(self):
         """none() takes only t, so parameters should be empty."""
diff --git a/tests/test_evaluate_2d.py b/tests/test_evaluate_2d.py
index 6e82737..6aefea4 100644
--- a/tests/test_evaluate_2d.py
+++ b/tests/test_evaluate_2d.py
@@ -249,6 +249,49 @@ def test_glp_time_dep_amplitude(self):
         F_idx = plan.opt_param_names.index("GLP_01_F")
         _perturb_theta(plan, model, theta, [A_idx, F_idx], [5.0, 0.3])
 
+    # Per-function parity for every non-expFun dynamics driver: each must
+    # match MCP at initial theta AND under perturbation of its params.
+    # expFun is already exercised pervasively via MonoExpPos; this list
+    # covers the remaining DynFuncKind enum/dispatch entries and their
+    # param counts (all decremented when y0 was removed). Perturbing more
+    # than one param where available also catches param-order bugs.
+    _DYN_DRIVER_PARITY_CASES = [
+        ("MonoStep", [("stepFun_01_A", 1.0)]),
+        (
+            "MonoSin",
+            [("sinFun_01_A", 0.5), ("sinFun_01_f", 0.02), ("sinFun_01_phi", 0.1)],
+        ),
+        ("MonoLin", [("linFun_01_m", 0.01)]),
+        ("MonoSinDivX", [("sinDivX_01_A", 0.5), ("sinDivX_01_f", 0.02)]),
+        ("MonoErf", [("erfFun_01_A", 0.5), ("erfFun_01_SD", 1.0)]),
+        ("MonoSqrt", [("sqrtFun_01_A", 0.1)]),
+    ]
+
+    #
+    @pytest.mark.parametrize(
+        ("dyn_model", "param_perturbations"),
+        _DYN_DRIVER_PARITY_CASES,
+        ids=[c[0] for c in _DYN_DRIVER_PARITY_CASES],
+    )
+    def test_dynamics_driver_parity(self, dyn_model, param_perturbations):
+        """evaluate_2d == interpreter for every lowerable dynamics driver."""
+
+        _file, model = _make_2d_model(
+            ["glp_only"],
+            [("GLP_01_A", [dyn_model])],
+        )
+        graph = build_graph(model)
+        assert can_lower_2d(graph)
+        plan = schedule_2d(graph)
+
+        theta = _compare_evaluator_vs_interpreter(model, plan)
+        perturb_idx = [
+            plan.opt_param_names.index(f"GLP_01_A_{suffix}")
+            for suffix, _ in param_perturbations
+        ]
+        perturb_deltas = [delta for _, delta in param_perturbations]
+        _perturb_theta(plan, model, theta, perturb_idx, perturb_deltas)
+
     def test_gls_time_dep_mixing(self):
         """GLS with mixing parameter m time-dependent."""
 
diff --git a/tests/test_functions_energy.py b/tests/test_functions_energy.py
index 2e03a15..9147871 100644
--- a/tests/test_functions_energy.py
+++ b/tests/test_functions_energy.py
@@ -368,6 +368,23 @@ def test_2d_broadcast_normalizes_per_slice(self):
         assert result.shape == (2, x.shape[-1])
         np.testing.assert_allclose(np.max(result, axis=-1), A[:, 0], rtol=1e-3)
 
+    #
+    def test_grid_independent_amplitude(self):
+        """Value at a fixed x must not depend on the sampled window/grid.
+
+        The analytic peak normalization (wofz at dx=0) replaced the old
+        max-over-grid normalization, which rescaled the whole profile when
+        x0 fell outside the energy window.
+        """
+
+        full_x = np.linspace(-10, 10, 2001)  # includes the peak at x0=0
+        idx_probe = np.argmin(np.abs(full_x - 2.0))
+        x_probe = full_x[idx_probe]
+        tail_x = np.linspace(x_probe, 10, 18)  # excludes the peak entirely
+        full = Voigt(full_x, A=4.0, x0=0.0, SD=1.0, W=1.0)
+        tail = Voigt(tail_x, A=4.0, x0=0.0, SD=1.0, W=1.0)
+        assert tail[0] == pytest.approx(full[idx_probe], rel=1e-9)
+
 
 #
 #
diff --git a/tests/test_functions_time.py b/tests/test_functions_time.py
index 07f2003..fe203e4 100644
--- a/tests/test_functions_time.py
+++ b/tests/test_functions_time.py
@@ -17,6 +17,7 @@
     sinDivX,
     sinFun,
     sqrtFun,
+    stepFun,
 )
 
 
@@ -38,47 +39,88 @@ def test_returns_zeros(self):
         assert result.shape == t.shape
 
 
+#
+#
+class TestStepFun:
+    #
+    def test_zero_before_t0(self):
+        t = make_time_axis()
+        result = stepFun(t, A=2.0, t0=0.0)
+        np.testing.assert_allclose(result[t < 0], 0.0)
+
+    #
+    def test_constant_after_t0(self):
+        t = make_time_axis()
+        A = 2.5
+        result = stepFun(t, A=A, t0=0.0)
+        np.testing.assert_allclose(result[t >= 0], A)
+
+    #
+    def test_shifted_t0(self):
+        t = make_time_axis()
+        result = stepFun(t, A=1.0, t0=5.0)
+        np.testing.assert_allclose(result[t < 5.0], 0.0)
+        np.testing.assert_allclose(result[t >= 5.0], 1.0)
+
+    #
+    def test_zero_amplitude(self):
+        """A=0 gives all zeros."""
+
+        t = make_time_axis()
+        result = stepFun(t, A=0.0, t0=0.0)
+        np.testing.assert_allclose(result, 0.0, atol=1e-15)
+
+    #
+    def test_offset_composition(self):
+        """expFun + stepFun sharing t0 reproduces a decay to a plateau."""
+
+        t = make_time_axis()
+        combined = expFun(t, A=5.0, tau=2.0, t0=0.0) + stepFun(t, A=1.0, t0=0.0)
+        np.testing.assert_allclose(combined[t < 0], 0.0)
+        assert combined[-1] == pytest.approx(1.0, abs=1e-6)
+
+
 #
 #
 class TestLinFun:
     #
     def test_zero_before_t0(self):
         t = make_time_axis()
-        result = linFun(t, m=2.0, t0=0.0, y0=1.0)
+        result = linFun(t, m=2.0, t0=0.0)
         np.testing.assert_allclose(result[t < 0], 0.0)
 
     #
     def test_value_at_t0(self):
         t = make_time_axis()
-        result = linFun(t, m=2.0, t0=0.0, y0=1.0)
+        result = linFun(t, m=2.0, t0=0.0)
         idx = np.argmin(np.abs(t - 0.0))
-        assert result[idx] == pytest.approx(1.0, abs=0.01)
+        assert result[idx] == pytest.approx(0.0, abs=0.01)
 
     #
     def test_slope(self):
         t = make_time_axis()
         m = 3.0
-        result = linFun(t, m=m, t0=0.0, y0=0.0)
+        result = linFun(t, m=m, t0=0.0)
         # At t=10, value should be m*10 = 30
         idx = np.argmin(np.abs(t - 10.0))
         assert result[idx] == pytest.approx(30.0, abs=0.1)
 
     #
-    def test_offset_t0(self):
+    def test_shifted_t0(self):
         t = make_time_axis()
-        result = linFun(t, m=1.0, t0=5.0, y0=2.0)
+        result = linFun(t, m=1.0, t0=5.0)
         # Zero before t0=5
         np.testing.assert_allclose(result[t < 5.0], 0.0)
-        # At t=10: m*(10-5) + y0 = 7
+        # At t=10: m*(10-5) = 5
         idx = np.argmin(np.abs(t - 10.0))
-        assert result[idx] == pytest.approx(7.0, abs=0.1)
+        assert result[idx] == pytest.approx(5.0, abs=0.1)
 
     #
     def test_zero_slope(self):
-        """m=0 with y0=0 gives all zeros."""
+        """m=0 gives all zeros."""
 
         t = make_time_axis()
-        result = linFun(t, m=0.0, t0=0.0, y0=0.0)
+        result = linFun(t, m=0.0, t0=0.0)
         np.testing.assert_allclose(result, 0.0, atol=1e-15)
 
     #
@@ -86,7 +128,7 @@ def test_monotonic_positive_slope(self):
         """Positive m: monotonically increasing after t0."""
 
         t = make_time_axis()
-        result = linFun(t, m=2.0, t0=0.0, y0=0.0)
+        result = linFun(t, m=2.0, t0=0.0)
         active = result[t >= 0]
         assert np.all(np.diff(active) >= -1e-12)
 
@@ -97,50 +139,50 @@ class TestExpFun:
     #
     def test_zero_before_t0(self):
         t = make_time_axis()
-        result = expFun(t, A=1.0, tau=5.0, t0=0.0, y0=0.0)
+        result = expFun(t, A=1.0, tau=5.0, t0=0.0)
         np.testing.assert_allclose(result[t < 0], 0.0)
 
     #
     def test_value_at_t0(self):
         t = make_time_axis()
-        result = expFun(t, A=3.0, tau=5.0, t0=0.0, y0=1.0)
+        result = expFun(t, A=3.0, tau=5.0, t0=0.0)
         idx = np.argmin(np.abs(t - 0.0))
-        assert result[idx] == pytest.approx(4.0, abs=0.01)  # A + y0
+        assert result[idx] == pytest.approx(3.0, abs=0.01)  # A
 
     #
-    def test_decay_to_y0(self):
-        """At t >> tau, value approaches y0."""
+    def test_decay_to_zero(self):
+        """At t >> tau, value approaches 0."""
 
         t = make_time_axis()
-        result = expFun(t, A=5.0, tau=2.0, t0=0.0, y0=1.0)
-        assert result[-1] == pytest.approx(1.0, abs=1e-6)
+        result = expFun(t, A=5.0, tau=2.0, t0=0.0)
+        assert result[-1] == pytest.approx(0.0, abs=1e-6)
 
     #
     def test_value_at_one_tau(self):
-        """At t = t0 + tau, value = A*exp(-1) + y0."""
+        """At t = t0 + tau, value = A*exp(-1)."""
 
         t = make_time_axis()
         tau = 5.0
-        result = expFun(t, A=1.0, tau=tau, t0=0.0, y0=0.0)
+        result = expFun(t, A=1.0, tau=tau, t0=0.0)
         idx = np.argmin(np.abs(t - tau))
         assert result[idx] == pytest.approx(np.exp(-1), abs=1e-3)
 
     #
     def test_negative_amplitude_rise(self):
-        """A < 0 gives a rise from y0+A toward y0."""
+        """A < 0 jumps to -|A| at t0 and rises toward 0."""
 
         t = make_time_axis()
-        result = expFun(t, A=-2.0, tau=5.0, t0=0.0, y0=0.0)
+        result = expFun(t, A=-2.0, tau=5.0, t0=0.0)
         idx_t0 = np.argmin(np.abs(t - 0.0))
         assert result[idx_t0] == pytest.approx(-2.0, abs=0.01)
         assert result[-1] == pytest.approx(0.0, abs=1e-3)
 
     #
     def test_zero_amplitude(self):
-        """A=0 with y0=0 gives all zeros."""
+        """A=0 gives all zeros."""
 
         t = make_time_axis()
-        result = expFun(t, A=0.0, tau=5.0, t0=0.0, y0=0.0)
+        result = expFun(t, A=0.0, tau=5.0, t0=0.0)
         np.testing.assert_allclose(result, 0.0, atol=1e-15)
 
     #
@@ -148,7 +190,7 @@ def test_monotonic_decay(self):
         """Positive A: monotonically decreasing after t0."""
 
         t = make_time_axis()
-        result = expFun(t, A=3.0, tau=5.0, t0=0.0, y0=0.0)
+        result = expFun(t, A=3.0, tau=5.0, t0=0.0)
         active = result[t >= 0]
         assert np.all(np.diff(active) <= 1e-12)
 
@@ -159,7 +201,7 @@ class TestSinFun:
     #
     def test_zero_before_t0(self):
         t = make_time_axis()
-        result = sinFun(t, A=1.0, f=0.5, phi=0.0, t0=0.0, y0=0.0)
+        result = sinFun(t, A=1.0, f=0.5, phi=0.0, t0=0.0)
         np.testing.assert_allclose(result[t < 0], 0.0)
 
     #
@@ -168,7 +210,7 @@ def test_frequency(self):
 
         t = make_time_axis()
         f = 0.5
-        result = sinFun(t, A=1.0, f=f, phi=0.0, t0=0.0, y0=0.0)
+        result = sinFun(t, A=1.0, f=f, phi=0.0, t0=0.0)
         # At t = 1/(4f), should be at maximum (A)
         idx = np.argmin(np.abs(t - 1 / (4 * f)))
         assert result[idx] == pytest.approx(1.0, abs=0.02)
@@ -178,29 +220,26 @@ def test_phase_shift(self):
         """Phi = pi/2 turns sin into cos (starts at maximum)."""
 
         t = make_time_axis()
-        result = sinFun(t, A=1.0, f=0.5, phi=np.pi / 2, t0=0.0, y0=0.0)
+        result = sinFun(t, A=1.0, f=0.5, phi=np.pi / 2, t0=0.0)
         idx_t0 = np.argmin(np.abs(t - 0.0))
         assert result[idx_t0] == pytest.approx(1.0, abs=0.02)
 
     #
     def test_zero_amplitude(self):
-        """A=0 should return y0 for t>=t0, 0 for t<t0."""
+        """A=0 gives all zeros."""
 
         t = make_time_axis()
-        result = sinFun(t, A=0.0, f=0.5, phi=0.0, t0=0.0, y0=2.0)
-        np.testing.assert_allclose(result[t < 0], 0.0)
-        np.testing.assert_allclose(result[t >= 0], 2.0)
+        result = sinFun(t, A=0.0, f=0.5, phi=0.0, t0=0.0)
+        np.testing.assert_allclose(result, 0.0, atol=1e-15)
 
     #
-    def test_offset(self):
-        """y0 shifts the oscillation center."""
+    def test_zero_at_t0(self):
+        """Phi = 0: oscillation starts at zero and centers on zero."""
 
         t = make_time_axis()
-        result = sinFun(t, A=1.0, f=0.5, phi=0.0, t0=0.0, y0=3.0)
-        # Mean of oscillation should be ~y0 over full cycles
-        # At t0, sin=0 so value = y0
+        result = sinFun(t, A=1.0, f=0.5, phi=0.0, t0=0.0)
         idx_t0 = np.argmin(np.abs(t - 0.0))
-        assert result[idx_t0] == pytest.approx(3.0, abs=0.02)
+        assert result[idx_t0] == pytest.approx(0.0, abs=0.02)
 
 
 #
@@ -209,7 +248,7 @@ class TestSinDivX:
     #
     def test_zero_before_t0(self):
         t = make_time_axis()
-        result = sinDivX(t, A=1.0, f=0.5, t0=0.0, y0=0.0)
+        result = sinDivX(t, A=1.0, f=0.5, t0=0.0)
         np.testing.assert_allclose(result[t < 0], 0.0)
 
     #
@@ -217,7 +256,7 @@ def test_decaying_envelope(self):
         """Amplitude should decrease over time (sinc envelope)."""
 
         t = make_time_axis()
-        result = sinDivX(t, A=1.0, f=0.5, t0=0.0, y0=0.0)
+        result = sinDivX(t, A=1.0, f=0.5, t0=0.0)
         active = result[t > 0.5]  # skip near t0 where sinc diverges
         peaks = np.abs(active[1:-1])[
             (active[1:-1] > active[:-2]) & (active[1:-1] > active[2:])
@@ -227,14 +266,13 @@ def test_decaying_envelope(self):
 
     #
     def test_value_at_t0(self):
-        """At t=t0, sinc(0)=1 so value should be A + y0."""
+        """At t=t0, sinc(0)=1 so value should be A."""
 
         t = make_time_axis()
         A = 1.5
-        y0 = 0.2
-        result = sinDivX(t, A=A, f=0.5, t0=0.0, y0=y0)
+        result = sinDivX(t, A=A, f=0.5, t0=0.0)
         idx_t0 = np.argmin(np.abs(t - 0.0))
-        assert result[idx_t0] == pytest.approx(A + y0, abs=0.01)
+        assert result[idx_t0] == pytest.approx(A, abs=0.01)
 
     #
     def test_first_zero_location(self):
@@ -243,27 +281,26 @@ def test_first_zero_location(self):
         t = make_time_axis()
         f = 0.5
         t0 = 0.0
-        y0 = 0.3
-        result = sinDivX(t, A=1.0, f=f, t0=t0, y0=y0)
+        result = sinDivX(t, A=1.0, f=f, t0=t0)
         t_zero = t0 + 1.0 / (2.0 * f)
         idx_zero = np.argmin(np.abs(t - t_zero))
-        assert result[idx_zero] == pytest.approx(y0, abs=0.01)
+        assert result[idx_zero] == pytest.approx(0.0, abs=0.01)
 
     #
     def test_zero_amplitude(self):
-        """A=0 with y0=0 gives all zeros."""
+        """A=0 gives all zeros."""
 
         t = make_time_axis()
-        result = sinDivX(t, A=0.0, f=0.5, t0=0.0, y0=0.0)
+        result = sinDivX(t, A=0.0, f=0.5, t0=0.0)
         np.testing.assert_allclose(result, 0.0, atol=1e-15)
 
     #
     def test_asymptotic_value(self):
-        """For t >> t0, sinc → 0 so value → y0."""
+        """For t >> t0, sinc → 0 so value → 0."""
 
         t = np.linspace(-10, 1000, 10000)
-        result = sinDivX(t, A=2.0, f=0.5, t0=0.0, y0=3.0)
-        assert result[-1] == pytest.approx(3.0, abs=0.01)
+        result = sinDivX(t, A=2.0, f=0.5, t0=0.0)
+        assert result[-1] == pytest.approx(0.0, abs=0.01)
 
 
 #
@@ -271,52 +308,52 @@ def test_asymptotic_value(self):
 class TestErfFun:
     #
     def test_midpoint_value(self):
-        """At t = t0, erf(0) = 0, so value = A/2 + y0."""
+        """At t = t0, erf(0) = 0, so value = A/2."""
 
         t = make_time_axis()
-        result = erfFun(t, A=4.0, SD=1.0, t0=10.0, y0=1.0)
+        result = erfFun(t, A=4.0, SD=1.0, t0=10.0)
         idx = np.argmin(np.abs(t - 10.0))
-        assert result[idx] == pytest.approx(3.0, abs=0.01)  # 4/2 + 1
+        assert result[idx] == pytest.approx(2.0, abs=0.01)  # A/2
 
     #
     def test_asymptotic_low(self):
-        """For t << t0, erf → -1, so value → y0."""
+        """For t << t0, erf → -1, so value → 0."""
 
         t = make_time_axis()
-        result = erfFun(t, A=4.0, SD=1.0, t0=10.0, y0=1.0)
-        assert result[0] == pytest.approx(1.0, abs=1e-6)
+        result = erfFun(t, A=4.0, SD=1.0, t0=10.0)
+        assert result[0] == pytest.approx(0.0, abs=1e-6)
 
     #
     def test_asymptotic_high(self):
-        """For t >> t0, erf → 1, so value → A + y0."""
+        """For t >> t0, erf → 1, so value → A."""
 
         t = make_time_axis()
-        result = erfFun(t, A=4.0, SD=1.0, t0=10.0, y0=1.0)
-        assert result[-1] == pytest.approx(5.0, abs=1e-3)
+        result = erfFun(t, A=4.0, SD=1.0, t0=10.0)
+        assert result[-1] == pytest.approx(4.0, abs=1e-3)
 
     #
     def test_monotonic_increase(self):
         """Error function rise should be monotonically increasing."""
 
         t = make_time_axis()
-        result = erfFun(t, A=2.0, SD=1.0, t0=10.0, y0=0.0)
+        result = erfFun(t, A=2.0, SD=1.0, t0=10.0)
         assert np.all(np.diff(result) >= -1e-12)
 
     #
     def test_zero_amplitude(self):
-        """A=0 with y0=0 gives all zeros."""
+        """A=0 gives all zeros."""
 
         t = make_time_axis()
-        result = erfFun(t, A=0.0, SD=1.0, t0=10.0, y0=0.0)
+        result = erfFun(t, A=0.0, SD=1.0, t0=10.0)
         np.testing.assert_allclose(result, 0.0, atol=1e-15)
 
     #
     def test_note_no_hard_t0_cutoff(self):
         """erfFun does NOT have a hard t0 cutoff — it's a smooth sigmoid.
-        Value at t << t0 approaches y0 but is never exactly zero."""
+        Value at t << t0 approaches zero but is never exactly zero."""
 
         t = make_time_axis()
-        result = erfFun(t, A=4.0, SD=1.0, t0=10.0, y0=0.0)
+        result = erfFun(t, A=4.0, SD=1.0, t0=10.0)
         assert result[0] == pytest.approx(0.0, abs=1e-6)
         # Slight non-zero values near t0 are expected
         idx_before = np.argmin(np.abs(t - 8.0))  # 2 SD before
@@ -328,55 +365,47 @@ def test_note_no_hard_t0_cutoff(self):
 class TestSqrtFun:
     #
     def test_zero_before_t0(self):
-        """
-        sqrtFun uses .clip(0) so t<t0 gives y0 (not zero)
-        Note sqrt(0)=0 so result = A*0 + y0).
-        """
-
         t = make_time_axis()
-        result = sqrtFun(t, A=1.0, t0=0.0, y0=0.0)
+        result = sqrtFun(t, A=1.0, t0=0.0)
         np.testing.assert_allclose(result[t < 0], 0.0, atol=1e-12)
 
+    #
+    def test_zero_before_shifted_t0(self):
+        t = make_time_axis()
+        result = sqrtFun(t, A=2.0, t0=5.0)
+        np.testing.assert_allclose(result[t < 5.0], 0.0, atol=1e-12)
+
     #
     def test_value_at_t0(self):
         t = make_time_axis()
-        result = sqrtFun(t, A=2.0, t0=0.0, y0=3.0)
+        result = sqrtFun(t, A=2.0, t0=0.0)
         idx = np.argmin(np.abs(t - 0.0))
-        assert result[idx] == pytest.approx(3.0, abs=0.01)  # A*sqrt(0) + y0
+        assert result[idx] == pytest.approx(0.0, abs=0.01)  # A*sqrt(0)
 
     #
     def test_known_value(self):
-        """At t=4 with t0=0: A*sqrt(4) + y0 = 2A + y0."""
+        """At t=4 with t0=0: A*sqrt(4) = 2A."""
 
         t = make_time_axis()
-        result = sqrtFun(t, A=3.0, t0=0.0, y0=1.0)
+        result = sqrtFun(t, A=3.0, t0=0.0)
         idx = np.argmin(np.abs(t - 4.0))
-        assert result[idx] == pytest.approx(7.0, abs=0.1)  # 3*2 + 1
+        assert result[idx] == pytest.approx(6.0, abs=0.1)  # 3*2
 
     #
     def test_monotonic_after_t0(self):
         t = make_time_axis()
-        result = sqrtFun(t, A=1.0, t0=0.0, y0=0.0)
+        result = sqrtFun(t, A=1.0, t0=0.0)
         active = result[t >= 0]
         assert np.all(np.diff(active) >= -1e-12)
 
     #
     def test_zero_amplitude(self):
-        """A=0 with y0=0 gives all zeros."""
+        """A=0 gives all zeros."""
 
         t = make_time_axis()
-        result = sqrtFun(t, A=0.0, t0=0.0, y0=0.0)
+        result = sqrtFun(t, A=0.0, t0=0.0)
         np.testing.assert_allclose(result, 0.0, atol=1e-15)
 
-    #
-    def test_offset_y0(self):
-        """y0 shifts everything (including before t0 via clip behavior)."""
-
-        t = make_time_axis()
-        result = sqrtFun(t, A=1.0, t0=0.0, y0=5.0)
-        # Before t0: A*sqrt(0) + y0 = y0
-        np.testing.assert_allclose(result[t <= 0], 5.0, atol=0.01)
-
 
 if __name__ == "__main__":
     pytest.main([__file__])
diff --git a/tests/test_graph_ir.py b/tests/test_graph_ir.py
index eeff438..060fdf8 100644
--- a/tests/test_graph_ir.py
+++ b/tests/test_graph_ir.py
@@ -452,15 +452,15 @@ def test_dynamics_param_nodes(self):
         )
         graph = build_graph(model)
 
-        # expFun has params: A, tau, t0, y0
-        # A and tau are vary=True, t0 and y0 are vary=False
+        # expFun has params: A, tau, t0
+        # A and tau are vary=True, t0 is vary=False
         dyn_A = _node_by_name(graph, "GLP_01_A_expFun_01_A")
         assert dyn_A is not None
         assert dyn_A.kind == NodeKind.OPT_PARAM
 
-        dyn_y0 = _node_by_name(graph, "GLP_01_A_expFun_01_y0")
-        assert dyn_y0 is not None
-        assert dyn_y0.kind == NodeKind.STATIC_PARAM
+        dyn_t0 = _node_by_name(graph, "GLP_01_A_expFun_01_t0")
+        assert dyn_t0 is not None
+        assert dyn_t0.kind == NodeKind.STATIC_PARAM
 
     #
     def test_dynamics_edges(self):
@@ -477,7 +477,7 @@ def test_dynamics_edges(self):
         assert len(trace_nodes) >= 1
         trace_nid = trace_nodes[0].id
         param_edges = _edges_to(graph, trace_nid, EdgeKind.PARAM_INPUT)
-        assert len(param_edges) == 4  # A, tau, t0, y0
+        assert len(param_edges) == 3  # A, tau, t0
 
         # PARAM_PLUS_TRACE has BASE_INPUT + TRACE_INPUT
         resolved = _node_by_name(graph, "GLP_01_A_resolved")
@@ -1307,16 +1307,16 @@ def test_profile_par_dynamics_params_in_graph(self):
         _file, model = _make_time_dep_profile_model()
         graph = build_graph(model)
 
-        # MonoExpPos has params: A, tau, t0, y0
+        # MonoExpPos has params: A, tau, t0
         # The profile par is GLP_01_A_pLinear_01_m, so dynamics params
         # are prefixed: GLP_01_A_pLinear_01_m_expFun_01_A, etc.
         dyn_A = _node_by_name(graph, "GLP_01_A_pLinear_01_m_expFun_01_A")
         assert dyn_A is not None
         assert dyn_A.kind == NodeKind.OPT_PARAM
 
-        dyn_y0 = _node_by_name(graph, "GLP_01_A_pLinear_01_m_expFun_01_y0")
-        assert dyn_y0 is not None
-        assert dyn_y0.kind == NodeKind.STATIC_PARAM
+        dyn_t0 = _node_by_name(graph, "GLP_01_A_pLinear_01_m_expFun_01_t0")
+        assert dyn_t0 is not None
+        assert dyn_t0.kind == NodeKind.STATIC_PARAM
 
     #
     def test_resolved_profile_par_in_sample_edges(self):
@@ -1992,7 +1992,7 @@ def test_dynamics_compiled(self):
         # Single substep in the group
         assert plan.dyn_group_indptr[1] - plan.dyn_group_indptr[0] == 1
         assert plan.dyn_sub_func_id[0] == int(DynFuncKind.EXPFUN)
-        assert plan.dyn_sub_n_params[0] == 4  # A, tau, t0, y0
+        assert plan.dyn_sub_n_params[0] == 3  # A, tau, t0
 
     #
     def test_dynamics_param_rows_valid(self):
@@ -2610,13 +2610,6 @@ def test_synthetic_graph_with_prefixed_expr_ref(self):
                 source_order=3,
                 value=0.0,
             ),
-            GraphNode(
-                id=4,
-                kind=NodeKind.STATIC_PARAM,
-                name="dyn_y0",
-                source_order=4,
-                value=0.0,
-            ),
             GraphNode(
                 id=5,
                 kind=NodeKind.DYNAMICS_TRACE,
@@ -2673,7 +2666,6 @@ def test_synthetic_graph_with_prefixed_expr_ref(self):
             GraphEdge(source=1, target=5, kind=EdgeKind.PARAM_INPUT, position=0),
             GraphEdge(source=2, target=5, kind=EdgeKind.PARAM_INPUT, position=1),
             GraphEdge(source=3, target=5, kind=EdgeKind.PARAM_INPUT, position=2),
-            GraphEdge(source=4, target=5, kind=EdgeKind.PARAM_INPUT, position=3),
             GraphEdge(source=0, target=6, kind=EdgeKind.BASE_INPUT),
             GraphEdge(source=5, target=6, kind=EdgeKind.TRACE_INPUT),
             # EXPR_REF: B_expr references A_base_resolved
@@ -2861,13 +2853,6 @@ def test_sparse_node_ids(self):
                 source_order=3,
                 value=0.0,
             ),
-            GraphNode(
-                id=104,
-                kind=NodeKind.STATIC_PARAM,
-                name="dyn_y0",
-                source_order=4,
-                value=0.0,
-            ),
             GraphNode(
                 id=105,
                 kind=NodeKind.DYNAMICS_TRACE,
@@ -2915,7 +2900,6 @@ def test_sparse_node_ids(self):
             GraphEdge(source=101, target=105, kind=EdgeKind.PARAM_INPUT, position=0),
             GraphEdge(source=102, target=105, kind=EdgeKind.PARAM_INPUT, position=1),
             GraphEdge(source=103, target=105, kind=EdgeKind.PARAM_INPUT, position=2),
-            GraphEdge(source=104, target=105, kind=EdgeKind.PARAM_INPUT, position=3),
             GraphEdge(source=100, target=106, kind=EdgeKind.BASE_INPUT),
             GraphEdge(source=105, target=106, kind=EdgeKind.TRACE_INPUT),
             GraphEdge(source=106, target=109, kind=EdgeKind.PARAM_INPUT, position=0),
diff --git a/tests/test_mcp_library.py b/tests/test_mcp_library.py
index caeeab5..b652557 100644
--- a/tests/test_mcp_library.py
+++ b/tests/test_mcp_library.py
@@ -279,7 +279,6 @@ def test_dynamics_model_creation(self):
                 "A": [2, True, 1, 1e2],
                 "tau": [5000, True, 1e3, 1e4],
                 "t0": [0, False, 0, 1],
-                "y0": [0, False, 0, 1],
             }
         )
 
@@ -289,7 +288,6 @@ def test_dynamics_model_creation(self):
                 "A": [5, True, 1, 1e2],
                 "tau": [1250, True, 1e2, 1e3],
                 "t0": [0, False, 0, 1],
-                "y0": [0, False, 0, 1],
             }
         )
 
@@ -318,7 +316,6 @@ def test_dynamics_parameter_handling(self):
                 "A": [1, True, 0, 5],
                 "tau": [2.5, True, 1, 10],
                 "t0": [0, False, 0, 1],
-                "y0": [0, False, 0, 1],
             }
         )
 
diff --git a/tests/test_model_parser.py b/tests/test_model_parser.py
index a363111..b9c5475 100644
--- a/tests/test_model_parser.py
+++ b/tests/test_model_parser.py
@@ -215,7 +215,6 @@ def test_simple_time_model(self):
         assert model.components[0].par_dict["A"] == [1, True, 0, 5]
         assert model.components[0].par_dict["tau"] == [2.5, True, 1, 10]
         assert model.components[0].par_dict["t0"] == [0, False, 0, 1]
-        assert model.components[0].par_dict["y0"] == [0, False, 0, 1]
 
     #
     def test_IRF_model(self):
@@ -239,7 +238,6 @@ def test_IRF_model(self):
         assert model.components[1].par_dict["A"] == [1, True, 0, 5]
         assert model.components[1].par_dict["tau"] == [2.5, True, 1, 10]
         assert model.components[1].par_dict["t0"] == [0, False, 0, 1]
-        assert model.components[1].par_dict["y0"] == [0, False, 0, 1]
 
     #
     def test_multi_cycle_expression_model(self):
@@ -271,7 +269,6 @@ def test_multi_cycle_expression_model(self):
         assert model.components[2].par_dict["A"] == ["-expFun_01_A"]
         assert model.components[2].par_dict["tau"] == ["expFun_01_tau"]
         assert model.components[2].par_dict["t0"] == [0, False, 0, 1]
-        assert model.components[2].par_dict["y0"] == [0, False, 0, 1]
 
         # Check lmfit parameters exist with prefixed names
         assert "parTEST_expFun_01_A" in model.lmfit_pars
@@ -351,7 +348,6 @@ def test_simple_2D_model(self):
         assert td_par_model.components[1].par_dict["A"] == [1, True, 0, 5]
         assert td_par_model.components[1].par_dict["tau"] == [2.5, True, 1, 10]
         assert td_par_model.components[1].par_dict["t0"] == [0, False, 0, 1]
-        assert td_par_model.components[1].par_dict["y0"] == [0, False, 0, 1]
         # end of time-dependent parameter model
         assert model.components[2].par_dict["F"] == [1.0, True, 0.75, 2.5]
         assert model.components[2].par_dict["m"] == [0.3, True, 0, 1]
diff --git a/tests/test_project_fit.py b/tests/test_project_fit.py
index bca00a0..e04e0b9 100644
--- a/tests/test_project_fit.py
+++ b/tests/test_project_fit.py
@@ -53,7 +53,6 @@ def _make_truth_file(*, amplitude=20.0, x0_shift=3.0, tau=5.0):
     model.lmfit_pars["GLP_01_x0_expFun_01_A"].value = x0_shift
     model.lmfit_pars["GLP_01_x0_expFun_01_tau"].value = tau
     model.lmfit_pars["GLP_01_x0_expFun_01_t0"].value = 0.0
-    model.lmfit_pars["GLP_01_x0_expFun_01_y0"].value = 0.0
 
     return file
 
@@ -144,11 +143,9 @@ def test_biexp_expr_t0_roundtrip(self):
             m.lmfit_pars["GLP_01_x0_expFun_01_A"].value = dx0_1
             m.lmfit_pars["GLP_01_x0_expFun_01_tau"].value = TRUE_TAU1
             m.lmfit_pars["GLP_01_x0_expFun_01_t0"].value = TRUE_T0
-            m.lmfit_pars["GLP_01_x0_expFun_01_y0"].value = 0.0
             m.lmfit_pars["GLP_01_x0_expFun_02_A"].value = dx0_2
             m.lmfit_pars["GLP_01_x0_expFun_02_tau"].value = TRUE_TAU2
             m.lmfit_pars["GLP_01_x0_expFun_02_t0"].value = TRUE_T0
-            m.lmfit_pars["GLP_01_x0_expFun_02_y0"].value = 0.0
             truth_files.append((tf, seed))
 
         # --- simulate and build fit files ---
@@ -276,7 +273,6 @@ def test_vary_levels_map(self):
         assert levels["GLP_01_x0_expFun_01_A"] == "file"
         assert levels["GLP_01_x0_expFun_01_tau"] == "project"
         assert levels["GLP_01_x0_expFun_01_t0"] == "static"
-        assert levels["GLP_01_x0_expFun_01_y0"] == "static"
 
     #
     def test_vary_levels_profile_with_dynamics(self):