diff --git a/src/easyscience/fitting/minimizers/minimizer_bumps.py b/src/easyscience/fitting/minimizers/minimizer_bumps.py
index 315d856c..500edce6 100644
--- a/src/easyscience/fitting/minimizers/minimizer_bumps.py
+++ b/src/easyscience/fitting/minimizers/minimizer_bumps.py
@@ -417,42 +417,45 @@ def mcmc_sample(
         weights : np.ndarray
             Flattened weight array.
         samples : int, default=10000
-            Number of raw samples to draw across all chains, before thinning.
-            A guaranteed minimum, not an exact count: DREAM advances in
-            blocks of 10 generations (one generation = one draw per chain)
-            and stops at the first block boundary at or past ``samples``.
+            Number of raw samples to draw across all chains, before
+            thinning. A guaranteed minimum, not an exact count: DREAM
+            advances in blocks of 10 generations (one generation = one
+            draw per chain) and stops at the first block boundary at or
+            past ``samples``.
         burn : int, default=2000
             Burn-in generations to discard. BUMPS counts ``burn`` in
-            generations while ``samples`` counts raw draws, so ``burn=500``
-            discards ``500 * n_chains`` raw samples.
+            generations while ``samples`` counts raw draws, so
+            ``burn=500`` discards ``500 * n_chains`` raw samples.
         thin : int, default=10
-            Thinning interval — only every ``thin``-th generation is stored.
+            Thinning interval — only every ``thin``-th generation is
+            stored.
         population : int | None, default=None
-            BUMPS DREAM population count per parameter (number of parallel
-            chains): BUMPS creates ``ceil(population * n_parameters)`` chains.
+            BUMPS DREAM population count per parameter (number of
+            parallel chains): BUMPS creates ``ceil(population *
+            n_parameters)`` chains.
         resume_state : MCMCDraw | None, default=None
             A BUMPS ``MCMCDraw`` state object from a previous
-            ``mcmc_sample()`` call (e.g. ``PosteriorResults.sampler_state``).
-            When provided, DREAM **continues** the saved chain instead of
-            starting cold.  The population, parameter count, and parameter
-            names must match the current model — a ``ValueError`` is raised
-            otherwise.
-
-            ``samples`` must be the **total** number of raw samples, not an
-            increment: to extend an existing chain of ``N`` raw samples by
-            ``M``, pass ``samples=N + M`` (DREAM keeps only the last
-            ``samples`` draws in its buffer). The `Sampler.extend` helper
-            computes this for you.
-
-            ``burn`` is forced to 0 on resume: a previously-converged chain is
-            never re-burned.
-
-            The ``population`` and ``initializer`` parameters
-            have **no effect** when ``resume_state`` is provided — they
-            are determined by the saved state.
-
-            Resuming against *different* data is undefined behaviour (the
-            chain's likelihood changes underneath it).
+            ``mcmc_sample()`` call (e.g.
+            ``PosteriorResults.sampler_state``). When provided, DREAM
+            **continues** the saved chain instead of starting cold.  The
+            population, parameter count, and parameter names must match
+            the current model — a ``ValueError`` is raised otherwise.
+
+            ``samples`` must be the **total** number of raw samples, not
+            an increment: to extend an existing chain of ``N`` raw
+            samples by ``M``, pass ``samples=N + M`` (DREAM keeps only
+            the last ``samples`` draws in its buffer). The
+            ``Sampler.extend`` helper computes this for you.
+
+            ``burn`` is forced to 0 on resume: a previously-converged
+            chain is never re-burned.
+
+            The ``population`` and ``initializer`` parameters have **no
+            effect** when ``resume_state`` is provided — they are
+            determined by the saved state.
+
+            Resuming against *different* data is undefined behaviour
+            (the chain's likelihood changes underneath it).
         sampler_kwargs : dict | None, default=None
             Additional keyword arguments forwarded to
             ``bumps.fitters.fit``.
@@ -475,9 +478,9 @@ def mcmc_sample(
         ------
         ValueError
             If the input shapes or weights are invalid, if
-            ``progress_callback`` is not callable, or if ``resume_state``
-            is incompatible with the current model (parameter count,
-            names/order, or population mismatch).
+            ``progress_callback`` is not callable, or if
+            ``resume_state`` is incompatible with the current model
+            (parameter count, names/order, or population mismatch).
         FitError
             If DREAM sampling was aborted by the user (via
             ``abort_test``).
@@ -608,28 +611,32 @@ def _validate_resume_state(
         population: int | None,
         burn: int,
     ) -> tuple[int, int]:
-        """Check that ``resume_state`` is compatible with ``problem`` and
+        """
+        Check that ``resume_state`` is compatible with ``problem`` and
         resolve the population and burn values to use when resuming.
 
         Parameters
         ----------
         problem : FitProblem
-            The freshly built BUMPS ``FitProblem`` for the current model.
+            The freshly built BUMPS ``FitProblem`` for the current
+            model.
         resume_state : MCMCDraw
             The saved chain state to resume from.
         population : int | None
             The caller-supplied population scale factor, or ``None``.
         burn : int
-            The caller-supplied burn-in, ignored (with a warning) on resume.
+            The caller-supplied burn-in, ignored (with a warning) on
+            resume.
 
         Returns
         -------
         tuple[int, int]
             ``(population, burn)`` to pass to DREAM. The population is
             returned as a **negative** number, which BUMPS'
-            ``initpop.generate`` reads as an absolute chain count, exactly
-            reproducing the saved state's population. ``burn`` is always 0:
-            a previously converged chain is never re-burned.
+            ``initpop.generate`` reads as an absolute chain count,
+            exactly reproducing the saved state's population. ``burn``
+            is always 0: a previously converged chain is never
+            re-burned.
 
         Raises
         ------
diff --git a/src/easyscience/fitting/sampler.py b/src/easyscience/fitting/sampler.py
index d3915115..3917a5bc 100644
--- a/src/easyscience/fitting/sampler.py
+++ b/src/easyscience/fitting/sampler.py
@@ -1,6 +1,8 @@
 # SPDX-FileCopyrightText: 2026 EasyScience contributors <https://github.com/easyscience>
 # SPDX-License-Identifier: BSD-3-Clause
-"""Bayesian MCMC sampling — the ``Sampler`` class and persistence helpers."""
+"""
+Bayesian MCMC sampling — the ``Sampler`` class and persistence helpers.
+"""
 
 from __future__ import annotations
 
@@ -38,7 +40,9 @@ def _data_fingerprint(
     y_list: list,
     w_list: list,
 ) -> str | None:
-    """Return a SHA-256 hex digest of concatenated (x|y|weights), or None."""
+    """
+    Return a SHA-256 hex digest of concatenated (x|y|weights), or None.
+    """
     try:
         h = hashlib.sha256()
         for arr in list(x_list) + list(y_list) + list(w_list):
@@ -53,7 +57,8 @@ def _validate_dataset_arrays(
     data: np.ndarray | list | tuple,
     allow_none_entries: bool = False,
 ) -> None:
-    """Check that ``data`` (an array or list of arrays) holds numeric,
+    """
+    Check that ``data`` (an array or list of arrays) holds numeric,
     at-least-1-D, non-empty arrays.
 
     Structural checks (array vs list, matching dataset counts) are in
@@ -101,12 +106,13 @@ def _validate_dataset_arrays(
 
 
 def _copy_data(data):
-    """Copy an array (or list of arrays) and mark the copies read-only.
+    """
+    Copy an array (or list of arrays) and mark the copies read-only.
 
-    The ``Sampler`` binds its data at construction; copying decouples the
-    bound data from the caller's arrays, and the read-only flag stops
-    in-place mutation of the copies, so the chain and the ``save()``
-    fingerprint always describe the data actually sampled.
+    The ``Sampler`` binds its data at construction; copying decouples
+    the bound data from the caller's arrays, and the read-only flag
+    stops in-place mutation of the copies, so the chain and the
+    ``save()`` fingerprint always describe the data actually sampled.
     """
     if data is None:
         return None
@@ -118,10 +124,13 @@ def _copy_data(data):
 
 
 def _validate_chain_path(path: str | os.PathLike, skip: int = 0) -> str:
-    """Validate the persistence arguments shared by the save/load functions.
+    """
+    Validate the persistence arguments shared by the save/load
+    functions.
 
-    Returns ``path`` coerced to ``str``. Raises ``TypeError`` if ``path`` is
-    not path-like and ``ValueError`` if ``skip`` is not a non-negative int.
+    Returns ``path`` coerced to ``str``. Raises ``TypeError`` if
+    ``path`` is not path-like and ``ValueError`` if ``skip`` is not a
+    non-negative int.
     """
     if not isinstance(path, (str, os.PathLike)):
         raise TypeError(f'path must be a str or os.PathLike, got {type(path).__name__}.')
@@ -131,14 +140,15 @@ def _validate_chain_path(path: str | os.PathLike, skip: int = 0) -> str:
 
 
 def _load_bumps_state(path: str, skip: int = 0) -> MCMCDraw:
-    """Read a BUMPS chain from disk, working around a BUMPS bug.
+    """
+    Read a BUMPS chain from disk, working around a BUMPS bug.
 
     ``bumps.dream.state.load_state`` reads the saved buffers with
-    ``numpy.loadtxt``, which collapses a single-row file to a 1-D array. A
-    short chain stores a single CR-weight update row, so ``load_state``'s
-    subsequent ``stats[:, 0]`` indexing raises ``IndexError: too many indices
-    for array``. We coerce each buffer read back to 2-D before ``load_state``
-    consumes it.
+    ``numpy.loadtxt``, which collapses a single-row file to a 1-D array.
+    A short chain stores a single CR-weight update row, so
+    ``load_state``'s subsequent ``stats[:, 0]`` indexing raises
+    ``IndexError: too many indices for array``. We coerce each buffer
+    read back to 2-D before ``load_state`` consumes it.
 
     Parameters
     ----------
@@ -170,28 +180,30 @@ def _loadtxt_2d(file, report: int = 0) -> np.ndarray:
 
 
 def load_chain(path: str | os.PathLike, skip: int = 0) -> tuple[MCMCDraw, list[str] | None, dict]:
-    """Reload a DREAM chain state saved by ``Sampler.save``.
+    """
+    Reload a DREAM chain state saved by ``Sampler.save``.
 
-    This is the standalone reader: unlike ``Sampler.load_state`` it needs no
-    fitter, model or data, so a saved chain can be inspected or post-processed
-    on a machine that does not have the model. Parameter names are restored
-    from the sidecar when available (schema versions 1 and 2), falling back to
-    the state's labels with the minimizer prefix stripped.
+    This is the standalone reader: unlike ``Sampler.load_state`` it
+    needs no fitter, model or data, so a saved chain can be inspected or
+    post-processed on a machine that does not have the model. Parameter
+    names are restored from the sidecar when available (schema versions
+    1 and 2), falling back to the state's labels with the minimizer
+    prefix stripped.
 
     Parameters
     ----------
     path : str | os.PathLike
         File path prefix used when saving.
     skip : int, default=0
-        Discard the first ``skip`` saved generations on load, forwarded to
-        ``bumps.dream.state.load_state``.
+        Discard the first ``skip`` saved generations on load, forwarded
+        to ``bumps.dream.state.load_state``.
 
     Returns
     -------
     tuple[MCMCDraw, list[str] | None, dict]
-        The reloaded BUMPS chain state, the parameter names (or ``None`` if
-        neither sidecar nor labels yielded them), and the raw sidecar dict
-        (empty if absent/unreadable).
+        The reloaded BUMPS chain state, the parameter names (or ``None``
+        if neither sidecar nor labels yielded them), and the raw sidecar
+        dict (empty if absent/unreadable).
 
     Raises
     ------
@@ -229,15 +241,18 @@ def load_chain(path: str | os.PathLike, skip: int = 0) -> tuple[MCMCDraw, list[s
 
 @dataclass
 class SamplingResults:
-    """Structured result of an MCMC sampling run (analogous to ``FitResults``).
+    """
+    Structured result of an MCMC sampling run (analogous to
+    ``FitResults``).
 
     Attributes
     ----------
     draws : np.ndarray
-        Posterior samples, shape ``(n_samples, n_params)``. For results from
-        ``Sampler.sample()``/``Sampler.extend()`` this is a *trimmed,
-        outlier-filtered* view of the chain, not the raw buffer, so it holds
-        fewer rows than ``samples / thin`` — see the notes on ``Sampler``.
+        Posterior samples, shape ``(n_samples, n_params)``. For results
+        from ``Sampler.sample()``/``Sampler.extend()`` this is a
+        *trimmed, outlier-filtered* view of the chain, not the raw
+        buffer, so it holds fewer rows than ``samples / thin`` — see the
+        notes on ``Sampler``.
     param_names : list[str]
         Parameter names (one per column of ``draws``).
     logp : np.ndarray
@@ -252,8 +267,10 @@ class SamplingResults:
     state: MCMCDraw
 
     def to_legacy_dict(self) -> dict:
-        """Return the legacy dict shape produced by the deprecated
-        ``mcmc_sample()`` APIs."""
+        """
+        Return the legacy dict shape produced by the deprecated
+        ``mcmc_sample()`` APIs.
+        """
         return {
             'draws': self.draws,
             'param_names': self.param_names,
@@ -263,58 +280,64 @@ def to_legacy_dict(self) -> dict:
 
 
 class Sampler:
-    """Bayesian MCMC sampler for one dataset, backed by a Fitter's BUMPS minimizer.
-
-    One ``Sampler`` instance represents one chain over one ``(x, y, weights)``
-    dataset. The data is bound at construction; ``sample()`` and ``extend()``
-    take no data arguments, so a chain can never be extended against different
-    data (undefined behaviour in BUMPS). The bound data is a defensive,
-    read-only copy of the caller's arrays, exposed via the ``x``, ``y`` and
-    ``weights`` properties — mutating the originals after construction has no
-    effect on the sampler, and there are deliberately no setters: to sample
-    different data, create a new ``Sampler``.
-
-    Construct directly with a configured ``Fitter`` (or ``MultiFitter``) whose
-    minimizer has been switched to ``AvailableMinimizers.Bumps``. **Running a
-    fit first is not required** — the ``Fitter`` supplies the model and the
-    minimizer, not a fit result, and sampling from the initial parameter values
-    works fine.
-
-    It is often worth fitting first anyway. DREAM seeds its whole starting
-    population inside a tiny ball around the parameters' *current* values
-    (BUMPS' default ``init='eps'``), so sampling from fitted values starts the
-    chain in the right region and shortens the burn-in needed to reach the
-    typical set. From a poor initial guess, expect to burn for longer.
-
-    The sampler is BUMPS/DREAM-specific for now: the BUMPS check in ``_run()``
-    is the seam where another backend would plug in.
+    """
+    Bayesian MCMC sampler for one dataset, backed by a Fitter's BUMPS
+    minimizer.
+
+    One ``Sampler`` instance represents one chain over one ``(x, y,
+    weights)`` dataset. The data is bound at construction; ``sample()``
+    and ``extend()`` take no data arguments, so a chain can never be
+    extended against different data (undefined behaviour in BUMPS). The
+    bound data is a defensive, read-only copy of the caller's arrays,
+    exposed via the ``x``, ``y`` and ``weights`` properties — mutating
+    the originals after construction has no effect on the sampler, and
+    there are deliberately no setters: to sample different data, create
+    a new ``Sampler``.
+
+    Construct directly with a configured ``Fitter`` (or ``MultiFitter``)
+    whose minimizer has been switched to ``AvailableMinimizers.Bumps``.
+    **Running a fit first is not required** — the ``Fitter`` supplies
+    the model and the minimizer, not a fit result, and sampling from the
+    initial parameter values works fine.
+
+    It is often worth fitting first anyway. DREAM seeds its whole
+    starting population inside a tiny ball around the parameters'
+    *current* values (BUMPS' default ``init='eps'``), so sampling from
+    fitted values starts the chain in the right region and shortens the
+    burn-in needed to reach the typical set. From a poor initial guess,
+    expect to burn for longer.
+
+    The sampler is BUMPS/DREAM-specific for now: the BUMPS check in
+    ``_run()`` is the seam where another backend would plug in.
 
     Parameters
     ----------
-    fitter : Fitter
-        A configured ``Fitter`` (or ``MultiFitter``) whose minimizer has been
-        switched to ``AvailableMinimizers.Bumps``.
+    fitter : 'Fitter'
+        A configured ``Fitter`` (or ``MultiFitter``) whose minimizer has
+        been switched to ``AvailableMinimizers.Bumps``.
     x : np.ndarray | list[np.ndarray]
-        Independent variable array (or list of arrays for ``MultiFitter``).
+        Independent variable array (or list of arrays for
+        ``MultiFitter``).
     y : np.ndarray | list[np.ndarray]
-        Dependent variable array (or list of arrays for ``MultiFitter``).
+        Dependent variable array (or list of arrays for
+        ``MultiFitter``).
     weights : np.ndarray | list[np.ndarray | None] | None, default=None
         Weight array (or list of arrays for ``MultiFitter``).
     vectorized : bool, default=False
         When ``True``, each x array may be multi-dimensional (e.g. an
         ``(N, M, 2)`` grid for a 2D model) and is left as-is.
     sampler_kwargs : dict | None, default=None
-        Per-instance default keyword arguments forwarded to the BUMPS DREAM
-        sampler on every run, merged with (and overridden by) per-call
-        ``sampler_kwargs``.
+        Per-instance default keyword arguments forwarded to the BUMPS
+        DREAM sampler on every run, merged with (and overridden by)
+        per-call ``sampler_kwargs``.
 
     Raises
     ------
     TypeError
-        If ``fitter`` is not Fitter-shaped (no ``minimizer``/``fit_function``),
-        if any dataset in ``x``/``y``/``weights`` is not a numeric array
-        (e.g. a string), or ``vectorized``/``sampler_kwargs`` have the wrong
-        type.
+        If ``fitter`` is not Fitter-shaped (no
+        ``minimizer``/``fit_function``), if any dataset in
+        ``x``/``y``/``weights`` is not a numeric array (e.g. a string),
+        or ``vectorized``/``sampler_kwargs`` have the wrong type.
     ValueError
         If ``x``, ``y`` and ``weights`` do not hold matching structures
         (all arrays, or lists of the same length), or any dataset is a
@@ -322,42 +345,44 @@ class Sampler:
 
     Notes
     -----
-    **The retained draws are a trimmed view, not the whole chain.** BUMPS'
-    DREAM sampler defaults to ``trim=True``: once sampling finishes it runs a
-    convergence-based burn-point detector over the chain and returns only the
-    portion after that point. ``state.draw()`` additionally drops chains
-    flagged as outliers. So ``results.draws`` is usually *smaller* than the
-    chain, and its length is not deterministic — the detector re-runs from
-    scratch on every ``sample()`` and ``extend()`` call and may place the burn
-    point differently each time. Read the count off the array; do not predict
-    it.
-
-    Nor is the chain itself exactly ``samples / thin`` rows: ``samples`` is a
-    guaranteed minimum, not an exact count. DREAM advances in blocks of 10
-    generations (one generation = one draw per chain) and only checks its
-    stopping condition between blocks, so the raw chain length is ``samples``
-    rounded up to a multiple of ``10 * n_chains``.
-
-    The trimming is only a *view* — nothing is dropped from the buffer. To
-    read the full untrimmed chain::
+    **The retained draws are a trimmed view, not the whole chain.**
+    BUMPS' DREAM sampler defaults to ``trim=True``: once sampling
+    finishes it runs a convergence-based burn-point detector over the
+    chain and returns only the portion after that point.
+    ``state.draw()`` additionally drops chains flagged as outliers. So
+    ``results.draws`` is usually *smaller* than the chain, and its
+    length is not deterministic — the detector re-runs from scratch on
+    every ``sample()`` and ``extend()`` call and may place the burn
+    point differently each time. Read the count off the array; do not
+    predict it.
+
+    Nor is the chain itself exactly ``samples / thin`` rows: ``samples``
+    is a guaranteed minimum, not an exact count. DREAM advances in
+    blocks of 10 generations (one generation = one draw per chain) and
+    only checks its stopping condition between blocks, so the raw chain
+    length is ``samples`` rounded up to a multiple of ``10 * n_chains``.
+
+    The trimming is only a *view* — nothing is dropped from the buffer.
+    To read the full untrimmed chain::
 
         results = sampler.sample(samples=10000, burn=500, thin=2)
         full = results.state.draw(portion=1.0, outliers=True)
         full.points  # (n_chain_rows, n_params)
         full.logp
 
-    To switch the automatic trimming off entirely, pass BUMPS' own ``trim``
-    option straight through ``sampler_kwargs``; ``results.draws`` then holds
-    the whole chain::
+    To switch the automatic trimming off entirely, pass BUMPS' own
+    ``trim`` option straight through ``sampler_kwargs``;
+    ``results.draws`` then holds the whole chain::
 
         sampler = Sampler(
             fitter, x, y, weights=w, sampler_kwargs={'trim': False}
         )
 
-    Note also that trimming does not survive a ``save()``/``load_state()``
-    round-trip: BUMPS does not persist the trim point, so a reloaded chain
-    comes back untrimmed and ``load_state()`` reports *more* draws than the
-    ``sample()`` call that created it. The chain itself is identical.
+    Note also that trimming does not survive a
+    ``save()``/``load_state()`` round-trip: BUMPS does not persist the
+    trim point, so a reloaded chain comes back untrimmed and
+    ``load_state()`` reports *more* draws than the ``sample()`` call
+    that created it. The chain itself is identical.
     """
 
     def __init__(
@@ -435,7 +460,9 @@ def weights(self) -> np.ndarray | list[np.ndarray | None] | None:
 
     @property
     def state(self) -> MCMCDraw | None:
-        """Raw BUMPS MCMCDraw state (None before first sample/load_state)."""
+        """
+        Raw BUMPS MCMCDraw state (None before first sample/load_state).
+        """
         return self._state
 
     @property
@@ -459,7 +486,9 @@ def logp(self) -> np.ndarray | None:
         return self._results.logp if self._results is not None else None
 
     def _fingerprint(self) -> str | None:
-        """SHA-256 fingerprint of the bound (x, y, weights) data, or None."""
+        """
+        SHA-256 fingerprint of the bound (x, y, weights) data, or None.
+        """
         x_list = list(self._x) if isinstance(self._x, (list, tuple)) else [self._x]
         y_list = list(self._y) if isinstance(self._y, (list, tuple)) else [self._y]
         if self._weights is None:
@@ -481,7 +510,8 @@ def _run(
         progress_callback: Callable[[dict], bool | None] | None,
         abort_test: Callable[[], bool] | None,
     ) -> SamplingResults:
-        """Shared sampling engine for ``sample()`` and ``extend()``.
+        """
+        Shared sampling engine for ``sample()`` and ``extend()``.
 
         Argument validation for ``samples``/``burn``/``thin`` lives in
         ``Bumps.mcmc_sample`` (single source of truth).
@@ -545,37 +575,43 @@ def sample(
         progress_callback: Callable[[dict], bool | None] | None = None,
         abort_test: Callable[[], bool] | None = None,
     ) -> SamplingResults:
-        """Run fresh Bayesian MCMC sampling on the bound data.
+        """
+        Run fresh Bayesian MCMC sampling on the bound data.
 
-        Calling ``sample()`` on a sampler that already holds a chain starts a
-        **fresh** chain — the previous state and results are replaced. Use
-        ``extend()`` to continue an existing chain.
+        Calling ``sample()`` on a sampler that already holds a chain
+        starts a **fresh** chain — the previous state and results are
+        replaced. Use ``extend()`` to continue an existing chain.
 
         Parameters
         ----------
         samples : int, default=10000
-            Number of raw samples to draw across all chains, before thinning.
-            This is a guaranteed minimum, not an exact count: DREAM advances
-            in blocks of 10 generations (one generation = one draw per chain)
-            and stops at the first block boundary at or past ``samples``.
+            Number of raw samples to draw across all chains, before
+            thinning. This is a guaranteed minimum, not an exact count:
+            DREAM advances in blocks of 10 generations (one generation =
+            one draw per chain) and stops at the first block boundary at
+            or past ``samples``.
         burn : int, default=2000
-            Burn-in generations to discard before collecting samples. Note
-            BUMPS counts ``burn`` in generations while ``samples`` counts raw
-            draws, so ``burn=500`` discards ``500 * n_chains`` raw samples.
+            Burn-in generations to discard before collecting samples.
+            Note BUMPS counts ``burn`` in generations while ``samples``
+            counts raw draws, so ``burn=500`` discards ``500 *
+            n_chains`` raw samples.
         thin : int, default=10
-            Thinning interval — only every ``thin``-th generation is kept,
-            which reduces autocorrelation between consecutive draws.
+            Thinning interval — only every ``thin``-th generation is
+            kept, which reduces autocorrelation between consecutive
+            draws.
         population : int | None, default=None
-            DREAM population **scale factor** (not an absolute chain count):
-            BUMPS creates ``ceil(population * n_parameters)`` parallel chains.
+            DREAM population **scale factor** (not an absolute chain
+            count): BUMPS creates ``ceil(population * n_parameters)``
+            parallel chains.
         sampler_kwargs : dict | None, default=None
             Additional keyword arguments forwarded to the BUMPS DREAM
             sampler (merged over the instance defaults).
         progress_callback : Callable[[dict], bool | None] | None, default=None
-            Optional callback invoked at each DREAM generation. The payload
-            dict includes ``iteration`` and ``sampling: True``.
+            Optional callback invoked at each DREAM generation. The
+            payload dict includes ``iteration`` and ``sampling: True``.
         abort_test : Callable[[], bool] | None, default=None
-            Optional callable that returns ``True`` to abort sampling early.
+            Optional callable that returns ``True`` to abort sampling
+            early.
 
         Returns
         -------
@@ -584,14 +620,16 @@ def sample(
 
         Notes
         -----
-        ``results.draws`` is a trimmed, outlier-filtered view of the chain,
-        so it is usually smaller than the chain and its length is not
-        predictable from ``samples``/``thin``. See the ``Sampler`` class notes
-        for how to read the full chain or switch trimming off.
+        ``results.draws`` is a trimmed, outlier-filtered view of the
+        chain, so it is usually smaller than the chain and its length is
+        not predictable from ``samples``/``thin``. See the ``Sampler``
+        class notes for how to read the full chain or switch trimming
+        off.
 
         Exceptions propagate from the sampling engine: ``ValueError`` if
-        ``samples``, ``burn``, or ``thin`` are invalid, and ``RuntimeError``
-        if the active minimizer is not a BUMPS instance.
+        ``samples``, ``burn``, or ``thin`` are invalid, and
+        ``RuntimeError`` if the active minimizer is not a BUMPS
+        instance.
         """
         if self._state is not None:
             global_object.log.getLogger('fitting').warning(
@@ -618,38 +656,43 @@ def extend(
         progress_callback: Callable[[dict], bool | None] | None = None,
         abort_test: Callable[[], bool] | None = None,
     ) -> SamplingResults:
-        """Continue the existing chain with additional samples.
+        """
+        Continue the existing chain with additional samples.
 
         DREAM stores draws in a fixed-size ring buffer sized to its
-        ``samples`` parameter; this method does the ring-buffer arithmetic
-        for you (``samples = stored_generations * population +
-        additional_samples``) so no existing draws are dropped from the
-        buffer, regardless of the thinning interval. Runs with ``burn=0`` —
-        re-burning a converged chain is usually a mistake, and BUMPS forces
-        it to 0 on resume in any case. The DREAM population is recovered from
-        the saved state and cannot be changed on extend.
+        ``samples`` parameter; this method does the ring-buffer
+        arithmetic for you (``samples = stored_generations * population
+        + additional_samples``) so no existing draws are dropped from
+        the buffer, regardless of the thinning interval. Runs with
+        ``burn=0`` — re-burning a converged chain is usually a mistake,
+        and BUMPS forces it to 0 on resume in any case. The DREAM
+        population is recovered from the saved state and cannot be
+        changed on extend.
 
         Parameters
         ----------
         additional_samples : int, default=5000
-            Number of additional DREAM samples to draw, in the same units
-            as ``samples`` in ``sample()``. This grows the ring buffer by
-            exactly ``additional_samples``; how many of the new draws become
-            *visible* in ``results.draws`` depends on the thinning interval
-            and on BUMPS' automatic trimming (see Notes).
+            Number of additional DREAM samples to draw, in the same
+            units as ``samples`` in ``sample()``. This grows the ring
+            buffer by exactly ``additional_samples``; how many of the
+            new draws become *visible* in ``results.draws`` depends on
+            the thinning interval and on BUMPS' automatic trimming (see
+            Notes).
         thin : int, default=10
             Thinning interval for the retained draws.
         total_samples : int | None, default=None
-            Advanced: total retained samples requested from the ring buffer,
-            **overriding** the ``additional_samples`` arithmetic. With
-            ``total_samples=N``, only the last N draws are retained.
+            Advanced: total retained samples requested from the ring
+            buffer, **overriding** the ``additional_samples``
+            arithmetic. With ``total_samples=N``, only the last N draws
+            are retained.
         sampler_kwargs : dict | None, default=None
             Additional keyword arguments forwarded to the BUMPS DREAM
             sampler (merged over the instance defaults).
         progress_callback : Callable[[dict], bool | None] | None, default=None
             Optional callback invoked at each DREAM generation.
         abort_test : Callable[[], bool] | None, default=None
-            Optional callable that returns ``True`` to abort sampling early.
+            Optional callable that returns ``True`` to abort sampling
+            early.
 
         Returns
         -------
@@ -665,16 +708,17 @@ def extend(
         Notes
         -----
         ``results.draws`` will **not** grow by exactly
-        ``additional_samples / thin``. Two things get in the way, neither of
-        them a re-applied burn-in: BUMPS re-runs its burn-point detector over
-        the whole extended chain and re-trims the returned view (so the
-        visible count can even shrink), and DREAM advances in blocks of 10
-        generations, so the raw growth is ``additional_samples`` rounded up
-        to a multiple of ``10 * n_chains`` (i.e. at least
-        ``additional_samples / thin`` retained rows). To see the chain
-        itself, compare ``results.state.draw(portion=1.0, outliers=True)``
-        before and after, or run with ``sampler_kwargs={'trim': False}``. See
-        the ``Sampler`` class notes.
+        ``additional_samples / thin``. Two things get in the way,
+        neither of them a re-applied burn-in: BUMPS re-runs its
+        burn-point detector over the whole extended chain and re-trims
+        the returned view (so the visible count can even shrink), and
+        DREAM advances in blocks of 10 generations, so the raw growth is
+        ``additional_samples`` rounded up to a multiple of ``10 *
+        n_chains`` (i.e. at least ``additional_samples / thin`` retained
+        rows). To see the chain itself, compare
+        ``results.state.draw(portion=1.0, outliers=True)`` before and
+        after, or run with ``sampler_kwargs={'trim': False}``. See the
+        ``Sampler`` class notes.
         """
         if self._state is None:
             raise RuntimeError('No chain to extend. Call sample() or load_state() first.')
@@ -698,14 +742,15 @@ def extend(
         )
 
     def save(self, path: str | os.PathLike) -> None:
-        """Persist the chain state and metadata to disk.
+        """
+        Persist the chain state and metadata to disk.
 
         Writes the BUMPS native files (``<path>-chain.mc.gz``,
         ``<path>-point.mc.gz`` and ``<path>-stats.mc.gz``) plus a
         ``<path>.params.json`` sidecar with the parameter names, the
-        easyscience version, and a fingerprint of the bound data (verified
-        with a warning on ``load_state()``). Use ``load_chain`` to read the
-        files back without a fitter.
+        easyscience version, and a fingerprint of the bound data
+        (verified with a warning on ``load_state()``). Use
+        ``load_chain`` to read the files back without a fitter.
 
         Parameters
         ----------
@@ -746,13 +791,14 @@ def save(self, path: str | os.PathLike) -> None:
             json.dump(sidecar, f, indent=2)
 
     def load_state(self, path: str | os.PathLike, skip: int = 0) -> SamplingResults:
-        """Load a previously saved chain into this sampler.
+        """
+        Load a previously saved chain into this sampler.
 
-        The sampler must be constructed with the same fitter and data used to
-        create the chain — ``extend()`` then continues the saved chain. If the
-        sidecar carries a data fingerprint and it does not match this
-        sampler's bound data, a warning is logged (extending a chain against
-        different data is undefined behaviour).
+        The sampler must be constructed with the same fitter and data
+        used to create the chain — ``extend()`` then continues the saved
+        chain. If the sidecar carries a data fingerprint and it does not
+        match this sampler's bound data, a warning is logged (extending
+        a chain against different data is undefined behaviour).
 
         Populates ``state`` and ``results`` (draws, log-posterior and
         parameter names) from the saved chain, so summaries and
@@ -763,8 +809,8 @@ def load_state(self, path: str | os.PathLike, skip: int = 0) -> SamplingResults:
         path : str | os.PathLike
             File path prefix used in ``save()``.
         skip : int, default=0
-            Discard the first ``skip`` saved generations on load. Useful for
-            trimming additional burn-in without re-sampling.
+            Discard the first ``skip`` saved generations on load. Useful
+            for trimming additional burn-in without re-sampling.
 
         Returns
         -------
@@ -780,11 +826,11 @@ def load_state(self, path: str | os.PathLike, skip: int = 0) -> SamplingResults:
 
         Notes
         -----
-        BUMPS does not persist its automatic trim point, so a reloaded chain
-        comes back **untrimmed**: this method reports more draws than the
-        ``sample()`` call that wrote the file, even though the chain is
-        identical. Use ``skip`` to discard leading generations explicitly.
-        See the ``Sampler`` class notes.
+        BUMPS does not persist its automatic trim point, so a reloaded
+        chain comes back **untrimmed**: this method reports more draws
+        than the ``sample()`` call that wrote the file, even though the
+        chain is identical. Use ``skip`` to discard leading generations
+        explicitly. See the ``Sampler`` class notes.
         """  # noqa: DOC502 -- raised in _validate_chain_path via load_chain
         state, param_names, sidecar = load_chain(path, skip=skip)
 
diff --git a/src/easyscience/variable/descriptor_array.py b/src/easyscience/variable/descriptor_array.py
index 35f1c022..34a55393 100644
--- a/src/easyscience/variable/descriptor_array.py
+++ b/src/easyscience/variable/descriptor_array.py
@@ -23,6 +23,9 @@
 
 from .descriptor_base import DescriptorBase
 from .descriptor_number import DescriptorNumber
+from .units import has_numeric_factor
+from .units import normalisation_target
+from .units import si_base_unit
 
 
 class DescriptorArray(DescriptorBase):
@@ -97,6 +100,7 @@ def __init__(
         except Exception as message:
             raise UnitError(message)
             # TODO: handle 1xn and nx1 arrays
+        self._remember_unit(unit)
 
         super().__init__(
             name=name,
@@ -107,12 +111,22 @@ def __init__(
             parent=parent,
         )
 
-        # Call convert_unit during initialization to ensure that the unit has no numbers in it, and to ensure unit consistency.
-        if self.unit is not None:
-            self.convert_unit(self._base_unit())
+        # Make sure no magnitude is left hiding inside the unit, e.g. that a unit of
+        # 'm/mm' becomes a dimensionless value scaled by 1000. This must not be recorded
+        # on the undo stack: it is part of building the object.
+        target_unit = normalisation_target(
+            self._array.unit, has_spelling=self._input_unit is not None
+        )
+        if target_unit is not None:
+            self.convert_unit(target_unit, record_undo=False)
+            # The normalised unit is scipp's choice, not the user's, so there is no
+            # spelling to remember for it.
+            self._remember_unit(None)
 
     @classmethod
-    def from_scipp(cls, name: str, full_value: Variable, **kwargs: Any) -> DescriptorArray:
+    def from_scipp(
+        cls, name: str, full_value: Variable, sources: tuple = (), **kwargs: Any
+    ) -> DescriptorArray:
         """
         Create a DescriptorArray from a scipp array.
 
@@ -122,6 +136,10 @@ def from_scipp(cls, name: str, full_value: Variable, **kwargs: Any) -> Descripto
             Name of the descriptor.
         full_value : Variable
             Value of the descriptor as a scipp variable.
+        sources : tuple, default=()
+            Operands of the operation which produced ``full_value``, if
+            any. Used to display the result with the unit spelling its
+            operands were given.
         **kwargs : Any
             Additional parameters for the descriptor.
 
@@ -140,7 +158,7 @@ def from_scipp(cls, name: str, full_value: Variable, **kwargs: Any) -> Descripto
         return cls(
             name=name,
             value=full_value.values,
-            unit=full_value.unit,
+            unit=cls._spelling_from_sources(full_value.unit, sources),
             variance=full_value.variances,
             dimensions=full_value.dims,
             **kwargs,
@@ -266,11 +284,20 @@ def unit(self) -> str:
         """
         Get the unit.
 
+        The unit is reported with the spelling it was given, rather than
+        with scipp's preferred name for it, so that an array created
+        with 'angstrom' does not report 'Å'. The remembered spelling is
+        only used while it still describes the array we hold; if the
+        array has since been converted or normalised, scipp's own name
+        is reported instead.
+
         Returns
         -------
         str
             Unit as a string.
         """
+        if self._input_unit_parsed is not None and self._input_unit_parsed == self._array.unit:
+            return self._input_unit
         return str(self._array.unit)
 
     @unit.setter
@@ -387,7 +414,7 @@ def error(self, error: Union[list, np.ndarray]) -> None:
         else:
             self._array.variances = None
 
-    def convert_unit(self, unit_str: str) -> None:
+    def convert_unit(self, unit_str: str, record_undo: bool = True) -> None:
         """
         Convert the value from one unit system to another.
 
@@ -395,6 +422,10 @@ def convert_unit(self, unit_str: str) -> None:
         ----------
         unit_str : str
             New unit in string form.
+        record_undo : bool, default=True
+            Whether to push the conversion onto the undo stack. False
+            while constructing the object, where there is nothing to
+            undo back to.
 
         Raises
         ------
@@ -408,7 +439,7 @@ def convert_unit(self, unit_str: str) -> None:
         new_unit = sc.Unit(unit_str)
 
         # Save the current state for undo/redo
-        old_array = self._array
+        old_state = (self._array, self._input_unit, self._input_unit_parsed)
 
         # Perform the unit conversion
         try:
@@ -416,19 +447,86 @@ def convert_unit(self, unit_str: str) -> None:
         except Exception as e:
             raise UnitError(f'Failed to convert unit: {e}') from e
 
-        # Define the setter function for the undo stack
-        def set_array(obj, scalar):
-            obj._array = scalar
+        self._array = new_array
+        self._remember_unit(unit_str)
 
-        # Push to undo stack
-        self._global_object.stack.push(
-            PropertyStack(
-                self, set_array, old_array, new_array, text=f'Convert unit to {unit_str}'
+        # Define the setter function for the undo stack
+        def set_unit_state(obj, state):
+            obj._array, obj._input_unit, obj._input_unit_parsed = state
+
+        if record_undo:
+            self._global_object.stack.push(
+                PropertyStack(
+                    self,
+                    set_unit_state,
+                    old_state,
+                    (self._array, self._input_unit, self._input_unit_parsed),
+                    text=f'Convert unit to {unit_str}',
+                )
             )
-        )
 
-        # Update the array
-        self._array = new_array
+    @staticmethod
+    def _spelling_from_sources(unit: sc.Unit, sources: tuple) -> Union[str, sc.Unit]:
+        """
+        Return an operand's spelling for ``unit``, if one of them has
+        the same unit.
+
+        An operation such as an addition, or a multiplication by a plain
+        number, leaves the unit untouched, and the result should be
+        displayed the way its operands were rather than falling back to
+        scipp's name for it. Only an exactly equal unit is used, so a
+        result can never be relabelled as something it is not.
+
+        Parameters
+        ----------
+        unit : sc.Unit
+            Unit of the result.
+        sources : tuple
+            Operands of the operation. Anything which is not a
+            descriptor, such as a plain number, is ignored.
+
+        Returns
+        -------
+        Union[str, sc.Unit]
+            The operand's spelling, or ``unit`` unchanged if no operand
+            offers one.
+        """
+        for source in sources:
+            parsed_unit = getattr(source, '_input_unit_parsed', None)
+            if parsed_unit is not None and parsed_unit == unit:
+                return source._input_unit
+        return unit
+
+    def _remember_unit(self, unit: Union[str, sc.Unit, None]) -> None:
+        """
+        Remember the spelling a unit was given with, for display
+        purposes.
+
+        Nothing is remembered for a unit which did not arrive as a
+        string, or for a dimensionless one: reporting 'one' or ''
+        instead of 'dimensionless' would break the comparisons against
+        'dimensionless' made throughout this class.
+
+        Parameters
+        ----------
+        unit : Union[str, sc.Unit, None]
+            Unit as it was supplied.
+        """
+        self._input_unit = None
+        self._input_unit_parsed = None
+        if not isinstance(unit, str) or not unit.strip():
+            return
+        if has_numeric_factor(unit.strip()):
+            # A spelling such as '10dm^2' is no better than what scipp would print.
+            return
+        try:
+            parsed_unit = sc.Unit(unit.strip())
+        except Exception:
+            return
+        if parsed_unit == sc.units.dimensionless:
+            return
+        self._input_unit = unit.strip()
+        self._input_unit_parsed = parsed_unit
 
     def __copy__(self) -> DescriptorArray:
         """Return a copy of the current DescriptorArray."""
@@ -481,7 +579,7 @@ def as_dict(self, skip: Optional[List[str]] = None) -> Dict[str, Any]:
         """
         raw_dict = super().as_dict(skip=skip)
         raw_dict['value'] = self._array.values
-        raw_dict['unit'] = str(self._array.unit)
+        raw_dict['unit'] = self.unit
         raw_dict['variance'] = self._array.variances
         raw_dict['dimensions'] = self._array.dims
         return raw_dict
@@ -591,7 +689,9 @@ def _apply_operation(
         else:
             return NotImplemented
 
-        descriptor_array = DescriptorArray.from_scipp(name=self.name, full_value=new_full_value)
+        descriptor_array = DescriptorArray.from_scipp(
+            name=self.name, full_value=new_full_value, sources=(self, other)
+        )
         descriptor_array.name = descriptor_array.unique_name
         return descriptor_array
 
@@ -887,7 +987,9 @@ def __pow__(self, other: Union[DescriptorNumber, numbers.Number]) -> DescriptorA
             raise message from None
         if np.any(np.isnan(new_value.values)):
             raise ValueError('The result of the exponentiation is not a number')
-        descriptor_number = DescriptorArray.from_scipp(name=self.name, full_value=new_value)
+        descriptor_number = DescriptorArray.from_scipp(
+            name=self.name, full_value=new_value, sources=(self, other)
+        )
         descriptor_number.name = descriptor_number.unique_name
         return descriptor_number
 
@@ -903,7 +1005,9 @@ def __rpow__(self, other: numbers.Number):
     def __neg__(self) -> DescriptorArray:
         """Negate all values in the DescriptorArray."""
         new_value = -self.full_value
-        descriptor_array = DescriptorArray.from_scipp(name=self.name, full_value=new_value)
+        descriptor_array = DescriptorArray.from_scipp(
+            name=self.name, full_value=new_value, sources=(self,)
+        )
         descriptor_array.name = descriptor_array.unique_name
         return descriptor_array
 
@@ -916,7 +1020,9 @@ def __abs__(self) -> DescriptorArray:
         DescriptorArray.
         """
         new_value = abs(self.full_value)
-        descriptor_array = DescriptorArray.from_scipp(name=self.name, full_value=new_value)
+        descriptor_array = DescriptorArray.from_scipp(
+            name=self.name, full_value=new_value, sources=(self,)
+        )
         descriptor_array.name = descriptor_array.unique_name
         return descriptor_array
 
@@ -927,7 +1033,7 @@ def __getitem__(self, a) -> DescriptorArray:
         Defer slicing to scipp.
         """
         descriptor = DescriptorArray.from_scipp(
-            name=self.name, full_value=self.full_value.__getitem__(a)
+            name=self.name, full_value=self.full_value.__getitem__(a), sources=(self,)
         )
         descriptor.name = descriptor.unique_name
         return descriptor
@@ -1027,7 +1133,7 @@ def trace(
             )
             constructor = DescriptorArray.from_scipp
 
-        descriptor = constructor(name=self.name, full_value=trace)
+        descriptor = constructor(name=self.name, full_value=trace, sources=(self,))
         descriptor.name = descriptor.unique_name
         return descriptor
 
@@ -1057,7 +1163,7 @@ def sum(
         else:
             constructor = DescriptorArray.from_scipp
 
-        descriptor = constructor(name=self.name, full_value=new_full_value)
+        descriptor = constructor(name=self.name, full_value=new_full_value, sources=(self,))
         descriptor.name = descriptor.unique_name
         return descriptor
 
@@ -1084,18 +1190,3 @@ def sum(
     #
     #     other = sc.array(dims=self._array.dims, values=other)
     #     new_full_value = operation(self._array, other)  # Let scipp handle operation for uncertainty propagation
-
-    def _base_unit(self) -> str:
-        """
-        Returns the base unit of the current array.
-
-        For example, if the unit is ``100m``, returns ``m``.
-        """
-        string = str(self._array.unit)
-        for i, letter in enumerate(string):
-            if letter == 'e':
-                if string[i : i + 2] not in ['e+', 'e-']:
-                    return string[i:]
-            elif letter not in ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '.', '+', '-']:
-                return string[i:]
-        return ''
diff --git a/src/easyscience/variable/descriptor_number.py b/src/easyscience/variable/descriptor_number.py
index 5bcbcb0d..aa8810d3 100644
--- a/src/easyscience/variable/descriptor_number.py
+++ b/src/easyscience/variable/descriptor_number.py
@@ -21,6 +21,9 @@
 from easyscience.global_object.undo_redo import property_stack
 
 from .descriptor_base import DescriptorBase
+from .units import has_numeric_factor
+from .units import normalisation_target
+from .units import si_base_unit
 
 
 # Why is this a decorator? Because otherwise we would need a flag on the convert_unit method to avoid
@@ -48,6 +51,24 @@ def wrapper(self, *args, **kwargs):
     return wrapper
 
 
+def _set_unit_state(obj: Any, state: tuple) -> None:
+    """
+    Restore a unit state on the undo stack.
+
+    A module level function rather than a closure, so that the call
+    dispatches to the ``_restore_unit_state`` of whichever subclass owns
+    the object.
+
+    Parameters
+    ----------
+    obj : Any
+        Object to restore the state on.
+    state : tuple
+        State captured by ``_unit_state``.
+    """
+    obj._restore_unit_state(state)
+
+
 class DescriptorNumber(DescriptorBase):
     """
     A ``Descriptor`` for Number values with units.
@@ -101,6 +122,7 @@ def __init__(
             self._scalar = sc.scalar(float(value), unit=unit, variance=variance)
         except Exception as message:
             raise UnitError(message)
+        self._remember_unit(unit)
         super().__init__(
             name=name,
             unique_name=unique_name,
@@ -109,13 +131,23 @@ def __init__(
             display_name=display_name,
             parent=parent,
         )
-
-        # Call convert_unit during initialization to ensure that the unit has no numbers in it, and to ensure unit consistency.
-        if self.unit is not None:
-            self._convert_unit(self._base_unit())
+        # Make sure no magnitude is left hiding inside the unit. This has to happen after
+        # super().__init__, so that subclasses which convert more than the scalar (such
+        # as Parameter, with its bounds) are fully constructed. It should not be recorded
+        # on the undo stack as it is part of building the object.
+        target_unit = normalisation_target(
+            self._scalar.unit, has_spelling=self._input_unit is not None
+        )
+        if target_unit is not None:
+            self._convert_unit(target_unit, record_undo=False)
+            # The normalised unit is scipp's choice, not the user's, so there is no
+            # spelling to remember for it.
+            self._remember_unit(None)
 
     @classmethod
-    def from_scipp(cls, name: str, full_value: Variable, **kwargs: Any) -> DescriptorNumber:
+    def from_scipp(
+        cls, name: str, full_value: Variable, sources: tuple = (), **kwargs: Any
+    ) -> DescriptorNumber:
         """
         Create a DescriptorNumber from a scipp constant.
 
@@ -125,6 +157,10 @@ def from_scipp(cls, name: str, full_value: Variable, **kwargs: Any) -> Descripto
             Name of the descriptor.
         full_value : Variable
             Value of the descriptor as a scipp scalar.
+        sources : tuple, default=()
+            Operands of the operation which produced ``full_value``, if
+            any. Used to display the result with the unit spelling its
+            operands were given.
         **kwargs : Any
             Additional parameters for the descriptor.
 
@@ -145,7 +181,7 @@ def from_scipp(cls, name: str, full_value: Variable, **kwargs: Any) -> Descripto
         return cls(
             name=name,
             value=full_value.value,
-            unit=full_value.unit,
+            unit=cls._spelling_from_sources(full_value.unit, sources),
             variance=full_value.variance,
             **kwargs,
         )
@@ -258,11 +294,19 @@ def unit(self) -> str:
         """
         Get the unit.
 
+        The unit is reported with the spelling it was given, rather than
+        with scipp's preferred name for it. The remembered spelling is
+        only used while it still describes the scalar we hold; if the
+        scalar has since been converted or normalised, scipp's own name
+        is reported instead.
+
         Returns
         -------
         str
             Unit as a string.
         """
+        if self._input_unit_parsed is not None and self._input_unit_parsed == self._scalar.unit:
+            return self._input_unit
         return str(self._scalar.unit)
 
     @unit.setter
@@ -358,7 +402,7 @@ def error(self, value: float) -> None:
 
     # When we convert units internally, we dont want to notify observers as this can cause infinite recursion.
     # Therefore the convert_unit method is split into two methods, a private internal method and a public method.
-    def _convert_unit(self, unit_str: str) -> None:
+    def _convert_unit(self, unit_str: str, record_undo: bool = True) -> None:
         """
         Convert the value from one unit system to another.
 
@@ -366,6 +410,10 @@ def _convert_unit(self, unit_str: str) -> None:
         ----------
         unit_str : str
             New unit in string form.
+        record_undo : bool, default=True
+            Whether to push the conversion onto the undo stack. False
+            while constructing the object, where there is nothing to
+            undo back to.
 
         Raises
         ------
@@ -373,33 +421,142 @@ def _convert_unit(self, unit_str: str) -> None:
             If ``unit_str`` is not a string.
         UnitError
             If the unit conversion fails.
-        """
+        """  # noqa: DOC503. UnitError re-raised after restoring the unit state
         if not isinstance(unit_str, str):
             raise TypeError(f'{unit_str=} must be a string representing a valid scipp unit')
         new_unit = sc.Unit(unit_str)
 
         # Save the current state for undo/redo
-        old_scalar = self._scalar
+        old_state = self._unit_state()
 
         # Perform the unit conversion
         try:
-            new_scalar = self._scalar.to(unit=new_unit)
+            self._apply_unit_conversion(new_unit)
+        except Exception:
+            self._restore_unit_state(old_state)
+            raise
+        self._remember_unit(unit_str)
+
+        if record_undo:
+            self._global_object.stack.push(
+                PropertyStack(
+                    self,
+                    _set_unit_state,
+                    old_state,
+                    self._unit_state(),
+                    text=f'Convert unit to {unit_str}',
+                )
+            )
+
+    def _apply_unit_conversion(self, new_unit: sc.Unit) -> None:
+        """
+        Convert everything that carries this object's unit to
+        ``new_unit``.
+
+        Subclasses that hold more than the scalar in the unit, such as
+        ``Parameter`` with its bounds, extend this so that the whole
+        conversion is a single undoable step.
+
+        Parameters
+        ----------
+        new_unit : sc.Unit
+            Unit to convert to.
+
+        Raises
+        ------
+        UnitError
+            If the unit conversion fails.
+        """
+        try:
+            self._scalar = self._scalar.to(unit=new_unit)
         except Exception as e:
             raise UnitError(f'Failed to convert unit: {e}') from e
 
-        # Define the setter function for the undo stack
-        def set_scalar(obj, scalar):
-            obj._scalar = scalar
+    @staticmethod
+    def _spelling_from_sources(unit: sc.Unit, sources: tuple) -> Union[str, sc.Unit]:
+        """
+        Return an operand's spelling for ``unit``, if one of them has
+        the same unit.
 
-        # Push to undo stack
-        self._global_object.stack.push(
-            PropertyStack(
-                self, set_scalar, old_scalar, new_scalar, text=f'Convert unit to {unit_str}'
-            )
-        )
+        An operation such as an addition, or a multiplication by a plain
+        number, leaves the unit untouched, and the result should be
+        displayed the way its operands were rather than falling back to
+        scipp's name for it. Only an exactly equal unit is used, so a
+        result can never be relabelled as something it is not.
+
+        Parameters
+        ----------
+        unit : sc.Unit
+            Unit of the result.
+        sources : tuple
+            Operands of the operation. Anything which is not a
+            descriptor, such as a plain number, is ignored.
 
-        # Update the scalar
-        self._scalar = new_scalar
+        Returns
+        -------
+        Union[str, sc.Unit]
+            The operand's spelling, or ``unit`` unchanged if no operand
+            offers one.
+        """
+        for source in sources:
+            parsed_unit = getattr(source, '_input_unit_parsed', None)
+            if parsed_unit is not None and parsed_unit == unit:
+                return source._input_unit
+        return unit
+
+    def _unit_state(self) -> tuple:
+        """
+        Capture everything a unit conversion changes, so that it can be
+        undone as one.
+
+        Returns
+        -------
+        tuple
+            Opaque state, to be passed back to ``_restore_unit_state``.
+        """
+        return (self._scalar, self._input_unit, self._input_unit_parsed)
+
+    def _restore_unit_state(self, state: tuple) -> None:
+        """
+        Restore state captured by ``_unit_state``.
+
+        Parameters
+        ----------
+        state : tuple
+            State to restore.
+        """
+        self._scalar, self._input_unit, self._input_unit_parsed = state
+
+    def _remember_unit(self, unit: str | sc.Unit | None) -> None:
+        """
+        Remember the spelling a unit was given with, for display
+        purposes.
+
+        Nothing is remembered for a unit which did not arrive as a
+        string, or for a dimensionless one: reporting 'one' or ''
+        instead of 'dimensionless' would break the comparisons against
+        'dimensionless' made throughout this class.
+
+        Parameters
+        ----------
+        unit : str | sc.Unit | None
+            Unit as it was supplied.
+        """
+        self._input_unit = None
+        self._input_unit_parsed = None
+        if not isinstance(unit, str) or not unit.strip():
+            return
+        if has_numeric_factor(unit.strip()):
+            # A spelling such as '10dm^2' is no better than what scipp would print.
+            return
+        try:
+            parsed_unit = sc.Unit(unit.strip())
+        except Exception:
+            return
+        if parsed_unit == sc.units.dimensionless:
+            return
+        self._input_unit = unit.strip()
+        self._input_unit_parsed = parsed_unit
 
     # When the user calls convert_unit, we want to notify observers of the change to propagate the change.
     @notify_observers
@@ -434,7 +591,7 @@ def __repr__(self) -> str:
             string += f'{self._scalar.value:.4f}'
             if self.variance:
                 string += f' \u00b1 {self.error:.4f}'
-        obj_unit = self._scalar.unit
+        obj_unit = self.unit
         if obj_unit == 'dimensionless':
             obj_unit = ''
         else:
@@ -447,7 +604,7 @@ def __repr__(self) -> str:
     def as_dict(self, skip: Optional[List[str]] = None) -> Dict[str, Any]:
         raw_dict = super().as_dict(skip=skip)
         raw_dict['value'] = self._scalar.value
-        raw_dict['unit'] = str(self._scalar.unit)
+        raw_dict['unit'] = self.unit
         raw_dict['variance'] = self._scalar.variance
         if hasattr(self, '_DescriptorNumber__serializer_id'):
             raw_dict['__serializer_id'] = self.__serializer_id
@@ -470,7 +627,9 @@ def __add__(self, other: Union[DescriptorNumber, numbers.Number]) -> DescriptorN
             other._convert_unit(original_unit)
         else:
             return NotImplemented
-        descriptor_number = DescriptorNumber.from_scipp(name=self.name, full_value=new_value)
+        descriptor_number = DescriptorNumber.from_scipp(
+            name=self.name, full_value=new_value, sources=(self, other)
+        )
         descriptor_number.name = descriptor_number.unique_name
         return descriptor_number
 
@@ -481,7 +640,9 @@ def __radd__(self, other: numbers.Number) -> DescriptorNumber:
             new_value = other + self.full_value
         else:
             return NotImplemented
-        descriptor_number = DescriptorNumber.from_scipp(name=self.name, full_value=new_value)
+        descriptor_number = DescriptorNumber.from_scipp(
+            name=self.name, full_value=new_value, sources=(self, other)
+        )
         descriptor_number.name = descriptor_number.unique_name
         return descriptor_number
 
@@ -502,7 +663,9 @@ def __sub__(self, other: Union[DescriptorNumber, numbers.Number]) -> DescriptorN
             other._convert_unit(original_unit)
         else:
             return NotImplemented
-        descriptor_number = DescriptorNumber.from_scipp(name=self.name, full_value=new_value)
+        descriptor_number = DescriptorNumber.from_scipp(
+            name=self.name, full_value=new_value, sources=(self, other)
+        )
         descriptor_number.name = descriptor_number.unique_name
         return descriptor_number
 
@@ -513,7 +676,9 @@ def __rsub__(self, other: numbers.Number) -> DescriptorNumber:
             new_value = other - self.full_value
         else:
             return NotImplemented
-        descriptor = DescriptorNumber.from_scipp(name=self.name, full_value=new_value)
+        descriptor = DescriptorNumber.from_scipp(
+            name=self.name, full_value=new_value, sources=(self, other)
+        )
         descriptor.name = descriptor.unique_name
         return descriptor
 
@@ -524,8 +689,9 @@ def __mul__(self, other: Union[DescriptorNumber, numbers.Number]) -> DescriptorN
             new_value = self.full_value * other.full_value
         else:
             return NotImplemented
-        descriptor_number = DescriptorNumber.from_scipp(name=self.name, full_value=new_value)
-        descriptor_number._convert_unit(descriptor_number._base_unit())
+        descriptor_number = DescriptorNumber.from_scipp(
+            name=self.name, full_value=new_value, sources=(self, other)
+        )
         descriptor_number.name = descriptor_number.unique_name
         return descriptor_number
 
@@ -534,7 +700,9 @@ def __rmul__(self, other: numbers.Number) -> DescriptorNumber:
             new_value = other * self.full_value
         else:
             return NotImplemented
-        descriptor_number = DescriptorNumber.from_scipp(name=self.name, full_value=new_value)
+        descriptor_number = DescriptorNumber.from_scipp(
+            name=self.name, full_value=new_value, sources=(self, other)
+        )
         descriptor_number.name = descriptor_number.unique_name
         return descriptor_number
 
@@ -549,8 +717,9 @@ def __truediv__(self, other: Union[DescriptorNumber, numbers.Number]) -> Descrip
             new_value = self.full_value / other.full_value
         else:
             return NotImplemented
-        descriptor_number = DescriptorNumber.from_scipp(name=self.name, full_value=new_value)
-        descriptor_number._convert_unit(descriptor_number._base_unit())
+        descriptor_number = DescriptorNumber.from_scipp(
+            name=self.name, full_value=new_value, sources=(self, other)
+        )
         descriptor_number.name = descriptor_number.unique_name
         return descriptor_number
 
@@ -561,7 +730,9 @@ def __rtruediv__(self, other: numbers.Number) -> DescriptorNumber:
             new_value = other / self.full_value
         else:
             return NotImplemented
-        descriptor_number = DescriptorNumber.from_scipp(name=self.name, full_value=new_value)
+        descriptor_number = DescriptorNumber.from_scipp(
+            name=self.name, full_value=new_value, sources=(self, other)
+        )
         descriptor_number.name = descriptor_number.unique_name
         return descriptor_number
 
@@ -582,7 +753,9 @@ def __pow__(self, other: Union[DescriptorNumber, numbers.Number]) -> DescriptorN
             raise message from None
         if np.isnan(new_value.value):
             raise ValueError('The result of the exponentiation is not a number')
-        descriptor_number = DescriptorNumber.from_scipp(name=self.name, full_value=new_value)
+        descriptor_number = DescriptorNumber.from_scipp(
+            name=self.name, full_value=new_value, sources=(self, other)
+        )
         descriptor_number.name = descriptor_number.unique_name
         return descriptor_number
 
@@ -599,26 +772,16 @@ def __rpow__(self, other: numbers.Number) -> numbers.Number:
 
     def __neg__(self) -> DescriptorNumber:
         new_value = -self.full_value
-        descriptor_number = DescriptorNumber.from_scipp(name=self.name, full_value=new_value)
+        descriptor_number = DescriptorNumber.from_scipp(
+            name=self.name, full_value=new_value, sources=(self,)
+        )
         descriptor_number.name = descriptor_number.unique_name
         return descriptor_number
 
     def __abs__(self) -> DescriptorNumber:
         new_value = abs(self.full_value)
-        descriptor_number = DescriptorNumber.from_scipp(name=self.name, full_value=new_value)
+        descriptor_number = DescriptorNumber.from_scipp(
+            name=self.name, full_value=new_value, sources=(self,)
+        )
         descriptor_number.name = descriptor_number.unique_name
         return descriptor_number
-
-    def _base_unit(self) -> str:
-        """
-        Extract the base unit from the unit string by removing numeric
-        components and scientific notation.
-        """
-        string = str(self._scalar.unit)
-        for i, letter in enumerate(string):
-            if letter == 'e':
-                if string[i : i + 2] not in ['e+', 'e-']:
-                    return string[i:]
-            elif letter not in ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '.', '+', '-']:
-                return string[i:]
-        return ''
diff --git a/src/easyscience/variable/parameter.py b/src/easyscience/variable/parameter.py
index 3ec654f7..f22ef351 100644
--- a/src/easyscience/variable/parameter.py
+++ b/src/easyscience/variable/parameter.py
@@ -231,6 +231,9 @@ def _update(self) -> None:
             )  # noqa: E501
             self._min.unit = temporary_parameter.unit
             self._max.unit = temporary_parameter.unit
+            # The scalar's unit was relabelled directly above, so the remembered
+            # spelling has to follow it.
+            self._remember_unit(temporary_parameter.unit)
 
             if self._desired_unit is not None:
                 self._convert_unit(self._desired_unit)
@@ -662,22 +665,44 @@ def error(self, value: float) -> None:
                 'This is a dependent parameter, its error cannot be set directly.'
             )
 
-    def _convert_unit(self, unit_str: str) -> None:
+    def _apply_unit_conversion(self, new_unit: sc.Unit) -> None:
         """
-        Perform unit conversion.
-
-        The value, max and min can change on unit change.
+        Convert the value and the bounds to ``new_unit``.
 
         Parameters
         ----------
-        unit_str : str
-            New unit.
+        new_unit : sc.Unit
+            Unit to convert to.
         """
-        super()._convert_unit(unit_str=unit_str)
-        new_unit = sc.Unit(unit_str)  # unit_str is tested in super method
+        super()._apply_unit_conversion(new_unit)
         self._min = self._min.to(unit=new_unit)
         self._max = self._max.to(unit=new_unit)
 
+    def _unit_state(self) -> tuple:
+        """
+        Capture the bounds along with the base class' unit state, so
+        that undoing a unit conversion returns the value and the bounds
+        together.
+
+        Returns
+        -------
+        tuple
+            Opaque state, to be passed back to ``_restore_unit_state``.
+        """
+        return (super()._unit_state(), self._min, self._max)
+
+    def _restore_unit_state(self, state: tuple) -> None:
+        """
+        Restore state captured by ``_unit_state``.
+
+        Parameters
+        ----------
+        state : tuple
+            State to restore.
+        """
+        base_state, self._min, self._max = state
+        super()._restore_unit_state(base_state)
+
     @notify_observers
     def convert_unit(self, unit_str: str) -> None:
         """
@@ -1067,7 +1092,11 @@ def __add__(self, other: Union[DescriptorNumber, Parameter, numbers.Number]) ->
         else:
             return NotImplemented
         parameter = Parameter.from_scipp(
-            name=self.name, full_value=new_full_value, min=min_value, max=max_value
+            name=self.name,
+            full_value=new_full_value,
+            min=min_value,
+            max=max_value,
+            sources=(self, other),
         )
         parameter.name = parameter.unique_name
         return parameter
@@ -1096,7 +1125,11 @@ def __radd__(self, other: Union[DescriptorNumber, numbers.Number]) -> Parameter:
         else:
             return NotImplemented
         parameter = Parameter.from_scipp(
-            name=self.name, full_value=new_full_value, min=min_value, max=max_value
+            name=self.name,
+            full_value=new_full_value,
+            min=min_value,
+            max=max_value,
+            sources=(self, other),
         )
         parameter.name = parameter.unique_name
         return parameter
@@ -1129,7 +1162,11 @@ def __sub__(self, other: Union[DescriptorNumber, Parameter, numbers.Number]) ->
         else:
             return NotImplemented
         parameter = Parameter.from_scipp(
-            name=self.name, full_value=new_full_value, min=min_value, max=max_value
+            name=self.name,
+            full_value=new_full_value,
+            min=min_value,
+            max=max_value,
+            sources=(self, other),
         )
         parameter.name = parameter.unique_name
         return parameter
@@ -1158,7 +1195,11 @@ def __rsub__(self, other: Union[DescriptorNumber, numbers.Number]) -> Parameter:
         else:
             return NotImplemented
         parameter = Parameter.from_scipp(
-            name=self.name, full_value=new_full_value, min=min_value, max=max_value
+            name=self.name,
+            full_value=new_full_value,
+            min=min_value,
+            max=max_value,
+            sources=(self, other),
         )
         parameter.name = parameter.unique_name
         return parameter
@@ -1168,7 +1209,7 @@ def __mul__(self, other: Union[DescriptorNumber, Parameter, numbers.Number]) ->
             new_full_value = self.full_value * other
             if other == 0:
                 descriptor_number = DescriptorNumber.from_scipp(
-                    name=self.name, full_value=new_full_value
+                    name=self.name, full_value=new_full_value, sources=(self, other)
                 )
                 descriptor_number.name = descriptor_number.unique_name
                 return descriptor_number
@@ -1181,7 +1222,7 @@ def __mul__(self, other: Union[DescriptorNumber, Parameter, numbers.Number]) ->
                 other.value == 0 and type(other) is DescriptorNumber
             ):  # Only return DescriptorNumber if other is strictly 0, i.e. not a parameter  # noqa: E501
                 descriptor_number = DescriptorNumber.from_scipp(
-                    name=self.name, full_value=new_full_value
+                    name=self.name, full_value=new_full_value, sources=(self, other)
                 )
                 descriptor_number.name = descriptor_number.unique_name
                 return descriptor_number
@@ -1206,9 +1247,12 @@ def __mul__(self, other: Union[DescriptorNumber, Parameter, numbers.Number]) ->
         min_value = min(combinations)
         max_value = max(combinations)
         parameter = Parameter.from_scipp(
-            name=self.name, full_value=new_full_value, min=min_value, max=max_value
+            name=self.name,
+            full_value=new_full_value,
+            min=min_value,
+            max=max_value,
+            sources=(self, other),
         )
-        parameter._convert_unit(parameter._base_unit())
         parameter.name = parameter.unique_name
         return parameter
 
@@ -1217,7 +1261,7 @@ def __rmul__(self, other: Union[DescriptorNumber, numbers.Number]) -> Parameter:
             new_full_value = other * self.full_value
             if other == 0:
                 descriptor_number = DescriptorNumber.from_scipp(
-                    name=self.name, full_value=new_full_value
+                    name=self.name, full_value=new_full_value, sources=(self, other)
                 )
                 descriptor_number.name = descriptor_number.unique_name
                 return descriptor_number
@@ -1228,7 +1272,7 @@ def __rmul__(self, other: Union[DescriptorNumber, numbers.Number]) -> Parameter:
             new_full_value = other.full_value * self.full_value
             if other.value == 0:
                 descriptor_number = DescriptorNumber.from_scipp(
-                    name=self.name, full_value=new_full_value
+                    name=self.name, full_value=new_full_value, sources=(self, other)
                 )
                 descriptor_number.name = descriptor_number.unique_name
                 return descriptor_number
@@ -1238,9 +1282,12 @@ def __rmul__(self, other: Union[DescriptorNumber, numbers.Number]) -> Parameter:
         min_value = min(combinations)
         max_value = max(combinations)
         parameter = Parameter.from_scipp(
-            name=self.name, full_value=new_full_value, min=min_value, max=max_value
+            name=self.name,
+            full_value=new_full_value,
+            min=min_value,
+            max=max_value,
+            sources=(self, other),
         )
-        parameter._convert_unit(parameter._base_unit())
         parameter.name = parameter.unique_name
         return parameter
 
@@ -1288,9 +1335,12 @@ def __truediv__(self, other: Union[DescriptorNumber, Parameter, numbers.Number])
         min_value = min(combinations)
         max_value = max(combinations)
         parameter = Parameter.from_scipp(
-            name=self.name, full_value=new_full_value, min=min_value, max=max_value
+            name=self.name,
+            full_value=new_full_value,
+            min=min_value,
+            max=max_value,
+            sources=(self, other),
         )
-        parameter._convert_unit(parameter._base_unit())
         parameter.name = parameter.unique_name
         return parameter
 
@@ -1303,7 +1353,7 @@ def __rtruediv__(self, other: Union[DescriptorNumber, numbers.Number]) -> Parame
             other_value = other
             if other_value == 0:
                 descriptor_number = DescriptorNumber.from_scipp(
-                    name=self.name, full_value=new_full_value
+                    name=self.name, full_value=new_full_value, sources=(self, other)
                 )
                 descriptor_number.name = descriptor_number.unique_name
                 return descriptor_number
@@ -1314,7 +1364,7 @@ def __rtruediv__(self, other: Union[DescriptorNumber, numbers.Number]) -> Parame
             other_value = other.value
             if other_value == 0:
                 descriptor_number = DescriptorNumber.from_scipp(
-                    name=self.name, full_value=new_full_value
+                    name=self.name, full_value=new_full_value, sources=(self, other)
                 )
                 descriptor_number.name = descriptor_number.unique_name
                 return descriptor_number
@@ -1337,9 +1387,12 @@ def __rtruediv__(self, other: Union[DescriptorNumber, numbers.Number]) -> Parame
         min_value = min(combinations)
         max_value = max(combinations)
         parameter = Parameter.from_scipp(
-            name=self.name, full_value=new_full_value, min=min_value, max=max_value
+            name=self.name,
+            full_value=new_full_value,
+            min=min_value,
+            max=max_value,
+            sources=(self, other),
         )
-        parameter._convert_unit(parameter._base_unit())
         parameter.name = parameter.unique_name
         return parameter
 
@@ -1366,7 +1419,7 @@ def __pow__(self, other: Union[DescriptorNumber, numbers.Number]) -> Parameter:
             raise ValueError('The result of the exponentiation is not a number')
         if exponent == 0:
             descriptor_number = DescriptorNumber.from_scipp(
-                name=self.name, full_value=new_full_value
+                name=self.name, full_value=new_full_value, sources=(self, other)
             )
             descriptor_number.name = descriptor_number.unique_name
             return descriptor_number
@@ -1392,7 +1445,11 @@ def __pow__(self, other: Union[DescriptorNumber, numbers.Number]) -> Parameter:
         min_value = min(combinations)
         max_value = max(combinations)
         parameter = Parameter.from_scipp(
-            name=self.name, full_value=new_full_value, min=min_value, max=max_value
+            name=self.name,
+            full_value=new_full_value,
+            min=min_value,
+            max=max_value,
+            sources=(self, other),
         )
         parameter.name = parameter.unique_name
         return parameter
@@ -1402,7 +1459,11 @@ def __neg__(self) -> Parameter:
         min_value = -self.max
         max_value = -self.min
         parameter = Parameter.from_scipp(
-            name=self.name, full_value=new_full_value, min=min_value, max=max_value
+            name=self.name,
+            full_value=new_full_value,
+            min=min_value,
+            max=max_value,
+            sources=(self,),
         )
         parameter.name = parameter.unique_name
         return parameter
@@ -1415,7 +1476,11 @@ def __abs__(self) -> Parameter:
         min_value = min(combinations)
         max_value = max(combinations)
         parameter = Parameter.from_scipp(
-            name=self.name, full_value=new_full_value, min=min_value, max=max_value
+            name=self.name,
+            full_value=new_full_value,
+            min=min_value,
+            max=max_value,
+            sources=(self,),
         )
         parameter.name = parameter.unique_name
         return parameter
diff --git a/src/easyscience/variable/units.py b/src/easyscience/variable/units.py
new file mode 100644
index 00000000..7b2eac5d
--- /dev/null
+++ b/src/easyscience/variable/units.py
@@ -0,0 +1,147 @@
+# SPDX-FileCopyrightText: 2026 EasyScience contributors <https://github.com/easyscience>
+# SPDX-License-Identifier: BSD-3-Clause
+
+"""
+Helpers for how units are stored and displayed.
+
+Scipp can display a unit in two ways that surprise a user, and they need
+different remedies:
+
+1. It prints a **numeric factor** because it has no name for a scaled unit, so
+   ``m/mm`` becomes ``1000`` and ``dm*m`` becomes ``0.1m^2``. The magnitude is then
+   hiding inside the unit. This module fixes that, by folding the factor back into
+   the value.
+2. It prints **a different name than the one that was written**, so ``angstrom``
+   becomes ``Å`` and ``nm*m/s`` becomes ``nGy*s``. That is a display concern only, and
+   is handled by the descriptors, which remember the spelling a unit arrived as. See
+   ``DescriptorNumber.unit``.
+
+The descriptors' ``unit`` property is therefore a *display* string;
+``_scalar.unit`` / ``_array.unit`` remains the source of truth for
+meaning. The two never disagree about what the unit is, only about how
+it is spelled.
+"""
+
+from __future__ import annotations
+
+import re
+from typing import Optional
+from typing import Union
+
+import scipp as sc
+
+# Scipp writes a scaled unit it cannot name as a numeric factor followed by the
+# dimensions, in any of the forms '1000', '0.1m^2', '2.57e-44*J^2' and '3.9e+43/J^2'.
+
+# ^: start of string
+# \s*: optional whitespace
+# (?P<number>[0-9]+(?:\.[0-9]+)?(?:[eE][+-]?[0-9]+)?): a number, possibly with decimal point and exponent, captured as 'number'
+# \s*: optional whitespace
+# (?P<separator>[*/])?: an optional separator, either '*' or '/', captured as 'separator'
+# \s*: optional whitespace
+# (?P<rest>.*): the rest of the string, captured as 'rest'
+_NUMERIC_FACTOR = re.compile(
+    r'^\s*(?P<number>[0-9]+(?:\.[0-9]+)?(?:[eE][+-]?[0-9]+)?)\s*(?P<separator>[*/])?\s*(?P<rest>.*)$'
+)
+
+
+def has_numeric_factor(unit: Union[str, sc.Unit]) -> bool:
+    """
+    Report whether a unit carries a numeric factor in the way it is
+    written.
+
+    Applied to a ``sc.Unit`` this asks how scipp displays it; applied to
+    a string it asks how the user wrote it, which is how a spelling such
+    as '10dm^2' is rejected while 'dm*m' is accepted.
+
+    Parameters
+    ----------
+    unit : Union[str, sc.Unit]
+        Unit to inspect.
+
+    Returns
+    -------
+    bool
+        True if the unit is written with a leading multiplier.
+    """
+    match = _NUMERIC_FACTOR.match(str(unit))
+    if match is None:
+        return False
+    if match.group('rest') == '' and match.group('separator') is None:
+        # A pure number, such as '1000'.
+        return True
+    # The leading '1' of '1/meV' is a placeholder for the numerator, not a factor.
+    return float(match.group('number')) != 1.0
+
+
+def si_base_unit(unit: sc.Unit) -> sc.Unit:
+    """
+    Return the unit with the same dimensions as ``unit`` but no
+    multiplier.
+
+    The dimensions are read from ``sc.Unit.to_dict()``, which stores the
+    base units and powers, which is more robust than parsing the string.
+
+    Parameters
+    ----------
+    unit : sc.Unit
+        Unit to reduce.
+
+    Returns
+    -------
+    sc.Unit
+        The corresponding SI base unit, or dimensionless if the unit has
+        no dimensions.
+    """
+    powers = unit.to_dict().get('powers', {})
+    if not powers:
+        return sc.units.dimensionless
+    return sc.Unit('*'.join(f'{base}^{power}' for base, power in powers.items()))
+
+
+def normalisation_target(unit: sc.Unit, has_spelling: bool = False) -> Optional[str]:
+    """
+    Return the unit to convert to so that no magnitude is left inside
+    the unit.
+
+    A unit which reduces to a pure number, such as 'm/mm', is always
+    folded away.
+
+    Beyond that, only units with nothing better to display are
+    converted. If the user's own spelling is being kept
+    (``has_spelling``), then nothing numeric is on show and there is
+    nothing to fix, so a parameter created with '1/meV^2' keeps both its
+    spelling and its value, even though scipp would print it as
+    '3.9e+43/J^2'.
+
+    Parameters
+    ----------
+    unit : sc.Unit
+        Unit to inspect.
+    has_spelling : bool, default=False
+        Whether a clean spelling for this unit is being displayed in its
+        place.
+
+    Returns
+    -------
+    Optional[str]
+        The unit to convert to, or None if there is nothing to fix, or
+        if no safe target could be determined.
+    """
+    if unit.to_dict().get('e_flag'):
+        # Offset units such as degC cannot be converted this way.
+        return None
+    if not unit.to_dict().get('powers'):
+        # A pure number, such as the 1000 that 'm/mm' reduces to.
+        return None if unit == sc.units.dimensionless else 'dimensionless'
+    if has_spelling or not has_numeric_factor(unit):
+        return None
+    try:
+        base_unit = si_base_unit(unit)
+        target = str(base_unit)
+        if sc.Unit(target) != base_unit:
+            # The target does not survive a round trip, so it is not safe to use.
+            return None
+    except Exception:
+        return None
+    return target
diff --git a/tests/unit/global_object/test_integration_comprehensive.py b/tests/unit/global_object/test_integration_comprehensive.py
index 61268954..a9ac4a17 100644
--- a/tests/unit/global_object/test_integration_comprehensive.py
+++ b/tests/unit/global_object/test_integration_comprehensive.py
@@ -422,6 +422,7 @@ def test_concurrent_access_simulation(self, clear_all):
         # Given
         global_obj = GlobalObject()
         results = []
+        created = []  # The map holds only weak references, so keep the objects alive
         errors = []
 
         def create_objects(thread_id, count=10):
@@ -429,6 +430,7 @@ def create_objects(thread_id, count=10):
             try:
                 for i in range(count):
                     param = Parameter(name=f'thread_{thread_id}_param_{i}', value=float(i))
+                    created.append(param)
                     results.append(param.unique_name)
                     time.sleep(0.001)  # Small delay to encourage race conditions
             except Exception as e:
diff --git a/tests/unit/variable/test_descriptor_array.py b/tests/unit/variable/test_descriptor_array.py
index 73bb0857..a7b7b4c2 100644
--- a/tests/unit/variable/test_descriptor_array.py
+++ b/tests/unit/variable/test_descriptor_array.py
@@ -226,19 +226,21 @@ def test_copy(self, descriptor: DescriptorArray):
         assert descriptor_copy._array.unit == descriptor._array.unit
 
     @pytest.mark.parametrize(
-        'unit_string, expected',
-        [('1e+9', 'dimensionless'), ('1000', 'dimensionless'), ('10dm^2', 'm^2')],
+        'unit_string, expected_unit, factor',
+        [
+            ('1e+9', 'dimensionless', 1e9),
+            ('1000', 'dimensionless', 1000.0),
+            ('10dm^2', 'm^2', 0.1),
+        ],
         ids=['scientific_notation', 'numbers', 'unit_prefix'],
     )
-    def test_base_unit(self, unit_string, expected):
+    def test_numeric_factor_is_folded_into_value(self, unit_string, expected_unit, factor):
         # When
         descriptor = DescriptorArray(name='name', value=[[1.0, 2.0], [3.0, 4.0]], unit=unit_string)
 
-        # Then
-        base_unit = descriptor._base_unit()
-
-        # Expect
-        assert base_unit == expected
+        # Expect: the magnitude ends up in the values, never inside the unit
+        assert descriptor.unit == expected_unit
+        assert descriptor.value == pytest.approx(np.array([[1.0, 2.0], [3.0, 4.0]]) * factor)
 
     @pytest.mark.parametrize(
         'test, expected, raises_warning',
@@ -662,7 +664,7 @@ def test_reverse_subtraction_dimensionless(
             (
                 DescriptorNumber('test', 1, 'kg', 10),
                 DescriptorArray(
-                    'test * name', [[1.0, 2.0], [3.0, 4.0]], 'kg*m', [[10.1, 40.2], [90.3, 160.4]]
+                    'test * name', [[1.0, 2.0], [3.0, 4.0]], 'm*kg', [[10.1, 40.2], [90.3, 160.4]]
                 ),
                 True,
             ),
@@ -795,7 +797,7 @@ def test_multiplication_dimensionless(
             (
                 DescriptorNumber('test', 1, 'kg', 10),
                 DescriptorArray(
-                    'test * name', [[1.0, 2.0], [3.0, 4.0]], 'kg*m', [[10.1, 40.2], [90.3, 160.4]]
+                    'test * name', [[1.0, 2.0], [3.0, 4.0]], 'm*kg', [[10.1, 40.2], [90.3, 160.4]]
                 ),
                 True,
             ),
diff --git a/tests/unit/variable/test_descriptor_number.py b/tests/unit/variable/test_descriptor_number.py
index 6abd3029..9bb64d6f 100644
--- a/tests/unit/variable/test_descriptor_number.py
+++ b/tests/unit/variable/test_descriptor_number.py
@@ -213,19 +213,21 @@ def test_copy(self, descriptor: DescriptorNumber):
         assert descriptor_copy._scalar.unit == descriptor._scalar.unit
 
     @pytest.mark.parametrize(
-        'unit_string, expected',
-        [('1e+9', 'dimensionless'), ('1000', 'dimensionless'), ('10dm^2', 'm^2')],
+        'unit_string, expected_unit, expected_value',
+        [
+            ('1e+9', 'dimensionless', 1e9),
+            ('1000', 'dimensionless', 1000.0),
+            ('10dm^2', 'm^2', 0.1),
+        ],
         ids=['scientific_notation', 'numbers', 'unit_prefix'],
     )
-    def test_base_unit(self, unit_string, expected):
+    def test_numeric_factor_is_folded_into_value(self, unit_string, expected_unit, expected_value):
         # When
         descriptor = DescriptorNumber(name='name', value=1, unit=unit_string)
 
-        # Then
-        base_unit = descriptor._base_unit()
-
-        # Expect
-        assert base_unit == expected
+        # Expect: the magnitude ends up in the value, never inside the unit
+        assert descriptor.unit == expected_unit
+        assert descriptor.value == pytest.approx(expected_value)
 
     @pytest.mark.parametrize(
         'test, expected',
diff --git a/tests/unit/variable/test_parameter.py b/tests/unit/variable/test_parameter.py
index 1ac05b45..95d6ee8d 100644
--- a/tests/unit/variable/test_parameter.py
+++ b/tests/unit/variable/test_parameter.py
@@ -1678,7 +1678,7 @@ def test_division_with_parameter_remaining_cases(self, first, second, expected):
             (
                 2,
                 Parameter('name / 2', 0.5, 'm', 0.0025, 0, 5),
-                Parameter('2 / name', 2, 'm**-1', 0.04, 0.2, np.inf),
+                Parameter('2 / name', 2, '1/m', 0.04, 0.2, np.inf),
             ),
         ],
         ids=['descriptor_number', 'number'],
diff --git a/tests/unit/variable/test_units.py b/tests/unit/variable/test_units.py
new file mode 100644
index 00000000..60bf754d
--- /dev/null
+++ b/tests/unit/variable/test_units.py
@@ -0,0 +1,252 @@
+# SPDX-FileCopyrightText: 2026 EasyScience contributors <https://github.com/easyscience>
+# SPDX-License-Identifier: BSD-3-Clause
+
+import pytest
+import scipp as sc
+from scipp import UnitError
+
+from easyscience import DescriptorNumber
+from easyscience import Parameter
+from easyscience import global_object
+from easyscience.variable.units import has_numeric_factor
+from easyscience.variable.units import normalisation_target
+from easyscience.variable.units import si_base_unit
+
+
+class TestUnitHelpers:
+    @pytest.mark.parametrize(
+        'unit_string, expected',
+        [
+            ('m', False),
+            ('mm', False),
+            ('meV', False),
+            ('m/s', False),
+            ('counts', False),
+            ('1/m', False),
+            ('1/meV', False),
+            ('1/J**2', False),
+            ('mm*mm', False),
+            ('nm*m/s', False),
+            ('degC', False),
+            ('deg', False),
+            ('1000', True),
+            ('1e+9', True),
+            ('m/mm', True),
+            ('dm*m', True),
+            ('10dm^2', True),
+            ('meV*meV', True),
+            ('1/meV**2', True),
+        ],
+    )
+    def test_has_numeric_factor(self, unit_string, expected):
+        assert has_numeric_factor(sc.Unit(unit_string)) is expected
+
+    @pytest.mark.parametrize(
+        'unit_string, expected',
+        [
+            ('m', 'm'),
+            ('mm', 'm'),
+            ('dm*m', 'm^2'),
+            ('m/mm', 'dimensionless'),
+            ('1000', 'dimensionless'),
+            ('nm*m/s', 'm^2/s'),
+            ('meV', 'J'),
+            ('counts/s', 'counts/s'),
+        ],
+    )
+    def test_si_base_unit(self, unit_string, expected):
+        assert str(si_base_unit(sc.Unit(unit_string))) == expected
+
+    @pytest.mark.parametrize('unit_string', ['degC', 'deg', 'm', 'mm', 'counts'])
+    def test_normalisation_target_leaves_clean_units_alone(self, unit_string):
+        # Offset units such as degC must not be touched either: they cannot be
+        # converted to their base unit this way.
+        assert normalisation_target(sc.Unit(unit_string)) is None
+
+
+class TestNumericFactorFolding:
+    @pytest.mark.parametrize(
+        'unit_string, expected_unit, expected_value',
+        [
+            ('1e+9', 'dimensionless', 1e9),
+            ('1000', 'dimensionless', 1000.0),
+            ('m/mm', 'dimensionless', 1000.0),
+            ('10dm^2', 'm^2', 0.1),
+        ],
+        ids=['scientific_notation', 'number', 'cancelling_prefixes', 'unit_prefix'],
+    )
+    def test_construction_folds_the_factor(self, unit_string, expected_unit, expected_value):
+        descriptor = DescriptorNumber(name='name', value=1.0, unit=unit_string)
+
+        assert descriptor.unit == expected_unit
+        assert descriptor.value == pytest.approx(expected_value)
+
+    @pytest.mark.parametrize(
+        'left, right, expected_unit, expected_value',
+        [
+            ('m', 'mm', 'dimensionless', 1000.0),
+            ('m', 'm', 'dimensionless', 1.0),
+            ('m', 's', 'm/s', 1.0),
+        ],
+        ids=['cancelling_prefixes', 'identical', 'unrelated'],
+    )
+    def test_division_folds_the_factor(self, left, right, expected_unit, expected_value):
+        result = DescriptorNumber('a', 1.0, left) / DescriptorNumber('b', 1.0, right)
+
+        assert result.unit == expected_unit
+        assert result.value == pytest.approx(expected_value)
+
+    @pytest.mark.parametrize(
+        'left, right, expected_unit, expected_value',
+        [
+            ('dm', 'm', 'm^2', 0.2),
+            ('mm', 'mm', 'mm^2', 2.0),
+            ('cm', 'm', 'dm^2', 2.0),
+        ],
+        ids=['scaled', 'already_named', 'already_named_prefix'],
+    )
+    def test_multiplication_folds_only_when_needed(
+        self, left, right, expected_unit, expected_value
+    ):
+        # A unit scipp can name is left exactly as it names it; only a unit displayed
+        # with a numeric factor is converted.
+        result = DescriptorNumber('a', 2.0, left) * DescriptorNumber('b', 1.0, right)
+
+        assert result.unit == expected_unit
+        assert result.value == pytest.approx(expected_value)
+
+    def test_bounds_are_folded_with_the_value(self):
+        parameter = Parameter(name='name', value=1.0, min=0.0, max=10.0, unit='m/mm')
+
+        assert parameter.unit == 'dimensionless'
+        assert parameter.value == pytest.approx(1000.0)
+        assert parameter.min == pytest.approx(0.0)
+        assert parameter.max == pytest.approx(10000.0)
+
+    def test_invalid_unit_still_raises(self):
+        with pytest.raises(UnitError):
+            DescriptorNumber(name='name', value=1.0, unit='not_a_unit')
+
+
+class TestUnitSpelling:
+    @pytest.mark.parametrize(
+        'unit_string',
+        ['angstrom', 'nm*m/s', '1/meV^2', 'dm*m', 'Hz'],
+    )
+    def test_unit_is_reported_as_written(self, unit_string):
+        descriptor = DescriptorNumber(name='name', value=1.0, unit=unit_string)
+
+        assert descriptor.unit == unit_string
+        # The spelling is a display concern only; the meaning is unchanged.
+        assert descriptor._scalar.unit == sc.Unit(unit_string)
+
+    def test_spellings_are_per_object(self):
+        # The property a global alias table could not provide: scipp permits only one
+        # alias per unit, so one of these two would have had to report the other's
+        # spelling.
+        angstrom = DescriptorNumber(name='a', value=1.0, unit='angstrom')
+        symbol = DescriptorNumber(name='b', value=1.0, unit='Å')
+
+        assert angstrom.unit == 'angstrom'
+        assert symbol.unit == 'Å'
+
+    @pytest.mark.parametrize('unit_string', ['', 'dimensionless', 'one'])
+    def test_dimensionless_is_always_reported_as_dimensionless(self, unit_string):
+        # Reporting '' or 'one' would break the comparisons against 'dimensionless'
+        # made throughout DescriptorNumber.
+        descriptor = DescriptorNumber(name='name', value=1.0, unit=unit_string)
+
+        assert descriptor.unit == 'dimensionless'
+        assert (descriptor + 1.0).value == pytest.approx(2.0)
+
+    def test_derived_units_use_scipps_spelling(self):
+        result = DescriptorNumber('a', 1.0, 'm') * DescriptorNumber('b', 1.0, 'm')
+
+        assert result.unit == 'm^2'
+
+    @pytest.mark.parametrize('unit_string', ['angstrom', '1/meV^2', 'dm*m'])
+    def test_unit_preserving_operations_keep_the_spelling(self, unit_string):
+        # An operation which leaves the unit alone should leave its spelling alone too,
+        # and must not disturb the value on the way.
+        a = DescriptorNumber('a', 2.0, unit_string)
+        b = DescriptorNumber('b', 3.0, unit_string)
+
+        assert (a + b).unit == unit_string
+        assert (a + b).value == pytest.approx(5.0, rel=0, abs=0)
+        assert (a * 2).unit == unit_string
+        assert (a * 2).value == pytest.approx(4.0, rel=0, abs=0)
+        assert (-a).unit == unit_string
+        assert abs(a).unit == unit_string
+
+    def test_a_new_unit_is_never_relabelled(self):
+        # Only an exactly equal unit is inherited from an operand.
+        a = DescriptorNumber('a', 2.0, 'dm')
+        b = DescriptorNumber('b', 1.0, 'm')
+
+        assert (a * b).unit == 'm^2'
+        assert (a * b).value == pytest.approx(0.2)
+
+    def test_bounds_follow_an_inherited_spelling(self):
+        parameter = Parameter(name='p', value=1.0, min=0.0, max=10.0, unit='angstrom')
+
+        result = parameter * 2
+
+        assert result.unit == 'angstrom'
+        assert result.max == pytest.approx(20.0)
+
+    def test_spelling_falls_back_when_the_scalar_changes_underneath(self):
+        descriptor = DescriptorNumber(name='name', value=1.0, unit='angstrom')
+        assert descriptor.unit == 'angstrom'
+
+        # Bypass convert_unit entirely, as Parameter._update does.
+        descriptor._scalar = descriptor._scalar.to(unit='m')
+
+        assert descriptor.unit == 'm'
+
+    def test_convert_unit_adopts_the_new_spelling(self):
+        descriptor = DescriptorNumber(name='name', value=1.0, unit='m')
+
+        descriptor.convert_unit('angstrom')
+
+        assert descriptor.unit == 'angstrom'
+        assert descriptor.value == pytest.approx(1e10)
+
+    def test_serialisation_round_trip_preserves_the_spelling(self):
+        descriptor = DescriptorNumber(name='name', value=1.0, unit='nm*m/s')
+
+        restored = DescriptorNumber.from_dict(descriptor.as_dict(skip=['unique_name']))
+
+        assert restored.unit == 'nm*m/s'
+        assert restored._scalar.unit == descriptor._scalar.unit
+
+
+class TestUnitUndoRedo:
+    def test_undo_restores_value_bounds_and_spelling_together(self):
+        global_object.stack.enabled = True
+        try:
+            parameter = Parameter(name='name', value=1.0, min=0.0, max=10.0, unit='m')
+            parameter.convert_unit('mm')
+            assert parameter.unit == 'mm'
+            assert parameter.value == pytest.approx(1000.0)
+            assert parameter.max == pytest.approx(10000.0)
+
+            global_object.stack.undo()
+
+            # One undo, and everything the conversion touched comes back with it.
+            assert parameter.unit == 'm'
+            assert parameter.value == pytest.approx(1.0)
+            assert parameter.min == pytest.approx(0.0)
+            assert parameter.max == pytest.approx(10.0)
+        finally:
+            global_object.stack.enabled = False
+
+    def test_construction_does_not_record_an_undo_entry(self):
+        global_object.stack.enabled = True
+        try:
+            global_object.stack.clear()
+            # A unit needing normalisation must not leave anything to undo.
+            DescriptorNumber(name='name', value=1.0, unit='m/mm')
+
+            assert global_object.stack.canUndo() is False
+        finally:
+            global_object.stack.enabled = False