diff --git a/phaser/engines/common/noise_models.py b/phaser/engines/common/noise_models.py
index c366f24..d9eb53a 100644
--- a/phaser/engines/common/noise_models.py
+++ b/phaser/engines/common/noise_models.py
@@ -37,7 +37,7 @@ def calc_loss(
         xp = get_array_module(model_wave, model_intensity, exp_patterns, mask)
         patterns = xp.maximum(exp_patterns, 0.0)
 
-        return ((
+        return (t.cast(numpy.floating,
             2. * xp.sum(mask * (
                 xp.sqrt(patterns + self.offset) - xp.sqrt(model_intensity + self.offset) - self.eps
             )**2) / self.var).astype(exp_patterns.dtype),
diff --git a/phaser/engines/common/output.py b/phaser/engines/common/output.py
index 7100d18..31f8b7c 100644
--- a/phaser/engines/common/output.py
+++ b/phaser/engines/common/output.py
@@ -186,7 +186,7 @@ def _plot_scan(state: ReconsState, out_path: Path, options: SaveOptions):
     ax.set_xlim(left, right)
     ax.set_ylim(bottom, top)
 
-    scan = to_numpy(state.scan)
+    scan = to_numpy(state.scan.data)
     i = numpy.arange(scan[..., 0].size)
     ax.scatter(scan[..., 1].ravel(), scan[..., 0].ravel(), c=i, cmap='plasma', s=0.5, edgecolors='none')
 
@@ -197,7 +197,7 @@ def _plot_scan(state: ReconsState, out_path: Path, options: SaveOptions):
 def _plot_tilt(state: ReconsState, out_path: Path, options: SaveOptions):
     from matplotlib import pyplot
 
-    if state.tilt is None:
+    if state.scan.tilt is None:
         logger = logging.getLogger(__name__)
         logger.warning("Tilt map (`state.tilt`) is missing, skipping `plot_tilt`")
         return
@@ -209,8 +209,8 @@ def _plot_tilt(state: ReconsState, out_path: Path, options: SaveOptions):
     ax.set_xlim(left, right)
     ax.set_ylim(bottom, top)
 
-    scan = to_numpy(state.scan)
-    tilt = to_numpy(state.tilt)
+    scan = to_numpy(state.scan.data)
+    tilt = to_numpy(state.scan.tilt)
     tilt = tilt[..., 1] + tilt[..., 0]*1.j
     max_tilt = max(numpy.max(numpy.abs(tilt)), 1.0)  # at least 1 mrad
     c = colorize_complex(tilt.ravel() / max_tilt, amp=True, rescale=False)
diff --git a/phaser/engines/common/position_correction.py b/phaser/engines/common/position_correction.py
index 2f0809f..724540a 100644
--- a/phaser/engines/common/position_correction.py
+++ b/phaser/engines/common/position_correction.py
@@ -41,8 +41,8 @@ def __init__(self, args: None, props: MomentumPositionSolverProps):
         self.momentum = props.momentum
 
     def init_state(self, sim: ReconsState) -> NDArray[numpy.floating]:
-        xp = get_array_module(sim.scan)
-        return xp.zeros_like(sim.scan)
+        xp = get_array_module(sim.scan.data)
+        return xp.zeros_like(sim.scan.data)
 
     def perform_update(
         self,
diff --git a/phaser/engines/common/regularizers.py b/phaser/engines/common/regularizers.py
index d8a94f2..97b91ea 100644
--- a/phaser/engines/common/regularizers.py
+++ b/phaser/engines/common/regularizers.py
@@ -250,7 +250,7 @@ def calc_loss_group(
         xp = get_array_module(sim.object.data)
 
         cost = xp.sum(xp.abs(sim.object.data - 1.0))
-        cost_scale = xp.array(group.shape[-1] / prod(sim.scan.shape[:-1]), dtype=cost.dtype)
+        cost_scale = xp.array(group.shape[-1] / prod(sim.scan.data.shape[:-1]), dtype=cost.dtype)
         return (cost * cost_scale * self.cost, state)
 
 
@@ -272,7 +272,7 @@ def calc_loss_group(
 
         cost = xp.sum(abs2(sim.object.data - 1.0))
 
-        cost_scale = xp.array(group.shape[-1] / prod(sim.scan.shape[:-1]), dtype=cost.dtype)
+        cost_scale = xp.array(group.shape[-1] / prod(sim.scan.data.shape[:-1]), dtype=cost.dtype)
         return (cost * cost_scale * self.cost, state)  # type: ignore
 
 
@@ -293,7 +293,7 @@ def calc_loss_group(
         xp = get_array_module(sim.object.data)
 
         cost = xp.sum(xp.abs(xp.angle(sim.object.data)))
-        cost_scale = xp.array(group.shape[-1] / prod(sim.scan.shape[:-1]), dtype=cost.dtype)
+        cost_scale = xp.array(group.shape[-1] / prod(sim.scan.data.shape[:-1]), dtype=cost.dtype)
         return (cost * cost_scale * self.cost, state)
 
 
@@ -319,7 +319,7 @@ def calc_loss_group(
             xp.abs(fft2(xp.prod(sim.object.data, axis=0)))
         )
         # scale cost by fraction of the total reconstruction in the group
-        cost_scale = xp.array(group.shape[-1] / prod(sim.scan.shape[:-1]), dtype=cost.dtype)
+        cost_scale = xp.array(group.shape[-1] / prod(sim.scan.data.shape[:-1]), dtype=cost.dtype)
 
         return (cost * cost_scale * self.cost, state)
 
@@ -351,7 +351,7 @@ def calc_loss_group(
         #)
         # scale cost by fraction of the total reconstruction in the group
         # TODO also scale by # of pixels or similar?
-        cost_scale = xp.array(group.shape[-1] / prod(sim.scan.shape[:-1]), dtype=cost.dtype)
+        cost_scale = xp.array(group.shape[-1] / prod(sim.scan.data.shape[:-1]), dtype=cost.dtype)
 
         return (cost * cost_scale * self.cost, state)
 
@@ -377,7 +377,7 @@ def calc_loss_group(
             xp.sum(abs2(xp.diff(sim.object.data, axis=-2)))
         )
         # scale cost by fraction of the total reconstruction in the group
-        cost_scale = xp.array(group.shape[-1] / prod(sim.scan.shape[:-1]), dtype=cost.dtype)
+        cost_scale = xp.array(group.shape[-1] / prod(sim.scan.data.shape[:-1]), dtype=cost.dtype)
 
         return (cost * cost_scale * self.cost, state)  # type: ignore
 
@@ -403,7 +403,7 @@ def calc_loss_group(
 
         cost = xp.sum(xp.abs(xp.diff(sim.object.data, axis=0)))
         # scale cost by fraction of the total reconstruction in the group
-        cost_scale = xp.array(group.shape[-1] / prod(sim.scan.shape[:-1]), dtype=cost.dtype)
+        cost_scale = xp.array(group.shape[-1] / prod(sim.scan.data.shape[:-1]), dtype=cost.dtype)
 
         return (cost * cost_scale * self.cost, state)
 
@@ -429,7 +429,7 @@ def calc_loss_group(
 
         cost = xp.sum(abs2(xp.diff(sim.object.data, axis=0)))
         # scale cost by fraction of the total reconstruction in the group
-        cost_scale = xp.array(group.shape[-1] / prod(sim.scan.shape[:-1]), dtype=cost.dtype)
+        cost_scale = xp.array(group.shape[-1] / prod(sim.scan.data.shape[:-1]), dtype=cost.dtype)
 
         return (cost * cost_scale * self.cost, state)  # type: ignore
 
@@ -519,7 +519,7 @@ def __init__(self, args: None, props: UnstructuredGaussianProps):
         self.attr_path = props.attr_path
 
     def init_state(self, sim: ReconsState) -> NDArray[numpy.floating]:
-        xp = get_array_module(sim.scan)
+        xp = get_array_module(sim.scan.data)
         try:
             self.getattr_nested(sim, self.attr_path)
         except AttributeError as e:
@@ -547,8 +547,8 @@ def setattr_nested(self, obj: t.Any, attr_path: str, value: t.Any):
     def apply_iter(self, sim: ReconsState, state: NDArray[numpy.floating]) -> t.Tuple[ReconsState, NDArray[numpy.floating]]:
         from scipy.spatial import KDTree
         obj_samp = sim.object.sampling
-        scan_flat = sim.scan.reshape(-1, 2)
-        scan_ndim = sim.scan.ndim - 1
+        scan_flat = sim.scan.data.reshape(-1, 2)
+        scan_ndim = sim.scan.data.ndim - 1
 
         attr = self.getattr_nested(sim, self.attr_path)
         vals = t.cast(NDArray[numpy.inexact], getattr(attr, 'data', attr))  # Extract raw array
diff --git a/phaser/engines/common/simulation.py b/phaser/engines/common/simulation.py
index b6c6778..189ca06 100644
--- a/phaser/engines/common/simulation.py
+++ b/phaser/engines/common/simulation.py
@@ -249,7 +249,7 @@ def cutout_group(
     """Returns (probe, obj) in the cutout region"""
     probes = state.probe.data
 
-    group_scan = state.scan[tuple(group)]
+    group_scan = state.scan.data[tuple(group)]
     group_obj = state.object.sampling.get_view_at_pos(state.object.data, group_scan, probes.shape[-2:])
     # group probes in real space
     # shape (len(group), 1, Ny, Nx)
diff --git a/phaser/engines/conventional/run.py b/phaser/engines/conventional/run.py
index 7e5fc98..704d0a4 100644
--- a/phaser/engines/conventional/run.py
+++ b/phaser/engines/conventional/run.py
@@ -51,7 +51,7 @@ def run_engine(args: EngineArgs, props: ConventionalEnginePlan) -> ReconsState:
 
     solver = props.solver(props)
     sim = solver.init(sim)
-    groups = GroupManager(sim.state.scan, props.grouping, props.compact, seed=seed)
+    groups = GroupManager(sim.state.scan.data, props.grouping, props.compact, seed=seed)
 
     calc_error_mask = mask_fraction_of_groups(len(groups), props.calc_error_fraction)
 
@@ -79,7 +79,7 @@ def run_engine(args: EngineArgs, props: ConventionalEnginePlan) -> ReconsState:
 
     # runs rescaling
     sim = solver.presolve(
-        sim, groups.iter(sim.state.scan),
+        sim, groups.iter(sim.state.scan.data),
         patterns=patterns, pattern_mask=pattern_mask,
         propagators=propagators
     )
@@ -95,7 +95,7 @@ def run_engine(args: EngineArgs, props: ConventionalEnginePlan) -> ReconsState:
         iter_shuffle_groups = shuffle_groups({'state': sim.state, 'niter': props.niter})
 
         sim, pos_update, group_errors = solver.run_iteration(
-            sim, groups.iter(sim.state.scan, i, iter_shuffle_groups),
+            sim, groups.iter(sim.state.scan.data, i, iter_shuffle_groups),
             patterns=patterns, pattern_mask=pattern_mask, propagators=propagators,
             update_object=update_object({'state': sim.state, 'niter': props.niter}),
             update_probe=update_probe({'state': sim.state, 'niter': props.niter}),
@@ -116,16 +116,16 @@ def run_engine(args: EngineArgs, props: ConventionalEnginePlan) -> ReconsState:
 
             # subtract mean position update
             pos_update -= xp.mean(pos_update, tuple(range(pos_update.ndim - 1)))
-            pos_update, position_solver_state = position_solver.perform_update(sim.state.scan, pos_update, position_solver_state)
+            pos_update, position_solver_state = position_solver.perform_update(sim.state.scan.data, pos_update, position_solver_state)
             # subtract mean again (this can change with momentum)
             pos_update -= xp.mean(pos_update, tuple(range(pos_update.ndim - 1)))
             pos_update_rms = float(xp.mean(xp.linalg.norm(pos_update, axis=-1, keepdims=True)))
             logger.info(f"Position update: mean {pos_update_rms}")
-            sim.state.scan += pos_update
-            assert_dtype(sim.state.scan, dtype)
+            sim.state.scan.data += pos_update
+            assert_dtype(sim.state.scan.data, dtype)
 
             # check positions are at least overlapping object
-            sim.state.object.sampling.check_scan(sim.state.scan, sim.state.probe.sampling.extent / 2.)
+            sim.state.object.sampling.check_scan(sim.state.scan.data, sim.state.probe.sampling.extent / 2.)
 
             progress['pos_update_rms'].iters.append(i + start_i)
             progress['pos_update_rms'].values.append(pos_update_rms)
@@ -139,7 +139,7 @@ def run_engine(args: EngineArgs, props: ConventionalEnginePlan) -> ReconsState:
                 progress[k].values.append(error)
 
         sim.state.progress = progress
-        observer.update_iteration(sim.state, i, props.niter, {'total_loss': error})
+        observer.update_iteration(sim.state, i, props.niter, {'total_loss': error} if error is not None else {})
 
     observer.finish_engine(sim.state)
     return sim.state
\ No newline at end of file
diff --git a/phaser/engines/conventional/solvers.py b/phaser/engines/conventional/solvers.py
index b442b35..117eef6 100644
--- a/phaser/engines/conventional/solvers.py
+++ b/phaser/engines/conventional/solvers.py
@@ -102,7 +102,7 @@ def run_iteration(
 
         new_obj_mag = xp.zeros_like(self.obj_mag)
         new_probe_mag = xp.zeros_like(self.probe_mag)
-        pos_update = xp.zeros_like(sim.state.scan, dtype=sim.dtype)
+        pos_update = xp.zeros_like(sim.state.scan.data, dtype=sim.dtype)
         iter_errors = []
 
         for (group_i, (group, group_patterns)) in enumerate(self.iter_patterns(groups, patterns, xp)):
@@ -171,7 +171,7 @@ def run_slice(slice_i: int, prop: t.Optional[NDArray[numpy.complexfloating]], st
         return (probe_mag, psi)
 
     props = tilt_propagators(sim.ky, sim.kx, sim.state, props,
-                             sim.state.tilt[tuple(group)] if sim.state.tilt is not None else None)
+                             sim.state.scan.tilt[tuple(group)] if sim.state.scan.tilt is not None else None)
     (probe_mag, psi) = slice_forwards(props, (probe_mag, psi), run_slice)
 
     # modeled and experimental intensity
@@ -238,7 +238,7 @@ def sim_slice(slice_i: int, prop: t.Optional[NDArray[numpy.complexfloating]], st
         return (group_probe_mag, psi)
 
     props = tilt_propagators(sim.ky, sim.kx, sim.state, props,
-                             sim.state.tilt[tuple(group)] if sim.state.tilt is not None else None)
+                             sim.state.scan.tilt[tuple(group)] if sim.state.scan.tilt is not None else None)
     (group_probe_mag, psi) = slice_forwards(props, (group_probe_mag, psi), sim_slice, jit_unroll_slices=jit_unroll_slices)
 
     new_obj_mag += group_obj_mag
@@ -377,7 +377,7 @@ def run_iteration(
         xp = sim.xp
 
         # TODO: ePIE position update
-        pos_update = xp.zeros_like(sim.state.scan)
+        pos_update = xp.zeros_like(sim.state.scan.data)
         iter_errors = []
 
         beta_object = process_schedule(self.plan.beta_object)({'state': sim.state, 'niter': self.engine_plan.niter})
@@ -432,7 +432,7 @@ def run_slice(slice_i: int, prop: t.Optional[NDArray[numpy.complexfloating]], ps
         return psi
 
     props = tilt_propagators(sim.ky, sim.kx, sim.state, props,
-                             sim.state.tilt[tuple(group)] if sim.state.tilt is not None else None)
+                             sim.state.scan.tilt[tuple(group)] if sim.state.scan.tilt is not None else None)
     psi = slice_forwards(props, psi, run_slice)
 
     # modeled and experimental intensity
@@ -479,7 +479,7 @@ def sim_slice(slice_i: int, prop: t.Optional[NDArray[numpy.complexfloating]], ps
         return psi
 
     props = tilt_propagators(sim.ky, sim.kx, sim.state, props,
-                             sim.state.tilt[tuple(group)] if sim.state.tilt is not None else None)
+                             sim.state.scan.tilt[tuple(group)] if sim.state.scan.tilt is not None else None)
     psi = slice_forwards(props, psi, sim_slice, jit_unroll_slices=jit_unroll_slices)
 
     model_wave = fft2(psi[-1] * group_obj[:, -1, None])
diff --git a/phaser/engines/gradient/run.py b/phaser/engines/gradient/run.py
index 2efbd46..018ecdf 100644
--- a/phaser/engines/gradient/run.py
+++ b/phaser/engines/gradient/run.py
@@ -26,6 +26,7 @@
 
 logger = logging.getLogger(__name__)
 _PER_ITER_VARS: t.FrozenSet[ReconsVar] = frozenset({'positions', 'tilt'})
+_PER_ITER_PATHS: t.FrozenSet[str] = frozenset({'initial'}) | _PER_ITER_VARS
 
 
 def process_solvers(
@@ -67,11 +68,12 @@ def process_solvers(
     )
 
 
-_PATH_MAP: t.Dict[t.Tuple[str, ...], ReconsVar] = {
+_PATH_MAP: t.Dict[t.Tuple[str, ...], str] = {
     ('object', 'data'): 'object',
     ('probe', 'data'): 'probe',
-    ('scan',): 'positions',
-    ('tilt',): 'tilt'
+    ('scan', 'data'): 'positions',
+    ('scan', 'tilt'): 'tilt',
+    ('scan', 'initial'): 'initial',  # not a solver variable, but need to apply group indexing
 }
 
 def _normalize_path(path: t.Tuple[tree.GetAttrKey, ...]) -> t.Tuple[str, ...]:
@@ -97,7 +99,7 @@ def f(path: t.Tuple[tree.GetAttrKey, ...], val: t.Any):
         if (var := _PATH_MAP.get(_normalize_path(path))):
             if var in vars:
                 return vars[var]
-            if var in _PER_ITER_VARS and val is not None and group is not None:
+            if var in _PER_ITER_PATHS and val is not None and group is not None:
                 return val[tuple(group)]
         return val
 
@@ -110,13 +112,13 @@ def apply_update(state: ReconsState, update: t.Dict[ReconsVar, numpy.ndarray]) -
     if 'object' in update:
         state.object.data += update['object']
     if 'tilt' in update:
-        state.tilt += update['tilt']
+        state.scan.tilt += update['tilt']
     if 'positions' in update:
         # subtract mean position update
         xp = get_array_module(update['positions'])
         update['positions'] -= xp.mean(update['positions'], tuple(range(update['positions'].ndim - 1)))
 
-        state.scan += update['positions']
+        state.scan.data += update['positions']
 
     return state
 
@@ -177,7 +179,7 @@ def run_engine(args: EngineArgs, props: GradientEnginePlan) -> ReconsState:
     }
     # shuffle_groups defaults to True for sparse groups, False for compact groups
     shuffle_groups = process_flag(props.shuffle_groups or not props.compact)
-    groups = GroupManager(state.scan, props.grouping, props.compact, seed)
+    groups = GroupManager(state.scan.data, props.grouping, props.compact, seed)
 
     observer.init_engine(
         state, recons_name=args['recons_name'],
@@ -210,7 +212,7 @@ def iter_patterns(groups: t.Iterable[NDArray[numpy.int_]]) -> t.Iterable[t.Tuple
 
     # runs rescaling
     rescale_factors = []
-    for (group_i, (group, group_patterns)) in enumerate(iter_patterns(groups.iter(state.scan))):
+    for (group_i, (group, group_patterns)) in enumerate(iter_patterns(groups.iter(state.scan.data))):
         group_rescale_factors = dry_run(
             state, group, propagators, group_patterns,
             xp=xp, dtype=dtype,
@@ -274,7 +276,7 @@ def iter_patterns(groups: t.Iterable[NDArray[numpy.int_]]) -> t.Iterable[t.Tuple
             for (solver, solver_state) in zip(iter_solvers, iter_solver_states)
         ]
 
-        for (group_i, (group, group_patterns)) in enumerate(iter_patterns(groups.iter(state.scan, i, iter_shuffle_groups))):
+        for (group_i, (group, group_patterns)) in enumerate(iter_patterns(groups.iter(state.scan.data, i, iter_shuffle_groups))):
             # prevent the loop running ahead of the GPU stream
             block_until_ready(losses_gpu['total_loss'])
 
@@ -341,8 +343,8 @@ def iter_patterns(groups: t.Iterable[NDArray[numpy.int_]]) -> t.Iterable[t.Tuple
 
         if 'positions' in iter_vars:
             # check positions are at least overlapping object
-            state.object.sampling.check_scan(state.scan, state.probe.sampling.extent / 2.)
-            assert_dtype(state.scan, dtype)
+            state.object.sampling.check_scan(state.scan.data, state.probe.sampling.extent / 2.)
+            assert_dtype(state.scan.data, dtype)
 
         state.progress = progress
         observer.update_iteration(state, i, props.niter, losses)
@@ -434,8 +436,8 @@ def run_model(
 ) -> t.Tuple[Float, t.Tuple[SolverStates, t.Dict[str, Float]]]:
     # apply vars to simulation
     sim = insert_vars(vars, sim, group)
-    group_scan = sim.scan
-    group_tilts = sim.tilt
+    group_scan = sim.scan.data
+    group_tilts = sim.scan.tilt
 
     (ky, kx) = sim.probe.sampling.recip_grid(dtype=dtype, xp=xp)
     xp = get_array_module(sim.probe.data)
@@ -489,7 +491,8 @@ def dry_run(
     dtype: t.Type[numpy.floating],
 ) -> NDArray[numpy.floating]:
     (ky, kx) = sim.probe.sampling.recip_grid(dtype=dtype, xp=xp)
-    group_scan = sim.scan[tuple(group)]
+    group_scan = sim.scan.data[tuple(group)]
+    group_tilt = sim.scan.tilt[tuple(group)] if sim.scan.tilt is not None else None
 
     probes = ifft2shift(sim.probe.data)
     group_obj = ifft2shift(sim.object.sampling.get_view_at_pos(sim.object.data, group_scan, probes.shape[-2:]))
@@ -501,7 +504,7 @@ def sim_slice(slice_i: int, prop: t.Optional[NDArray[numpy.complexfloating]], ps
             return ifft2(fft2(psi * group_obj[:, slice_i, None], shift=False) * prop[:, None], shift=False)
         return psi * group_obj[:, slice_i, None]
 
-    t_props = tilt_propagators(ky, kx, sim, props, sim.tilt[tuple(group)] if sim.tilt is not None else None)
+    t_props = tilt_propagators(ky, kx, sim, props, group_tilt)
     model_wave = fft2(slice_forwards(t_props, probes, sim_slice), shift=False)
     model_intensity = xp.sum(abs2(model_wave), axis=(1, -2, -1))
     exp_intensity = xp.sum(group_patterns, axis=(-2, -1))
diff --git a/phaser/execute.py b/phaser/execute.py
index 62fdae8..59374c7 100644
--- a/phaser/execute.py
+++ b/phaser/execute.py
@@ -5,16 +5,17 @@
 import sys
 import typing as t
 
+from frozendict import frozendict
 import numpy
 import pane
 
 from phaser.types import EarlyTermination
 from phaser.utils.num import Device, cast_array_module, get_array_module, get_backend_devices, get_backend_module, set_default_device, to_device, xp_is_jax, Sampling, to_complex_dtype, xp_is_torch
 from phaser.utils.object import ObjectSampling
-from phaser.utils.misc import unwrap
+from phaser.utils.misc import freeze, unwrap
 from .hooks import EngineHook, Hook, ObjectHook, RawData
 from .plan import GradientEnginePlan, ReconsPlan, EnginePlan, ScanHook, ProbeHook, TiltHook
-from .state import Patterns, ReconsState, PartialReconsState, IterState, PreparedRecons
+from .state import Patterns, ReconsState, PartialReconsState, IterState, PreparedRecons, ScanState
 from .observer import Observer, LoggingObserver, PatienceObserver, SaveObserver, ObserverSet
 from .version import version_info
 
@@ -135,11 +136,13 @@ def _normalize_scan_shape(
     Normalizes 'patterns' and 'state' to share a common scan shape.
 
     Requires that there are an equal number of patterns and scan positions.
-    Reshapes 'state.scan' and 'patterns' to match shape, choosing the highest
-    dimensional shape of the two. 'state.tilt' is reshaped as well.
+    Reshapes 'state.scan' (scan and tilt) and 'patterns' to match shape,
+    choosing the highest dimensional shape of the two.
     """
+    scan = state.scan
+
     patterns_shape = patterns.patterns.shape[:-2]
-    scan_shape = state.scan.shape[:-1]
+    scan_shape = scan.data.shape[:-1]
 
     n_patterns = math.prod(patterns_shape)
     n_scan = math.prod(scan_shape)
@@ -150,14 +153,22 @@ def _normalize_scan_shape(
     new_shape = scan_shape if len(scan_shape) > len(patterns_shape) else patterns_shape
 
     patterns.patterns = patterns.patterns.reshape((*new_shape, *patterns.patterns.shape[-2:]))
-    state.scan = state.scan.reshape((*new_shape, 2))
+    scan.data = scan.data.reshape((*new_shape, 2))
+    scan.initial = scan.initial.reshape((*new_shape, 2))
 
-    if state.tilt is not None:
-        n_tilt = math.prod(state.tilt.shape[:-1])
+    if scan.tilt is not None:
+        n_tilt = math.prod(scan.tilt.shape[:-1])
         if n_tilt != n_patterns:
             raise ValueError(f"# of tilt positions {n_scan} doesn't match # of patterns {n_patterns}")
 
-        state.tilt = state.tilt.reshape((*new_shape, 2))
+        scan.tilt = scan.tilt.reshape((*new_shape, 2))
+
+    # also normalize raster_rows and raster_cols metadata
+    scan.meta = frozendict(
+        scan.meta,
+        **({'raster_rows': freeze(numpy.array(scan.meta['raster_rows']).reshape(new_shape))} if 'raster_rows' in scan.meta else {}),
+        **({'raster_cols': freeze(numpy.array(scan.meta['raster_cols']).reshape(new_shape))} if 'raster_cols' in scan.meta else {}),
+    )
 
     return patterns, state
 
@@ -309,27 +320,27 @@ def initialize_reconstruction(
 
     if init_state.scan is not None and plan.init.scan is None:
         logging.info("Re-using scan from initial state...")
-        scan = init_state.scan
+        scan = init_state.scan.copy()
+        scan.data = scan.data.astype(dtype)
+        scan.initial = scan.initial.astype(dtype)
     else:
         logging.info("Initializing scan...")
         scan = pane.from_data(scan_hook, ScanHook)(  # type: ignore
             {'dtype': dtype, 'seed': seed, 'xp': xp}
         )
 
-    if init_state.tilt is not None and plan.init.tilt is None:
+    if init_state.scan is not None and init_state.scan.tilt is not None and plan.init.tilt is None:
         logging.info("Re-using tilt from initial state...")
-        tilt = init_state.tilt
+        scan.tilt = init_state.scan.tilt.astype(dtype)
     elif tilt_hook is not None:
         logging.info("Initializing tilt...")
-        tilt = pane.from_data(tilt_hook, TiltHook)(  # type: ignore
-            {'dtype': dtype, 'xp': xp, 'shape': scan.shape[:-1]}
+        scan.tilt = pane.from_data(tilt_hook, TiltHook)(  # type: ignore
+            {'dtype': dtype, 'xp': xp, 'shape': scan.data.shape[:-1]}
         )
-    else:
-        tilt = None
 
     obj_pad_px: float = plan.engines[0].obj_pad_px if len(plan.engines) > 0 else 5.0  # type: ignore
     obj_sampling = ObjectSampling.from_scan(
-        scan, sampling.sampling, sampling.extent / 2. + obj_pad_px * sampling.sampling
+        scan.data, sampling.sampling, sampling.extent / 2. + obj_pad_px * sampling.sampling
     )
 
     if init_state.object is not None and plan.init.object is None:
@@ -352,7 +363,6 @@ def initialize_reconstruction(
         probe=probe,
         object=obj,
         scan=scan,
-        tilt=tilt,
         wavelength=wavelength
     )
     state = state.to_xp(xp)  # TODO: figure out why this isn't already the case
@@ -414,7 +424,7 @@ def prepare_for_engine(patterns: Patterns, state: ReconsState, xp: t.Any, engine
         obj_sampling = obj_sampling.with_sampling(state.probe.sampling.sampling)
 
     obj_sampling_pad = obj_sampling.expand_to_scan(
-        state.scan, state.probe.sampling.extent / 2. + engine.obj_pad_px * state.probe.sampling.sampling
+        state.scan.data, state.probe.sampling.extent / 2. + engine.obj_pad_px * state.probe.sampling.sampling
     )
 
     if obj_sampling_pad != obj_sampling:
@@ -448,9 +458,9 @@ def prepare_for_engine(patterns: Patterns, state: ReconsState, xp: t.Any, engine
 
     if isinstance(engine, GradientEnginePlan):
         solver_vars = set(itertools.chain.from_iterable(engine.solvers.keys()))
-        if 'tilt' in solver_vars and state.tilt is None:
+        if 'tilt' in solver_vars and state.scan.tilt is None:
             logging.info("Creating new, zeroed tilt map...")
-            state.tilt = xp.zeros_like(state.scan)
+            state.scan.tilt = xp.zeros_like(state.scan.data)
 
     return patterns, state
 
diff --git a/phaser/hooks/__init__.py b/phaser/hooks/__init__.py
index f861270..c2554f8 100644
--- a/phaser/hooks/__init__.py
+++ b/phaser/hooks/__init__.py
@@ -1,82 +1,89 @@
-from pathlib import Path
 import typing as t
+from pathlib import Path
 
 import numpy
-from numpy.typing import NDArray, DTypeLike
-import pane.annotations as annotations
+from numpy.typing import DTypeLike, NDArray
+from pane import annotations
 from typing_extensions import NotRequired
 
-from ..types import Dataclass, Slices, Aberration
+from ..types import Aberration, Dataclass, Slices
 from .hook import Hook
 
 if t.TYPE_CHECKING:
     from phaser.utils.num import Sampling
     from phaser.utils.object import ObjectSampling
-    from ..state import ObjectState, ProbeState, ReconsState, Patterns  # noqa: F401
+
     from ..execute import Observer
+    from ..state import (  # noqa: F401
+        ObjectState,
+        Patterns,
+        ProbeState,
+        ReconsState,
+        ScanState,
+    )
 
 
 class RawData(t.TypedDict):
     patterns: NDArray[numpy.floating]
     mask: NDArray[numpy.floating]
     sampling: 'Sampling'
-    wavelength: NotRequired[t.Optional[float]]
-    scan_hook: NotRequired[t.Union[t.Dict[str, t.Any], None]]
-    tilt_hook: NotRequired[t.Union[t.Dict[str, t.Any], None]]
-    probe_hook: NotRequired[t.Union[t.Dict[str, t.Any], None]]
-    seed: NotRequired[t.Optional[object]]
+    wavelength: NotRequired[float | None]
+    scan_hook: NotRequired[dict[str, t.Any] | None]
+    tilt_hook: NotRequired[dict[str, t.Any] | None]
+    probe_hook: NotRequired[dict[str, t.Any] | None]
+    seed: NotRequired[object | None]
 
 
 class LoadEmpadProps(Dataclass):
     path: Path
 
-    diff_step: t.Optional[float] = None
-    kv: t.Optional[float] = None
-    adu: t.Optional[float] = None
-    det_flips: t.Optional[t.Tuple[bool, bool, bool]] = None
+    diff_step: float | None = None
+    kv: float | None = None
+    adu: float | None = None
+    det_flips: tuple[bool, bool, bool] | None = None
 
 
 class LoadGatanProps(Dataclass):
     path: Path
     
-    diff_step: t.Optional[float] = None
-    kv: t.Optional[float] = None
-    adu: t.Optional[float] = None
+    diff_step: float | None = None
+    kv: float | None = None
+    adu: float | None = None
 
 class LoadNionProps(Dataclass):
     path: Path
 
     diff_step: float
-    detector_rotation_offset: t.Optional[float] = None
+    detector_rotation_offset: float | None = None
 
 class LoadManualProps(Dataclass, kw_only=True):
     path: Path
 
-    det_shape: t.Optional[t.Tuple[int, int]] = None
+    det_shape: tuple[int, int] | None = None
     """Detector shape `(ny, nx)` (after flips are applied). Required when loading raw binary files, optional otherwise."""
-    dtype: t.Optional[str] = None
+    dtype: str | None = None
     """Numpy dtype to load (e.g. 'float32'). Applies only when loading raw binary files."""
     gap: int = 0
     """Gap (in bytes) between patterns in the file. Applies only when loading raw binary files."""
     offset: int = 0
     """Offset (in bytes) before start of patterns in the file. Applies only when loading raw binary files."""
 
-    key: t.Optional[str] = None
+    key: str | None = None
     """Key to load from HDF5 or mat file (ex. 'raw.patterns.data')"""
 
     diff_step: float
     # TODO: post-validate (one of kv or wavelength must be specified)
-    kv: t.Optional[float] = None
-    wavelength: t.Optional[float] = None
-    adu: t.Optional[float] = None
+    kv: float | None = None
+    wavelength: float | None = None
+    adu: float | None = None
     """Detector ADU, representing the single-particle signal. Used to scale patterns."""
 
-    det_flips: t.Optional[t.Tuple[bool, bool, bool]] = None
+    det_flips: tuple[bool, bool, bool] | None = None
     fftshifted: bool = False
     """Whether patterns are fftshifted (zero-frequency in corner of array)"""
 
 class RawDataHook(Hook[None, RawData]):
-    known = {
+    known: t.ClassVar = {
         'empad': ('phaser.hooks.io.empad:load_empad', LoadEmpadProps),
         'gatan': ('phaser.hooks.io.gatan:load_gatan', LoadGatanProps, ('rsciio',)),
         'nion': ('phaser.hooks.io.nion:load_nion', LoadNionProps),
@@ -87,19 +94,19 @@ class RawDataHook(Hook[None, RawData]):
 class ProbeHookArgs(t.TypedDict):
     sampling: 'Sampling'
     wavelength: float
-    seed: t.Optional[object]
+    seed: object | None
     dtype: DTypeLike
     xp: t.Any
 
 
 class FocusedProbeProps(Dataclass):
-    defocus: t.Optional[float] = None  # defocus, + is overfocus [A]
-    conv_angle: t.Optional[float] = None  # semiconvergence angle [mrad]
+    defocus: float | None = None  # defocus, + is overfocus [A]
+    conv_angle: float | None = None  # semiconvergence angle [mrad]
     aberrations: t.Sequence[Aberration] = ()
 
 
 class ProbeHook(Hook[ProbeHookArgs, 'ProbeState']):
-    known = {
+    known: t.ClassVar = {
         'focused': ('phaser.hooks.probe:focused_probe', FocusedProbeProps),
     }
 
@@ -107,8 +114,8 @@ class ProbeHook(Hook[ProbeHookArgs, 'ProbeState']):
 class ObjectHookArgs(t.TypedDict):
     sampling: 'ObjectSampling'
     wavelength: float
-    slices: t.Optional[Slices]
-    seed: t.Optional[object]
+    slices: Slices | None
+    seed: object | None
     dtype: DTypeLike
     xp: t.Any
 
@@ -118,26 +125,26 @@ class RandomObjectProps(Dataclass):
 
 
 class ObjectHook(Hook[ObjectHookArgs, 'ObjectState']):
-    known = {
+    known: t.ClassVar = {
         'random': ('phaser.hooks.object:random_object', RandomObjectProps),
     }
 
 
 class ScanHookArgs(t.TypedDict):
-    seed: t.Optional[object]
+    seed: object | None
     dtype: DTypeLike
     xp: t.Any
 
 
 class RasterScanProps(Dataclass):
-    shape: t.Optional[t.Tuple[int, int]] = None  # ny, nx (total shape)
-    step_size: t.Union[None, float, t.Tuple[float, float]] = None  # A
-    rotation: t.Optional[float] = None     # degrees CCW
-    affine: t.Optional[t.Annotated[NDArray[numpy.floating], annotations.shape((2, 2))]] = None
+    shape: tuple[int, int] | None = None  # ny, nx (total shape)
+    step_size: None | float | tuple[float, float] = None  # A
+    rotation: float | None = None     # degrees CCW
+    affine: t.Annotated[NDArray[numpy.floating], annotations.shape((2, 2))] | None = None
 
 
-class ScanHook(Hook[ScanHookArgs, NDArray[numpy.floating]]):
-    known = {
+class ScanHook(Hook[ScanHookArgs, 'ScanState']):
+    known: t.ClassVar = {
         'raster': ('phaser.hooks.scan:raster_scan', RasterScanProps),
     }
 
@@ -145,7 +152,7 @@ class ScanHook(Hook[ScanHookArgs, NDArray[numpy.floating]]):
 class TiltHookArgs(t.TypedDict):
     dtype: DTypeLike
     xp: t.Any
-    shape: t.Tuple[int, ...]  # To match raster scan shape
+    shape: tuple[int, ...]  # To match raster scan shape
 
 
 class GlobalTiltProps(Dataclass):
@@ -162,7 +169,7 @@ class CustomTiltProps(Dataclass):
 
 
 class TiltHook(Hook[TiltHookArgs, NDArray[numpy.floating]]):
-    known = {
+    known: t.ClassVar = {
         'global': ('phaser.hooks.tilt:generate_global_tilt', GlobalTiltProps),
         'custom': ('phaser.hooks.tilt:load_custom_tilt', CustomTiltProps),
     }
@@ -171,7 +178,7 @@ class TiltHook(Hook[TiltHookArgs, NDArray[numpy.floating]]):
 class PostInitArgs(t.TypedDict):
     data: 'Patterns'
     state: 'ReconsState'
-    seed: t.Optional[object]
+    seed: object | None
     dtype: DTypeLike
     xp: t.Any
 
@@ -188,15 +195,15 @@ class BinProps(Dataclass):
 
 
 class CropDataProps(Dataclass):
-    crop: t.Tuple[
+    crop: tuple[
         # y_i, y_f, x_i, x_f
-        t.Optional[int], t.Optional[int], t.Optional[int], t.Optional[int],
+        int | None, int | None, int | None, int | None,
     ] 
 
 
 class PoissonProps(Dataclass):
-    scale: t.Optional[float] = None
-    gaussian: t.Optional[float] = 1.0e-3
+    scale: float | None = None
+    gaussian: float | None = 1.0e-3
 
 
 class DropNanProps(Dataclass):
@@ -208,7 +215,7 @@ class DiffractionAlignProps(Dataclass):
 
 
 class PostLoadHook(Hook[RawData, RawData]):
-    known = {
+    known: t.ClassVar = {
         'crop_data': ('phaser.hooks.preprocessing:crop_data', CropDataProps),
         'poisson': ('phaser.hooks.preprocessing:add_poisson_noise', PoissonProps),
         'scale': ('phaser.hooks.preprocessing:scale_patterns', ScaleProps),
@@ -217,8 +224,8 @@ class PostLoadHook(Hook[RawData, RawData]):
     }
 
 
-class PostInitHook(Hook[PostInitArgs, t.Tuple['Patterns', 'ReconsState']]):
-    known = {
+class PostInitHook(Hook[PostInitArgs, tuple['Patterns', 'ReconsState']]):
+    known: t.ClassVar = {
         'drop_nans': ('phaser.hooks.preprocessing:drop_nan_patterns', DropNanProps),
         'diffraction_align': ('phaser.hooks.preprocessing:diffraction_align', DiffractionAlignProps),
     }
@@ -227,7 +234,7 @@ class PostInitHook(Hook[PostInitArgs, t.Tuple['Patterns', 'ReconsState']]):
 class EngineArgs(t.TypedDict):
     data: 'Patterns'
     state: 'ReconsState'
-    dtype: t.Type[numpy.floating]
+    dtype: type[numpy.floating]
     xp: t.Any
     recons_name: str
     observer: 'Observer'
@@ -235,4 +242,4 @@ class EngineArgs(t.TypedDict):
 
 
 class EngineHook(Hook[EngineArgs, 'ReconsState']):
-    known = {}  # filled in by plan.py
+    known: t.ClassVar = {}  # filled in by plan.py
diff --git a/phaser/hooks/object.py b/phaser/hooks/object.py
index 2dac79d..7371311 100644
--- a/phaser/hooks/object.py
+++ b/phaser/hooks/object.py
@@ -3,11 +3,12 @@
 
 from phaser.utils.num import cast_array_module, to_complex_dtype
 from phaser.utils.object import random_phase_object
-from ..state import ObjectState
+
+from ..state import PixelatedObjectState
 from . import ObjectHookArgs, RandomObjectProps
 
 
-def random_object(args: ObjectHookArgs, props: RandomObjectProps) -> ObjectState:
+def random_object(args: ObjectHookArgs, props: RandomObjectProps) -> PixelatedObjectState:
     sampling = args['sampling']
 
     if args['slices'] is not None:
@@ -22,4 +23,4 @@ def random_object(args: ObjectHookArgs, props: RandomObjectProps) -> ObjectState
         dtype=to_complex_dtype(args['dtype']),
         xp=cast_array_module(args['xp'])
     )
-    return ObjectState(sampling, obj, thicknesses)
\ No newline at end of file
+    return PixelatedObjectState(sampling, obj, thicknesses)
\ No newline at end of file
diff --git a/phaser/hooks/preprocessing.py b/phaser/hooks/preprocessing.py
index 14ba218..4b46967 100644
--- a/phaser/hooks/preprocessing.py
+++ b/phaser/hooks/preprocessing.py
@@ -5,12 +5,22 @@
 import numpy
 from numpy.typing import NDArray
 
+from phaser.state import Patterns, ReconsState, ScanState
 from phaser.types import cast_length
-from phaser.utils.num import get_array_module, cast_array_module, to_numpy, Sampling
-from phaser.utils.misc import create_rng, create_sparse_groupings
 from phaser.utils.image import affine_transform
-from phaser.state import Patterns, ReconsState
-from . import RawData, PostInitArgs, PoissonProps, ScaleProps, DropNanProps, CropDataProps, OffsetProps, BinProps
+from phaser.utils.misc import create_rng, create_sparse_groupings, freeze
+from phaser.utils.num import Sampling, cast_array_module, get_array_module, to_numpy
+
+from . import (
+    BinProps,
+    CropDataProps,
+    DropNanProps,
+    OffsetProps,
+    PoissonProps,
+    PostInitArgs,
+    RawData,
+    ScaleProps,
+)
 
 logger = logging.getLogger(__name__)
 
@@ -24,12 +34,11 @@ def crop_data(raw_data: RawData, props: CropDataProps) -> RawData:
                  f" {0 if x_i is None else x_i}:{raw_data['patterns'].shape[1] if x_f is None else x_f}")
     raw_data['patterns'] = raw_data['patterns'][slice(y_i, y_f), slice(x_i, x_f)]
 
-    if (scan_hook := raw_data.get('scan_hook', None)) is not None:
-        if scan_hook['type'] == 'raster':
-            raw_data['scan_hook'] = {
-                **scan_hook,
-                'shape': raw_data['patterns'].shape[:2],
-            }
+    if (scan_hook := raw_data.get('scan_hook', None)) is not None and scan_hook['type'] == 'raster':
+        raw_data['scan_hook'] = {
+            **scan_hook,
+            'shape': raw_data['patterns'].shape[:2],
+        }
 
     return raw_data
 
@@ -84,9 +93,18 @@ def add_poisson_noise(raw_data: RawData, props: PoissonProps) -> RawData:
 def drop_nan_patterns(args: PostInitArgs, props: DropNanProps) -> t.Tuple[Patterns, ReconsState]:
     xp = get_array_module(args['data'].patterns)
 
+    scan = args['state'].scan
+
     # flatten scan, tilt, and patterns
-    scan = args['state'].scan.reshape(-1, 2)
-    tilt = None if args['state'].tilt is None else args['state'].tilt.reshape(-1, 2)
+    scan_arr = scan.data.reshape(-1, 2)
+    initial_arr = scan.initial.reshape(-1, 2)
+    scan_meta = dict(scan.meta)
+    if 'raster_rows' in scan_meta:
+        scan_meta['raster_rows'] = numpy.array(scan_meta['raster_rows']).ravel()
+    if 'raster_cols' in scan_meta:
+        scan_meta['raster_cols'] = numpy.array(scan_meta['raster_cols']).ravel()
+
+    tilt_arr = None if scan.tilt is None else scan.tilt.reshape(-1, 2)
     patterns = args['data'].patterns.reshape(-1, *args['data'].patterns.shape[-2:])
 
     fraction_nan = xp.sum(xp.isnan(patterns), axis=(-1, -2)) / xp.prod(patterns.shape[-2:])
@@ -97,23 +115,30 @@ def drop_nan_patterns(args: PostInitArgs, props: DropNanProps) -> t.Tuple[Patter
         logger.info(f"Dropping {n}/{patterns.shape[0]} patterns which are at least {props.threshold:.1%} NaN values")
         patterns = patterns[~mask]
 
-        if scan.shape[0] == xp.size(mask):
+        if scan_arr.shape[0] == xp.size(mask):
             # apply mask to scan as well
-            scan = scan[~mask]
-        elif scan.shape[0] != patterns.shape[0]:
-            raise ValueError(f"# of scan positions {scan.shape[0]} doesn't match # of patterns"
+            scan_arr = scan_arr[~mask]
+            initial_arr = initial_arr[~mask]
+            if 'raster_rows' in scan_meta:
+                scan_meta['raster_rows'] = scan_meta['raster_rows'][~mask]
+            if 'raster_cols' in scan_meta:
+                scan_meta['raster_cols'] = scan_meta['raster_cols'][~mask]
+        elif scan_arr.shape[0] != patterns.shape[0]:
+            raise ValueError(f"# of scan positions {scan_arr.shape[0]} doesn't match # of patterns"
                              f" before ({mask.size}) or after ({patterns.shape[0]}) filtering")
-        # otherwise, we assume the mask has already been applied to the scan
-
-        if tilt is not None:
-            if tilt.shape[0] == mask.size:
-                tilt = tilt[~mask]
-            elif tilt.shape[0] != patterns.shape[0]:
-                raise ValueError(f"# of tilt positions {tilt.shape[0]} doesn't match # of patterns"
+        # otherwise, we assume the mask has already been applied to the scan (and metadata)
+
+        # tilt can come from an alternate source, so we need to check it separately
+        if tilt_arr is not None:
+            if tilt_arr.shape[0] == mask.size:
+                tilt_arr = tilt_arr[~mask]
+            elif tilt_arr.shape[0] != patterns.shape[0]:
+                raise ValueError(f"# of tilt positions {tilt_arr.shape[0]} doesn't match # of patterns"
                                 f" before ({mask.size}) or after ({patterns.shape[0]}) filtering")
 
-    args['state'].scan = scan
-    args['state'].tilt = tilt
+    args['state'].scan = ScanState(
+        scan_arr, initial_arr, tilt_arr, freeze(scan_meta)
+    )
     args['data'].patterns = patterns
 
     return (args['data'], args['state'])
diff --git a/phaser/hooks/probe.py b/phaser/hooks/probe.py
index d39e513..32835c0 100644
--- a/phaser/hooks/probe.py
+++ b/phaser/hooks/probe.py
@@ -1,11 +1,12 @@
 import logging
 
 from phaser.utils.optics import make_focused_probe
-from ..state import ProbeState
-from . import ProbeHookArgs, FocusedProbeProps
 
+from ..state import PixelatedProbeState
+from . import FocusedProbeProps, ProbeHookArgs
 
-def focused_probe(args: ProbeHookArgs, props: FocusedProbeProps) -> ProbeState:
+
+def focused_probe(args: ProbeHookArgs, props: FocusedProbeProps) -> PixelatedProbeState:
     logger = logging.getLogger(__name__)
 
     if props.conv_angle is None:
@@ -24,4 +25,4 @@ def focused_probe(args: ProbeHookArgs, props: FocusedProbeProps) -> ProbeState:
         ky, kx, args['wavelength'],
         props.conv_angle, defocus=props.defocus, aberrations=props.aberrations
     )
-    return ProbeState(sampling, probe)
\ No newline at end of file
+    return PixelatedProbeState(sampling, probe)
\ No newline at end of file
diff --git a/phaser/hooks/scan.py b/phaser/hooks/scan.py
index a0b932f..a811fe7 100644
--- a/phaser/hooks/scan.py
+++ b/phaser/hooks/scan.py
@@ -1,14 +1,17 @@
 
 import logging
+
 import numpy
-from numpy.typing import NDArray
+from frozendict import frozendict
 
+from phaser.state import ScanState
 from phaser.utils.num import cast_array_module
 from phaser.utils.scan import make_raster_scan
-from . import ScanHookArgs, RasterScanProps
+
+from . import RasterScanProps, ScanHookArgs
 
 
-def raster_scan(args: ScanHookArgs, props: RasterScanProps) -> NDArray[numpy.floating]:
+def raster_scan(args: ScanHookArgs, props: RasterScanProps) -> ScanState:
     xp = cast_array_module(args['xp'])
     logger = logging.getLogger(__name__)
 
@@ -33,5 +36,13 @@ def raster_scan(args: ScanHookArgs, props: RasterScanProps) -> NDArray[numpy.flo
         props.shape, step_size, rot, affine,
         dtype=args['dtype'], xp=xp,
     )
-
-    return scan
+    ii, jj = numpy.indices(props.shape, dtype=numpy.int64)
+    assert ii.shape == jj.shape == scan.shape[:-1]
+
+    return ScanState(
+        scan, scan.copy(), tilt=None, meta=frozendict(
+            type='raster',
+            raster_rows=tuple(map(tuple, ii.tolist())),
+            raster_cols=tuple(map(tuple, jj.tolist())),
+        )
+    )
diff --git a/phaser/state.py b/phaser/state.py
index 6450788..15b8672 100644
--- a/phaser/state.py
+++ b/phaser/state.py
@@ -1,6 +1,7 @@
 import typing as t
 
 import numpy
+from frozendict import frozendict
 from numpy.typing import NDArray
 from typing_extensions import Self
 
@@ -15,7 +16,7 @@
 
 
 @tree_dataclass
-class Patterns():
+class Patterns:
     patterns: NDArray[numpy.floating]
     """Raw diffraction patterns, with 0-frequency sample in corner"""
     pattern_mask: NDArray[numpy.floating]
@@ -28,7 +29,7 @@ def to_numpy(self) -> Self:
 
 
 @tree_dataclass
-class IterState():
+class IterState:
     engine_num: int
     """Engine number. 1-indexed (0 means before any reconstruction)."""
     engine_iter: int
@@ -36,9 +37,9 @@ class IterState():
     total_iter: int
     """Total iteration number. 1-indexed (0 means before any iterations)."""
 
-    n_engine_iters: t.Optional[int] = None
+    n_engine_iters: int | None = None
     """Total number of iterations in this engine."""
-    n_total_iters: t.Optional[int] = None
+    n_total_iters: int | None = None
     """Total number of iterations in the reconstruction."""
 
     def to_numpy(self) -> Self:
@@ -57,13 +58,16 @@ def empty() -> 'IterState':
         return IterState(0, 0, 0)
 
 
-@tree_dataclass(static_fields=('sampling',))
-class ProbeState():
+@tree_dataclass(static_fields=('sampling', 'meta', 'ty'))
+class PixelatedProbeState:
     sampling: Sampling
     """Probe coordinate system. See `Sampling` for more details."""
     data: NDArray[numpy.complexfloating]
     """Probe wavefunction, in realspace. Shape (modes, y, x)"""
 
+    meta: frozendict[str, t.Any] = field(default_factory=frozendict)
+    ty: t.Literal['pixelated'] = 'pixelated'
+
     def resample(
         self, new_samp: Sampling,
         rotation: float = 0.0,
@@ -76,16 +80,16 @@ def resample(
             order=order,
             mode=mode,
         )
-        return self.__class__(new_samp, new_data)
+        return self.__class__(new_samp, new_data, self.meta)
 
     def to_xp(self, xp: t.Any) -> Self:
         return self.__class__(
-            self.sampling, xp.asarray(self.data)
+            self.sampling, xp.asarray(self.data), self.meta
         )
 
     def to_numpy(self) -> Self:
         return self.__class__(
-            self.sampling, to_numpy(self.data)
+            self.sampling, to_numpy(self.data), self.meta
         )
 
     def copy(self) -> Self:
@@ -93,8 +97,12 @@ def copy(self) -> Self:
         return copy.deepcopy(self)
 
 
-@tree_dataclass(static_fields=('sampling',))
-class ObjectState():
+# discriminated union of probe state types
+ProbeState: t.TypeAlias = PixelatedProbeState
+
+
+@tree_dataclass(static_fields=('sampling', 'meta', 'ty'))
+class PixelatedObjectState:
     sampling: ObjectSampling
     """Object coordinate system. See `ObjectSampling` for more details."""
     data: NDArray[numpy.complexfloating]
@@ -105,14 +113,17 @@ class ObjectState():
     Length < 2 for single slice, equal to the number of slices otherwise.
     """
 
+    meta: frozendict[str, t.Any] = field(default_factory=frozendict)
+    ty: t.Literal['pixelated'] = 'pixelated'
+
     def to_xp(self, xp: t.Any) -> Self:
         return self.__class__(
-            self.sampling, xp.asarray(self.data), xp.asarray(self.thicknesses)
+            self.sampling, xp.asarray(self.data), xp.asarray(self.thicknesses), self.meta,
         )
 
     def to_numpy(self) -> Self:
         return self.__class__(
-            self.sampling, to_numpy(self.data), to_numpy(self.thicknesses)
+            self.sampling, to_numpy(self.data), to_numpy(self.thicknesses), self.meta,
         )
 
     def zs(self) -> NDArray[numpy.floating]:
@@ -126,11 +137,47 @@ def copy(self) -> Self:
         return copy.deepcopy(self)
 
 
+# discriminated union of object state types
+ObjectState: t.TypeAlias = PixelatedObjectState
+
+
+@tree_dataclass(static_fields=('meta',))
+class ScanState:
+    data: NDArray[numpy.floating]
+    """Scan coordinates (y, x), in length units. Shape (..., 2)"""
+    initial: NDArray[numpy.floating]
+    """Inital scan coordinates (y, x), in length units."""
+    tilt: NDArray[numpy.floating] | None = None
+    """Tilt angles (y, x) per scan position, in mrad. Shape (..., 2)"""
+
+    meta: frozendict[str, t.Any] = field(default_factory=frozendict)
+
+    def to_xp(self, xp: t.Any) -> Self:
+        return self.__class__(
+            xp.asarray(self.data),
+            xp.asarray(self.initial),
+            None if self.tilt is None else xp.asarray(self.tilt),
+            self.meta,
+        )
+
+    def to_numpy(self) -> Self:
+        return self.__class__(
+            to_numpy(self.data),
+            to_numpy(self.initial),
+            None if self.tilt is None else to_numpy(self.tilt),
+            self.meta,
+        )
+
+    def copy(self) -> Self:
+        import copy
+        return copy.deepcopy(self)
+
+
 @tree_dataclass
 class ProgressState:
-    iters: t.List[int] = field(default_factory=list)
+    iters: list[int] = field(default_factory=list)
     """Iterations error measurements were taken at."""
-    values: t.List[float] = field(default_factory=list)
+    values: list[float] = field(default_factory=list)
     """Detector error measurements at those iterations"""
 
     def copy(self) -> Self:
@@ -145,19 +192,16 @@ class ReconsState:
 
     probe: ProbeState
     object: ObjectState
-    scan: NDArray[numpy.floating]
-    """Scan coordinates (y, x), in length units. Shape (..., 2)"""
-    tilt: t.Optional[NDArray[numpy.floating]] = None
-    """Tilt angles (y, x) per scan position, in mrad. Shape (..., 2)"""
-    progress: t.Dict[str, ProgressState] = field(default_factory=dict)
+    scan: ScanState
+
+    progress: dict[str, ProgressState] = field(default_factory=dict)
 
     def to_xp(self, xp: t.Any) -> Self:
         return self.__class__(
             iter=self.iter,
             probe=self.probe.to_xp(xp),
             object=self.object.to_xp(xp),
-            scan=xp.asarray(self.scan),
-            tilt=None if self.tilt is None else xp.asarray(self.tilt),
+            scan=self.scan.to_xp(xp),
             progress=self.progress,
             wavelength=self.wavelength,
         )
@@ -167,8 +211,7 @@ def to_numpy(self) -> Self:
             iter=self.iter.to_numpy(),
             probe=self.probe.to_numpy(),
             object=self.object.to_numpy(),
-            scan=to_numpy(self.scan),
-            tilt=None if self.tilt is None else to_numpy(self.tilt),
+            scan=self.scan.to_numpy(),
             progress=self.progress,
             wavelength=float(self.wavelength),
         )
@@ -189,23 +232,20 @@ def read_hdf5(file: 'HdfLike') -> 'ReconsState':
 
 @tree_dataclass(kw_only=True, static_fields=('progress',))
 class PartialReconsState:
-    iter: t.Optional[IterState] = None
-    wavelength: t.Optional[Float] = None
+    iter: IterState | None = None
+    wavelength: Float | None = None
 
-    probe: t.Optional[ProbeState] = None
-    object: t.Optional[ObjectState] = None
-    scan: t.Optional[NDArray[numpy.floating]] = None
-    """Scan coordinates (y, x), in length units. Shape (..., 2)"""
-    tilt: t.Optional[NDArray[numpy.floating]] = None
-    progress: t.Optional[t.Dict[str, ProgressState]] = None
+    probe: ProbeState | None = None
+    object: ObjectState | None = None
+    scan: ScanState | None = None
+    progress: dict[str, ProgressState] | None = None
 
     def to_numpy(self) -> Self:
         return self.__class__(
             iter=self.iter.to_numpy() if self.iter is not None else None,
             probe=self.probe.to_numpy() if self.probe is not None else None,
             object=self.object.to_numpy() if self.object is not None else None,
-            scan=to_numpy(self.scan) if self.scan is not None else None,
-            tilt=to_numpy(self.tilt) if self.tilt is not None else None,
+            scan=self.scan.to_numpy() if self.scan is not None else None,
             wavelength=float(self.wavelength) if self.wavelength is not None else None,
             progress=self.progress,
         )
@@ -222,8 +262,8 @@ def to_complete(self) -> ReconsState:
             wavelength=t.cast(Float, self.wavelength),
             probe=t.cast(ProbeState, self.probe),
             object=t.cast(ObjectState, self.object),
-            scan=t.cast(NDArray[numpy.floating], self.scan),
-            tilt=self.tilt, progress=progress, iter=iter,
+            scan=t.cast(ScanState, self.scan),
+            progress=progress, iter=iter,
         )
 
     def write_hdf5(self, file: 'HdfLike'):
@@ -264,4 +304,19 @@ def with_observer(self, observer: t.Union['Observer', t.Iterable['Observer']]) -
         else:
             observers.extend(observer)
 
-        return self.__class__(self.patterns, self.state, self.name, ObserverSet(observers))
\ No newline at end of file
+        return self.__class__(self.patterns, self.state, self.name, ObserverSet(observers))
+
+
+__all__ = [
+    'IterState',
+    'ObjectState',
+    'PartialReconsState',
+    'Patterns',
+    'PixelatedObjectState',
+    'PixelatedProbeState',
+    'PreparedRecons',
+    'ProbeState',
+    'ProgressState',
+    'ReconsState',
+    'ScanState',
+]
\ No newline at end of file
diff --git a/phaser/utils/io.py b/phaser/utils/io.py
index 6890807..b3ed4fa 100644
--- a/phaser/utils/io.py
+++ b/phaser/utils/io.py
@@ -1,21 +1,33 @@
 import contextlib
-from pathlib import Path
+import json
 import typing as t
+from pathlib import Path
 
+import h5py
 import numpy
+from frozendict import frozendict
 from numpy.typing import NDArray
-import h5py
 
+from phaser.state import (
+    IterState,
+    ObjectState,
+    PartialReconsState,
+    PixelatedObjectState,
+    PixelatedProbeState,
+    ProbeState,
+    ProgressState,
+    ReconsState,
+    ScanState,
+)
 from phaser.utils.num import Sampling, to_numpy
 from phaser.utils.object import ObjectSampling
-from phaser.state import ReconsState, IterState, ProbeState, ObjectState, ProgressState, PartialReconsState
-
+from phaser.utils.misc import freeze
 
-HdfLike: t.TypeAlias = t.Union[h5py.File, str, Path]
+HdfLike: t.TypeAlias = h5py.File | str | Path
 OpenMode: t.TypeAlias = t.Literal['r', 'r+', 'w', 'w-', 'x', 'a']
 DTypeT = t.TypeVar('DTypeT', bound=numpy.generic)
 
-_DTYPE_CATEGORIES: t.Dict[t.Type[numpy.generic], t.Type[numpy.generic]] = {
+_DTYPE_CATEGORIES: dict[type[numpy.generic], type[numpy.generic]] = {
     numpy.bool_: numpy.bool_,
     numpy.float32: numpy.floating,
     numpy.float64: numpy.floating,
@@ -37,7 +49,7 @@
     numpy.unsignedinteger: numpy.unsignedinteger,
 }
 
-_CATEGORY_MIN_DTYPE: t.Dict[t.Type[numpy.generic], t.Type[numpy.generic]] = {
+_CATEGORY_MIN_DTYPE: dict[type[numpy.generic], type[numpy.generic]] = {
     numpy.bool_: numpy.bool_,
     numpy.inexact: numpy.float32,
     numpy.floating: numpy.float32,
@@ -70,7 +82,7 @@ def __enter__(self) -> Path:
 
         return self.out_dir
 
-    def __exit__(self, exc_type: t.Optional[type], exc_value: t.Optional[BaseException], tb: t.Any):
+    def __exit__(self, exc_type: type[BaseException] | None, exc_value: BaseException | None, tb: object):
         if exc_value is None and self.any_output:
             # create finished file
             (self.out_dir / 'finished').touch(mode=0o664)
@@ -94,61 +106,108 @@ def hdf5_read_state(file: HdfLike) -> PartialReconsState:
     file = open_hdf5(file, 'r')
 
     ty = _hdf5_read_string(file, 'type')
-    version = _hdf5_read_string(file, 'version')
+    version_str = _hdf5_read_string(file, 'version')
     if ty != 'phaser_state':
         raise ValueError(f"While reading file '{file.filename}':\nExpected a file of type 'phaser_state', instead got type '{ty}'")
 
-    if _parse_version(version) > (0, 1):
-        raise ValueError(f"While reading file '{file.filename}':\nUnsupported file version '{version}'. Maximum supported version is '0.1'.")
+    version = _parse_version(version_str)
+    # VERSION 0.2:
+    # - moved 'scan' into its own group
+    # - added 'meta' to probe, object, scan (backwards compatible)
+    # - added 'type' to probe, object (backwards compatible)
+    if version > (0, 2):
+        raise ValueError(f"While reading file '{file.filename}':\nUnsupported file version '{version_str}'. Maximum supported version is '0.2'.")
+    read_scan_as_group = version >= (0, 2)
 
     wavelength = _hdf5_read_scalar(file, 'wavelength', numpy.float64) if 'wavelength' in file else None
 
     probe = hdf5_read_probe_state(_assert_group(file['probe'])) if 'probe' in file else None
     obj = hdf5_read_object_state(_assert_group(file['object'])) if 'object' in file else None
     iter = hdf5_read_iter_state(_assert_group(file['iter'])) if 'iter' in file else IterState.empty()
-    scan = numpy.asarray(_hdf5_read_dataset(file, 'scan', numpy.float64)) if 'scan' in file else None
-    tilt = numpy.asarray(_hdf5_read_dataset(file, 'tilt', numpy.float64)) if 'tilt' in file else None
 
-    if tilt is not None and scan is not None:
-        assert tilt.shape == scan.shape
+    if 'scan' not in file:
+        scan = None
+    elif read_scan_as_group:  # new behavior
+        scan = hdf5_read_scan_state(_assert_group(file['scan']))
+    else:  # old behavior
+        scan_arr = _hdf5_read_array(file, 'scan', numpy.float64)
+        tilt_arr = _hdf5_read_array(file, 'tilt', numpy.float64, nullable=True)
+        if tilt_arr is not None:
+            assert tilt_arr.shape == scan_arr.shape
+
+        # use current scan as initial
+        scan = ScanState(scan_arr, initial=scan_arr.copy(), tilt=tilt_arr)
+
     progress = hdf5_read_progress_state(_assert_group(file['progress'])) if 'progress' in file else None
 
     return PartialReconsState(
         wavelength=wavelength, iter=iter, probe=probe,
-        object=obj, scan=scan, tilt=tilt, progress=progress
+        object=obj, scan=scan, progress=progress
     )
 
 
 def hdf5_read_probe_state(group: h5py.Group) -> ProbeState:
-    probes = _hdf5_read_dataset(group, 'data', numpy.complexfloating)
+    ty = _hdf5_read_string(group, 'type', nullable=True)
+    if ty not in (None, 'pixelated'):
+        # currently we only support pixelated probes
+        raise ValueError(f"While reading file '{group.file.filename}':\nUnsupported probe type '{ty}'")
+
+    probes = _hdf5_read_array(group, 'data', numpy.complexfloating)
     assert probes.ndim == 3
 
-    extent = _hdf5_read_dataset_shape(group, 'extent', numpy.float64, (2,))
+    extent = _hdf5_read_array_shape(group, 'extent', numpy.float64, (2,))
     (n_y, n_x) = probes.shape[-2:]
 
-    return ProbeState(
+    meta = hdf5_read_meta(group)
+
+    return PixelatedProbeState(
         Sampling((n_y, n_x), extent=(extent[0], extent[1])),
-        data=probes
+        data=probes, meta=meta,
+    )
+
+
+def hdf5_read_scan_state(group: h5py.Group) -> ScanState:
+    scan = _hdf5_read_array(group, 'data', numpy.float64)
+
+    tilt = _hdf5_read_array(group, 'tilt', numpy.float64, nullable=True)
+    initial = _hdf5_read_array(group, 'initial', numpy.float64, nullable=True)
+
+    if tilt is not None:
+        assert tilt.shape == scan.shape
+    if initial is not None:
+        assert initial.shape == scan.shape
+
+    meta = hdf5_read_meta(group)
+
+    return ScanState(
+        scan, initial=scan.copy() if initial is None else initial, tilt=tilt, meta=meta,
     )
 
 
 def hdf5_read_object_state(group: h5py.Group) -> ObjectState:
-    obj = numpy.asarray(_hdf5_read_dataset(group, 'data', numpy.complexfloating))
+    ty = _hdf5_read_string(group, 'type', nullable=True)
+    if ty not in (None, 'pixelated'):
+        # currently we only support pixelated objects
+        raise ValueError(f"While reading file '{group.file.filename}':\nUnsupported object type '{ty}'")
+
+    obj = _hdf5_read_array(group, 'data', numpy.complexfloating)
     (n_z, n_y, n_x) = obj.shape
-    
-    thicknesses = numpy.asarray(_hdf5_read_dataset(group, 'thicknesses', numpy.floating))
+
+    thicknesses = _hdf5_read_array(group, 'thicknesses', numpy.floating)
     assert thicknesses.ndim == 1
     assert thicknesses.size == n_z if n_z > 1 else thicknesses.size in (0, 1)
 
-    sampling = _hdf5_read_dataset_shape(group, 'sampling', numpy.float64, (2,))
-    corner = _hdf5_read_dataset_shape(group, 'corner', numpy.float64, (2,))
+    sampling = _hdf5_read_array_shape(group, 'sampling', numpy.float64, (2,))
+    corner = _hdf5_read_array_shape(group, 'corner', numpy.float64, (2,))
 
-    region_min = _hdf5_read_dataset_shape(group, 'region_min', numpy.float64, (2,)) if 'region_min' in group else None
-    region_max = _hdf5_read_dataset_shape(group, 'region_max', numpy.float64, (2,)) if 'region_max' in group else None
+    region_min = _hdf5_read_array_shape(group, 'region_min', numpy.float64, (2,), nullable=True)
+    region_max = _hdf5_read_array_shape(group, 'region_max', numpy.float64, (2,), nullable=True)
 
-    return ObjectState(
+    meta = hdf5_read_meta(group)
+
+    return PixelatedObjectState(
         ObjectSampling((n_y, n_x), sampling, corner, region_min, region_max),
-        data=obj, thicknesses=thicknesses
+        data=obj, thicknesses=thicknesses, meta=meta,
     )
 
 
@@ -162,24 +221,24 @@ def hdf5_read_iter_state(group: h5py.Group) -> IterState:
     )
 
 
-def hdf5_read_progress_state(group: h5py.Group) -> t.Dict[str, ProgressState]:
+def hdf5_read_progress_state(group: h5py.Group) -> dict[str, ProgressState]:
     if 'iters' in group and 'detector_errors' in group:
         # read old-style, convert to new style
-        iters = numpy.asarray(_hdf5_read_dataset(group, 'iters', numpy.int64))
-        values = numpy.asarray(_hdf5_read_dataset(group, 'detector_errors', numpy.float64))
+        iters = _hdf5_read_array(group, 'iters', numpy.int64)
+        values = _hdf5_read_array(group, 'detector_errors', numpy.float64)
         assert iters.ndim == values.ndim == 1
         assert iters.shape == values.shape
 
         return {'total_loss': ProgressState(iters.tolist(), values.tolist())}
 
     # read new-style
-    d: t.Dict[str, ProgressState] = {}
+    d: dict[str, ProgressState] = {}
 
-    for (k, group) in group.items():
-        if not isinstance(group, h5py.Group):
+    for (k, inner) in group.items():
+        if not isinstance(inner, h5py.Group):
             continue
-        iters = numpy.asarray(_hdf5_read_dataset(group, 'iters', numpy.int64))
-        values = numpy.asarray(_hdf5_read_dataset(group, 'values', numpy.float64))
+        iters = _hdf5_read_array(inner, 'iters', numpy.int64)
+        values = _hdf5_read_array(inner, 'values', numpy.float64)
         assert iters.ndim == values.ndim == 1
         assert iters.shape == values.shape
 
@@ -188,10 +247,10 @@ def hdf5_read_progress_state(group: h5py.Group) -> t.Dict[str, ProgressState]:
     return d
 
 
-def hdf5_write_state(state: t.Union[ReconsState, PartialReconsState], file: HdfLike):
+def hdf5_write_state(state: ReconsState | PartialReconsState, file: HdfLike):
     file = open_hdf5(file, 'w')  # overwrite if existing
-    file.create_dataset('type', (), h5py.string_dtype(), "phaser_state")
-    file.create_dataset('version', (), h5py.string_dtype(), "0.1")
+    _hdf5_write_string(file, 'type', "phaser_state")
+    _hdf5_write_string(file, 'version', "0.2")
     file.create_dataset('wavelength', (), numpy.float64, state.wavelength)
 
     if state.probe is not None:
@@ -199,9 +258,7 @@ def hdf5_write_state(state: t.Union[ReconsState, PartialReconsState], file: HdfL
     if state.object is not None:
         hdf5_write_object_state(state.object, file.create_group("object"))
     if state.scan is not None:
-        file.create_dataset('scan', data=to_numpy(state.scan).astype(numpy.float64))
-    if state.tilt is not None:
-        file.create_dataset('tilt', data=to_numpy(state.tilt).astype(numpy.float64))
+        hdf5_write_scan_state(state.scan, file.create_group("scan"))
     if state.iter is not None:
         hdf5_write_iter_state(state.iter, file.create_group("iter"))
     if state.progress is not None:
@@ -209,7 +266,11 @@ def hdf5_write_state(state: t.Union[ReconsState, PartialReconsState], file: HdfL
 
 
 def hdf5_write_probe_state(state: ProbeState, group: h5py.Group):
+    # we only support pixelated probes currently
+    assert state.ty == 'pixelated'
     assert state.data.ndim == 3
+
+    _hdf5_write_string(group, 'type', state.ty)
     dataset = group.create_dataset('data', data=to_numpy(state.data))
     dataset.dims[0].label = 'mode'
     dataset.dims[1].label = 'y'
@@ -218,15 +279,21 @@ def hdf5_write_probe_state(state: ProbeState, group: h5py.Group):
     group.create_dataset('sampling', data=state.sampling.sampling.astype(numpy.float64))
     group.create_dataset('extent', data=state.sampling.extent.astype(numpy.float64))
 
+    hdf5_write_meta(group, state.meta)
+
 
 def hdf5_write_object_state(state: ObjectState, group: h5py.Group):
+    # we only support pixelated objects currently
+    assert state.ty == 'pixelated'
     assert state.data.ndim == 3
     assert state.thicknesses.ndim == 1
-    n_z = state.data.shape[0]
 
+    n_z = state.data.shape[0]
     thick = to_numpy(state.thicknesses)
     assert thick.ndim == 1
     assert thick.size == n_z if n_z > 1 else thick.size in (0, 1)
+
+    _hdf5_write_string(group, 'type', state.ty)
     group.create_dataset('thicknesses', data=thick)
     zs = group.create_dataset('zs', data=to_numpy(state.zs()))
     zs.make_scale("z")
@@ -244,6 +311,27 @@ def hdf5_write_object_state(state: ObjectState, group: h5py.Group):
     _hdf5_write_nullable_dataset(group, 'region_min', state.sampling.region_min, numpy.float64)
     _hdf5_write_nullable_dataset(group, 'region_max', state.sampling.region_max, numpy.float64)
 
+    hdf5_write_meta(group, state.meta)
+
+
+def hdf5_write_scan_state(state: ScanState, group: h5py.Group):
+    assert state.data.ndim >= 2
+    assert state.initial.shape == state.data.shape
+    if state.tilt is not None:
+        assert state.tilt.shape == state.data.shape
+
+    dataset = group.create_dataset('data', data=to_numpy(state.data).astype(numpy.float64))
+    dataset.dims[dataset.ndim - 1].label = 'yx'
+
+    dataset = group.create_dataset('initial', data=to_numpy(state.initial).astype(numpy.float64))
+    dataset.dims[dataset.ndim - 1].label = 'yx'
+
+    if state.tilt is not None:
+        dataset = group.create_dataset('tilt', data=to_numpy(state.tilt).astype(numpy.float64))
+        dataset.dims[dataset.ndim - 1].label = 'yx'
+
+    hdf5_write_meta(group, state.meta)
+
 
 def hdf5_write_iter_state(state: IterState, group: h5py.Group):
     group.create_dataset("engine_num", (), numpy.uint64, data=state.engine_num)
@@ -251,7 +339,7 @@ def hdf5_write_iter_state(state: IterState, group: h5py.Group):
     group.create_dataset("total_iter", (), numpy.uint64, data=state.total_iter)
 
 
-def hdf5_write_progress_state(state: t.Dict[str, ProgressState], group: h5py.Group):
+def hdf5_write_progress_state(state: dict[str, ProgressState], group: h5py.Group):
     for (k, v) in state.items():
         subgroup = group.require_group(k)
 
@@ -262,82 +350,178 @@ def hdf5_write_progress_state(state: t.Dict[str, ProgressState], group: h5py.Gro
         dataset.dims[0].attach_scale(iters)
 
 
-def _parse_version(version: str) -> t.Tuple[int, ...]:
+def hdf5_read_meta(group: h5py.Group, path: str = 'meta') -> frozendict[str, t.Any]:
+    if path not in group:
+        return frozendict()
+
+    dataset = group[path]
+    if not isinstance(dataset, h5py.Dataset):
+        raise TypeError(f"While reading '{group.file.filename}:\n"
+                        f"Expected a dataset at path '{group.name}/{path}', instead found {type(dataset)}")
+
+    if h5py.check_string_dtype(dataset.dtype) is None:
+        raise TypeError(f"While reading '{group.file.filename}:\n"
+                        f"Expected a string dataset at path '{group.name}/{path}', instead found {dataset.dtype}") from None
+
+    s = dataset[()]
+
+    if isinstance(s, numpy.ndarray):
+        assert s.shape == ()
+
+    try:
+        s = bytes(s).decode('utf-8')
+    except (TypeError, UnicodeDecodeError) as e:
+        e.add_note(f"While reading '{group.file.filename}', invalid string dataset '{group.name}/{path}'")
+        raise
+
+    try:
+        d = freeze(json.loads(s))
+        assert isinstance(d, frozendict)  # config toplevel should be dict
+        return d
+    except (json.JSONDecodeError, AssertionError) as e:
+        e.add_note(f"While reading '{group.file.filename}', invalid JSON dataset '{group.name}/{path}'")
+        raise
+
+
+def hdf5_write_meta(group: h5py.Group, meta: frozendict[str, t.Any], path: str = 'meta'):
+    if len(meta):
+        s = json.dumps(meta, indent='', ensure_ascii=False).encode('utf-8')
+        group.create_dataset(path, data=s, dtype=h5py.string_dtype('utf-8'))
+
+
+def _parse_version(version: str) -> tuple[int, ...]:
     try:
         return tuple(map(int, version.split(".")))
     except ValueError:
         raise ValueError(f"Unable to parse version '{version}'") from None
 
 
-def _assert_group(group: t.Union[h5py.Group, h5py.Dataset, h5py.Datatype]) -> h5py.Group:
+def _assert_group(group: h5py.Group | h5py.Dataset | h5py.Datatype) -> h5py.Group:
     if isinstance(group, h5py.Group):
         return group
     raise ValueError(f"While reading '{group.file.filename}':\n"
                      f"Expected a group at path '{group.name}', instead found {type(group)}.")
 
 
-def _hdf5_read_dataset(group: h5py.Group, path: str, dtype: t.Type[DTypeT]) -> t.Union[DTypeT, NDArray[DTypeT]]:
+@t.overload
+def _hdf5_read_dataset(group: h5py.Group, path: str, dtype: type[DTypeT], nullable: t.Literal[False] = ...) -> DTypeT | NDArray[DTypeT]: ...
+@t.overload
+def _hdf5_read_dataset(group: h5py.Group, path: str, dtype: type[DTypeT], nullable: bool = ...) -> DTypeT | NDArray[DTypeT] | None: ...
+
+def _hdf5_read_dataset(group: h5py.Group, path: str, dtype: type[DTypeT], nullable: bool = False) -> DTypeT | NDArray[DTypeT] | None:
     dtype_category = _DTYPE_CATEGORIES[dtype]
 
     if path not in group:
+        if nullable:
+            return None
         raise ValueError(f"While reading '{group.file.filename}':\n"
-                         f"Path '{group.name}{path}' not found.")
+                         f"Path '{group.name}/{path}' not found.")
 
     dataset = group[path]
 
     if not isinstance(dataset, h5py.Dataset):
-        raise ValueError(f"While reading '{group.file.filename}':\n"
-                         f"Expected a dataset at path '{group.name}{path}', instead found {type(dataset)}.")
+        raise TypeError(f"While reading '{group.file.filename}':\n"
+                        f"Expected a dataset at path '{group.name}/{path}', instead found {type(dataset)}.")
 
     if not numpy.issubdtype(dataset.dtype, dtype_category):
         raise ValueError(f"While reading '{group.file.filename}':\n"
-                         f"Expected a dataset of dtype '{dtype_category}' at path '{group.name}{path}', instead found {dataset.dtype}.")
+                         f"Expected a dataset of dtype '{dtype_category}' at path '{group.name}/{path}', instead found {dataset.dtype}.")
+
+    if dataset.shape is None:
+        if nullable:
+            return None
+        raise ValueError(f"While reading '{group.file.filename}':\n"
+                         f"Dataset at path '{group.name}/{path}' is empty.")
 
     # ensure promotion is correct. eg dtype = numpy.floating promotes with numpy.float32
     out_dtype = numpy.promote_types(dataset.dtype, _CATEGORY_MIN_DTYPE.get(dtype_category, dtype))
     return dataset[()].astype(out_dtype)
 
 
-def _hdf5_read_dataset_shape(group: h5py.Group, path: str, dtype: t.Type[DTypeT], shape: t.Tuple[int, ...]) -> NDArray[DTypeT]:
-    arr = numpy.asarray(_hdf5_read_dataset(group, path, dtype))
+@t.overload
+def _hdf5_read_array(group: h5py.Group, path: str, dtype: type[DTypeT], nullable: t.Literal[False] = ...) -> NDArray[DTypeT]: ...
+@t.overload
+def _hdf5_read_array(group: h5py.Group, path: str, dtype: type[DTypeT], nullable: bool = ...) -> NDArray[DTypeT] | None: ...
+
+def _hdf5_read_array(group: h5py.Group, path: str, dtype: type[DTypeT], nullable: bool = False) -> NDArray[DTypeT] | None:
+    arr = _hdf5_read_dataset(group, path, dtype, nullable=nullable)
+    if arr is not None:
+        return numpy.asarray(arr)
+
+
+@t.overload
+def _hdf5_read_array_shape(
+    group: h5py.Group, path: str, dtype: type[DTypeT], shape: tuple[int, ...], nullable: t.Literal[False] = ...
+) -> NDArray[DTypeT]: ...
+@t.overload
+def _hdf5_read_array_shape(
+    group: h5py.Group, path: str, dtype: type[DTypeT], shape: tuple[int, ...], nullable: bool = ...
+) -> NDArray[DTypeT] | None: ...
+
+def _hdf5_read_array_shape(
+    group: h5py.Group, path: str, dtype: type[DTypeT], shape: tuple[int, ...], nullable: bool = False
+) -> NDArray[DTypeT] | None:
+    arr = _hdf5_read_dataset(group, path, dtype, nullable)
+    if arr is None:
+        return None
+    arr = numpy.asarray(arr)
+
     if arr.shape != shape:
         raise ValueError(f"While reading '{group.file.filename}':\n"
-                         f"Expected a dataset of shape '{shape}' at path '{group.name}{path}', instead got shape {arr.shape}.")
+                         f"Expected a dataset of shape '{shape}' at path '{group.name}/{path}', instead got shape {arr.shape}.")
     return arr
 
 
-def _hdf5_read_scalar(group: h5py.Group, path: str, dtype: t.Type[DTypeT]) -> DTypeT:
+def _hdf5_read_scalar(group: h5py.Group, path: str, dtype: type[DTypeT]) -> DTypeT:
     arr = _hdf5_read_dataset(group, path, dtype)
     if isinstance(arr, numpy.ndarray):
-        raise ValueError(f"While reading '{group.file.filename}':\n"
+        raise ValueError(f"While reading '{group.file.filename}':\n"  # noqa: TRY004
                          f"Expected a scalar dataset, instead got shape {arr.shape}.")
     return arr
 
 
-def _hdf5_read_string(group: h5py.Group, path: str) -> str:
+@t.overload
+def _hdf5_read_string(group: h5py.Group, path: str, nullable: t.Literal[False] = ...) -> str: ...
+
+@t.overload
+def _hdf5_read_string(group: h5py.Group, path: str, nullable: bool = ...) -> str | None: ...
+
+def _hdf5_read_string(group: h5py.Group, path: str, nullable: bool = False) -> str | None:
     if path not in group:
+        if nullable:
+            return None
         raise ValueError(f"While reading '{group.file.filename}':\n"
-                         f"Path '{group.name}{path}' not found.")
+                         f"Path '{group.name}/{path}' not found.")
 
     dataset = group[path]
 
     if not isinstance(dataset, h5py.Dataset):
+        raise TypeError(f"While reading '{group.file.filename}':\n"
+                         f"Expected a string at path '{group.name}/{path}', instead found {type(dataset)}.")
+
+    if dataset.shape is None:
+        if nullable:
+            return None
         raise ValueError(f"While reading '{group.file.filename}':\n"
-                         f"Expected a string at path '{group.name}{path}', instead found {type(dataset)}.")
+                         f"Dataset at path '{group.name}/{path}' is empty.")
 
     dataset = dataset[()]
     if not isinstance(dataset, bytes):
-        raise ValueError(f"While reading '{group.file.filename}':\n"
-                         f"Expected a scalar string at path '{group.name}{path}', instead found {dataset} (type {type(dataset)}).")
+        raise TypeError(f"While reading '{group.file.filename}':\n"
+                         f"Expected a scalar string at path '{group.name}/{path}', instead found {dataset} (type {type(dataset)}).")
 
     try:
         return dataset.decode('utf-8')
     except ValueError:
         raise ValueError(f"While reading '{group.file.filename}':\n"
-                         f"Invalid string at path '{group.name}{path}")
+                         f"Invalid string at path '{group.name}/{path}")
+
+
+def _hdf5_write_string(group: h5py.Group, name: str, data: str):
+    group.create_dataset(name, (), h5py.string_dtype('utf-8'), data)
 
 
-def _hdf5_write_nullable_dataset(group: h5py.Group, name: str, data: t.Optional[numpy.ndarray], dtype: t.Any):
+def _hdf5_write_nullable_dataset(group: h5py.Group, name: str, data: numpy.ndarray | None, dtype: t.Any):
     if data is not None:
         group.create_dataset(name, data=to_numpy(data).astype(dtype))
     else:
@@ -345,12 +529,12 @@ def _hdf5_write_nullable_dataset(group: h5py.Group, name: str, data: t.Optional[
 
 
 def tiff_write_opts(
-    sampling: t.Union[Sampling, ObjectSampling],
-    corner: t.Optional[NDArray[numpy.floating]] = None, *,
+    sampling: Sampling | ObjectSampling,
+    corner: NDArray[numpy.floating] | None = None, *,
     unit: t.Literal['angstrom'] = 'angstrom',  # other units not yet supported
     n_slices: int = 1,
-    zs: t.Union[t.Sequence[float], NDArray[numpy.floating], None] = None,
-) -> t.Dict[str, t.Any]:
+    zs: t.Sequence[float] | NDArray[numpy.floating] | None = None,
+) -> dict[str, t.Any]:
     if corner is None:
         corner = sampling.corner
 
@@ -383,11 +567,11 @@ def tiff_write_opts(
 
 
 def tiff_write_opts_recip(
-    sampling: t.Union[Sampling, ObjectSampling], *,
+    sampling: Sampling | ObjectSampling, *,
     unit: t.Literal['1/angstrom'] = '1/angstrom',  # other units not yet supported
     n_slices: int = 1,
-    zs: t.Union[t.Sequence[float], NDArray[numpy.floating], None] = None,
-) -> t.Dict[str, t.Any]:
+    zs: t.Sequence[float] | NDArray[numpy.floating] | None = None,
+) -> dict[str, t.Any]:
     z_dict = {}
     if zs is not None:
         n_slices = len(zs)
@@ -417,8 +601,11 @@ def tiff_write_opts_recip(
 
 
 __all__ = [
+    'HdfLike',
+    'OpenMode',
+    'hdf5_read_state',
+    'hdf5_write_state',
     'open_hdf5',
-    'hdf5_read_state', 'hdf5_write_state',
-    'tiff_write_opts', 'tiff_write_opts_recip',
-    'HdfLike', 'OpenMode',
+    'tiff_write_opts',
+    'tiff_write_opts_recip',
 ]
\ No newline at end of file
diff --git a/phaser/utils/misc.py b/phaser/utils/misc.py
index cdd118f..0f48cf2 100644
--- a/phaser/utils/misc.py
+++ b/phaser/utils/misc.py
@@ -1,12 +1,14 @@
 import functools
 import math
-from types import ModuleType
 import typing as t
+from collections.abc import Mapping, Sequence
+from collections.abc import Set as AbstractSet
+from types import ModuleType
 
 import numpy
+from frozendict import frozendict
+from numpy.random import PCG64, BitGenerator, Generator, SeedSequence
 from numpy.typing import NDArray
-from numpy.random import SeedSequence, PCG64, BitGenerator, Generator
-
 
 T = t.TypeVar('T')
 
@@ -222,6 +224,50 @@ def unwrap(val: t.Optional[T]) -> T:
     return val
 
 
+def freeze(obj: object) -> t.Any:
+    """Attempt to freeze an object, making it immutable."""
+    return _freeze(obj, set())
+
+
+def _freeze(obj: object, stack: set[int]) -> t.Any:
+    if id(obj) in stack:
+        raise ValueError(f"Cannot freeze self-referential object {type(obj).__name__}")
+
+    stack.add(id(obj))
+    try:
+        # handle numpy types
+        if isinstance(obj, numpy.generic):
+            return _freeze(obj.item(), stack)
+        if isinstance(obj, numpy.ndarray):
+            return _freeze(obj.tolist(), stack)
+
+        # mappings
+        if isinstance(obj, Mapping):
+            return frozendict((_freeze(k, stack), _freeze(v, stack)) for (k, v) in obj.items())
+
+        # sets
+        if isinstance(obj, AbstractSet):
+            return frozenset(_freeze(v, stack) for v in obj)
+
+        # sequences
+        if isinstance(obj, Sequence):
+            # immutable sequence types
+            if isinstance(obj, (str, bytes, range)):
+                return obj
+            # byte arrays
+            if isinstance(obj, (bytearray, memoryview)):
+                return bytes(obj)
+            return tuple(_freeze(v, stack) for v in obj)
+
+        try:
+            hash(obj)
+            return obj
+        except TypeError as e:
+            raise TypeError(f"Don't know how to freeze type '{type(obj)}'") from e
+    finally:
+        stack.discard(id(obj))
+
+
 class _MockModule:
     def __init__(self, module: ModuleType, rewrites: t.Dict[str, t.Callable], wrap: t.Callable):
         self._inner: ModuleType = module
@@ -269,5 +315,5 @@ def inner(*args, **kwargs):
     'create_rng', 'create_rng_group',
     'create_sparse_groupings', 'create_compact_groupings',
     'mask_fraction_of_groups', 'FloatKey',
-    'unwrap',
+    'unwrap', 'freeze',
 ]
diff --git a/phaser/utils/scan.py b/phaser/utils/scan.py
index b559f1d..24f173e 100644
--- a/phaser/utils/scan.py
+++ b/phaser/utils/scan.py
@@ -7,21 +7,21 @@
 import numpy
 from numpy.typing import ArrayLike, DTypeLike, NDArray
 
-from .num import get_array_module, cast_array_module, NumT
+from .num import FloatT, cast_array_module, get_array_module
 
 
 @t.overload
-def make_raster_scan(shape: t.Tuple[int, int], scan_step: ArrayLike,  # pyright: ignore[reportOverlappingOverload]
-                     rotation: float = 0., affine: t.Union[None, ArrayLike] = None, *, dtype: NumT, xp: t.Any = None) -> NDArray[NumT]:
+def make_raster_scan(shape: tuple[int, int], scan_step: ArrayLike,
+                     rotation: float = 0., affine: None | ArrayLike = None, *, dtype: FloatT, xp: t.Any = None) -> NDArray[FloatT]:
     ...
 
 @t.overload
-def make_raster_scan(shape: t.Tuple[int, int], scan_step: ArrayLike,
-                     rotation: float = 0., affine: t.Union[None, ArrayLike] = None, *, dtype: t.Optional[DTypeLike] = None, xp: t.Any = None) -> NDArray[numpy.floating]:
+def make_raster_scan(shape: tuple[int, int], scan_step: ArrayLike,
+                     rotation: float = 0., affine: None | ArrayLike = None, *, dtype: DTypeLike | None = None, xp: t.Any = None) -> NDArray[numpy.floating]:
     ...
 
-def make_raster_scan(shape: t.Tuple[int, int], scan_step: ArrayLike,
-                     rotation: float = 0., affine: t.Union[None, ArrayLike] = None, *, dtype: t.Any = None, xp: t.Any = None) -> NDArray[numpy.number]:
+def make_raster_scan(shape: tuple[int, int], scan_step: ArrayLike,
+                     rotation: float = 0., affine: None | ArrayLike = None, *, dtype: t.Any = None, xp: t.Any = None) -> NDArray[numpy.floating]:
     """
     Make a raster scan, centered around the origin.
 
@@ -56,7 +56,7 @@ def make_raster_scan(shape: t.Tuple[int, int], scan_step: ArrayLike,
         mat = xp2.asarray([[numpy.cos(theta), -numpy.sin(theta)], [numpy.sin(theta), numpy.cos(theta)]], dtype=dtype)
         pts = (pts @ mat.T)
 
-    return t.cast(NDArray[numpy.number], pts)
+    return t.cast(NDArray[numpy.floating], pts)
 
 
 __all__ = [
diff --git a/phaser/web/views.py b/phaser/web/views.py
index 28e0989..7e254a2 100644
--- a/phaser/web/views.py
+++ b/phaser/web/views.py
@@ -34,9 +34,11 @@ def _probes_view(cache: Cache, params: t.Mapping[str, t.Any]) -> t.Any:
 
 
 def _positions_view(cache: Cache, params: t.Mapping[str, t.Any]) -> t.Any:
-    # `ReconsState.scan`, wire-form as sent (shape (..., 2), in length units). Passed
+    # `ScanState.data`, wire-form as sent (shape (..., 2), in length units). Passed
     # through verbatim like `_probes_view`; the client flattens the leading axes.
-    return cache.raw.get('scan')
+    # `initial` and `tilt` ride along in the same payload, unused by this view.
+    scan = cache.raw.get('scan')
+    return scan['data'] if scan is not None else None
 
 
 def _recip_probes(cache: Cache) -> t.Any:
diff --git a/pyproject.toml b/pyproject.toml
index 5098e8d..2fd681f 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -34,10 +34,11 @@ dependencies = [
     "h5py~=3.8",
     "pyyaml>=5.3.1",
     "click~=8.1.0",
+    "frozendict~=2.4.7",
     "rich>=12.0.0,<15",
     "tifffile>=2023.8.25",
     "optree>=0.13.0",
-    "py-pane==0.11.6",
+    "py-pane==0.11.7",
     "typing_extensions~=4.7",
     "platformdirs>=4.8.0",
     "PyYAML>=6.0",
diff --git a/tests/input_files/state_v0.1.h5 b/tests/input_files/state_v0.1.h5
new file mode 100644
index 0000000..3f6f44b
--- /dev/null
+++ b/tests/input_files/state_v0.1.h5
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ec8b5b121f3d063f624634fcda6c0b7d4aa43965fa89efa8cf98e6ef0b769a80
+size 23928
diff --git a/tests/input_files/state_v0.1_no_tilt.h5 b/tests/input_files/state_v0.1_no_tilt.h5
new file mode 100644
index 0000000..f872206
--- /dev/null
+++ b/tests/input_files/state_v0.1_no_tilt.h5
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:964ef83e2db64fce238a0b205ee39b579e16fc4da9df883db3becc786dcfe6d7
+size 23152
diff --git a/tests/test_engines.py b/tests/test_engines.py
new file mode 100644
index 0000000..f226aba
--- /dev/null
+++ b/tests/test_engines.py
@@ -0,0 +1,82 @@
+import numpy
+import pytest
+from numpy.testing import assert_array_equal
+
+from phaser.execute import execute_engine, initialize_reconstruction
+from phaser.hooks import RawData
+from phaser.plan import ReconsPlan
+from phaser.state import ScanState
+from phaser.utils.num import Sampling, get_backend_module
+
+from .utils import make_recons_state
+
+SCAN_SHAPE = (8, 8)
+DET_SHAPE = (32, 32)
+
+
+def load_random(args, props) -> RawData:
+    rng = numpy.random.default_rng(0x9E3D)
+
+    return {
+        'patterns': rng.random((*SCAN_SHAPE, *DET_SHAPE), dtype=numpy.float32) + 1e-3,
+        'mask': numpy.ones(DET_SHAPE, dtype=numpy.float32),
+        'sampling': Sampling(DET_SHAPE, sampling=(1.0, 1.0)),
+        'wavelength': 0.0251,
+        'scan_hook': {'type': 'raster', 'shape': SCAN_SHAPE, 'step_size': (0.6, 0.6)},
+        'probe_hook': {'type': 'focused', 'conv_angle': 25.0, 'defocus': 0.0},
+        'seed': None,
+    }
+
+
+def test_position_update_preserves_initial():
+    plan = ReconsPlan.from_data({
+        'name': 'test',
+        'backend': 'numpy',
+        'raw_data': 'tests.test_engines:load_random',
+        'engines': [{
+            'type': 'conventional',
+            'probe_modes': 1,
+            'niter': 2,
+            'grouping': 16,
+            'noise_model': {'type': 'amplitude'},
+            'solver': {'type': 'lsqml'},
+            'position_solver': {'type': 'momentum', 'step_size': 1e-2},
+            'update_positions': True,
+            'iter_constraints': [],
+            'group_constraints': [],
+        }],
+    })
+
+    recons = initialize_reconstruction(plan)
+    initial = recons.state.scan.initial.copy()
+
+    for engine in plan.engines:
+        recons = execute_engine(recons, engine)
+
+    scan = recons.state.scan
+    assert not numpy.allclose(scan.data, initial)
+    assert_array_equal(scan.initial, initial)
+
+
+@pytest.mark.jax
+def test_gradient_group_indexing():
+    try:
+        get_backend_module('jax')
+    except ValueError as e:
+        pytest.skip(str(e))
+
+    from phaser.engines.gradient.run import extract_vars, insert_vars
+
+    flat = numpy.arange(32.).reshape(16, 2)
+    scan = ScanState(flat, flat + 100., flat * 0.01)
+    group = numpy.array([[1, 3, 5, 7]])
+
+    (vars, stripped) = extract_vars(make_recons_state(scan), {'positions'}, group)
+    assert_array_equal(vars['positions'], flat[[1, 3, 5, 7]])
+
+    grouped = insert_vars(vars, stripped, group)
+
+    # `initial` and `tilt` are indexed alongside `data`
+    assert_array_equal(grouped.scan.data, flat[[1, 3, 5, 7]])
+    assert_array_equal(grouped.scan.initial, flat[[1, 3, 5, 7]] + 100.)
+    assert_array_equal(grouped.scan.tilt, flat[[1, 3, 5, 7]] * 0.01)
diff --git a/tests/test_initialization.py b/tests/test_initialization.py
index b8494b4..a931e9a 100644
--- a/tests/test_initialization.py
+++ b/tests/test_initialization.py
@@ -1,17 +1,29 @@
 # type: ignore
 
-import re
 import logging
+import re
+import typing as t
 
 import numpy
-import pane
 import pytest
+from numpy.testing import assert_allclose, assert_array_equal
+from numpy.typing import NDArray
 
-from phaser.utils.num import Sampling
-from phaser.hooks import RawData
+import pane
+from phaser.execute import (
+    _normalize_scan_shape,
+    initialize_reconstruction,
+    load_raw_data,
+)
+from phaser.hooks import DropNanProps, RasterScanProps, RawData
+from phaser.hooks.preprocessing import drop_nan_patterns
+from phaser.hooks.scan import raster_scan
 from phaser.plan import ReconsPlan
-from phaser.execute import load_raw_data, initialize_reconstruction
-from phaser.state import PartialReconsState, ProbeState
+from phaser.state import PartialReconsState, Patterns, ProbeState, ScanState
+from phaser.utils.misc import freeze
+from phaser.utils.num import Sampling
+
+from .utils import make_recons_state
 
 
 def load_empty(args, props) -> RawData:
@@ -110,7 +122,7 @@ def test_load_raw_data_prev_state(caplog):
     }
 
     probe_state = ProbeState(Sampling((64, 64), sampling=(1.0, 1.0)), numpy.zeros((64, 64), dtype=numpy.complex64))
-    scan_state = numpy.zeros((32, 32, 2))
+    scan_state = ScanState(numpy.zeros((32, 32, 2)), numpy.zeros((32, 32, 2)))
 
     xp = numpy
     with caplog.at_level(logging.WARNING):
@@ -133,7 +145,7 @@ def test_load_raw_data_prev_state(caplog):
     # probe from state overrides probe from raw data
     assert numpy.all(numpy.isclose(recons.state.probe.data, probe_state.data))
     # but scan should be modeled
-    assert ~numpy.all(numpy.isclose(recons.state.scan, scan_state))
+    assert ~numpy.all(numpy.isclose(recons.state.scan.data, scan_state.data))
 
     plan['init'] = {
         'scan': {},
@@ -150,7 +162,7 @@ def test_load_raw_data_prev_state(caplog):
 
     # both should be modeled
     assert ~numpy.all(numpy.isclose(recons.state.probe.data, probe_state.data))
-    assert ~numpy.all(numpy.isclose(recons.state.scan, scan_state))
+    assert ~numpy.all(numpy.isclose(recons.state.scan.data, scan_state.data))
 
 
 def test_load_3d_raw_data():
@@ -180,5 +192,200 @@ def test_load_3d_raw_data():
     })
     recons = initialize_reconstruction(plan)
 
-    assert recons.state.scan.shape == (*scan_shape, 2)
-    assert recons.patterns.patterns.shape == (*scan_shape, *det_shape)
\ No newline at end of file
+    assert recons.state.scan.data.shape == (*scan_shape, 2)
+    assert recons.patterns.patterns.shape == (*scan_shape, *det_shape)
+
+def _plan_with_scan(**init: object) -> ReconsPlan:
+    return ReconsPlan.from_data({
+        'name': 'test',
+        'raw_data': {
+            'type': 'tests.test_initialization:load_empty',
+            'scan_shape': (8, 8),
+            'det_shape': (32, 32),
+        },
+        'engines': [],
+        'init': init,
+    })
+
+
+def test_initialize_scan_initial_independent():
+    plan = _plan_with_scan(
+        scan={'type': 'raster', 'shape': (8, 8), 'step_size': (1., 1.)},
+        probe={'type': 'focused', 'conv_angle': 20.0, 'defocus': 100.0},
+    )
+    scan = initialize_reconstruction(plan, xp=numpy).state.scan
+
+    assert not numpy.shares_memory(scan.data, scan.initial)
+
+    initial = scan.initial.copy()
+    scan.data += 1.
+    assert numpy.array_equal(scan.initial, initial)
+
+
+def test_initialize_reuses_scan_state():
+    plan = _plan_with_scan()
+
+    probe = ProbeState(
+        Sampling((32, 32), sampling=(1.0, 1.0)),
+        numpy.zeros((1, 32, 32), dtype=numpy.complex64),
+    )
+    # a previous state, as read back from disk
+    prev = ScanState(
+        numpy.arange(128.).reshape(8, 8, 2),
+        numpy.arange(128.).reshape(8, 8, 2) + 100.,
+        numpy.arange(128.).reshape(8, 8, 2) * 0.01,
+    )
+
+    scan = initialize_reconstruction(plan, xp=numpy, init_state=PartialReconsState(
+        wavelength=1.0, probe=probe, scan=prev,
+    )).state.scan
+
+    # `initial` is carried over, not reset to `data`
+    assert_allclose(scan.data, prev.data, rtol=1e-6)
+    assert_allclose(scan.initial, prev.initial, rtol=1e-6)
+    assert_allclose(scan.tilt, prev.tilt, rtol=1e-6)
+
+    assert scan.data.dtype == scan.initial.dtype == scan.tilt.dtype == numpy.float32
+
+    # the caller's state is left untouched
+    assert prev.data.dtype == numpy.float64
+    assert not numpy.shares_memory(scan.data, prev.data)
+    assert not numpy.shares_memory(scan.initial, prev.initial)
+
+
+@pytest.mark.parametrize('flat_scan', (True, False))
+def test_normalize_scan_shape_keeps_initial(flat_scan: bool):
+    flat = numpy.arange(128.).reshape(64, 2)
+    shape = (64, 2) if flat_scan else (8, 8, 2)
+    scan = ScanState(flat.reshape(shape), (flat + 100.).reshape(shape), (flat * 0.01).reshape(shape))
+
+    patterns_shape = (8, 8, 4, 4) if flat_scan else (64, 4, 4)
+    patterns = Patterns(
+        numpy.zeros(patterns_shape, dtype=numpy.float32),
+        numpy.ones((4, 4), dtype=numpy.float32),
+    )
+
+    (patterns, state) = _normalize_scan_shape(patterns, make_recons_state(scan))
+
+    assert patterns.patterns.shape == (8, 8, 4, 4)
+    assert state.scan.data.shape == (8, 8, 2)
+    assert state.scan.initial.shape == (8, 8, 2)
+    assert state.scan.tilt.shape == (8, 8, 2)
+
+    assert numpy.array_equal(state.scan.initial, (flat + 100.).reshape(8, 8, 2))
+
+
+def test_drop_nan_patterns_filters_initial():
+    patterns = numpy.zeros((4, 4, 2, 2), dtype=numpy.float32)
+    patterns[0, 0] = numpy.nan
+    patterns[2, 3] = numpy.nan
+
+    flat = numpy.arange(32.).reshape(16, 2)
+    scan = ScanState(flat.reshape(4, 4, 2), flat.reshape(4, 4, 2) + 100., flat.reshape(4, 4, 2) * 0.01)
+
+    (data, state) = drop_nan_patterns({
+        'data': Patterns(patterns, numpy.ones((2, 2), dtype=numpy.float32)),
+        'state': make_recons_state(scan),
+        'seed': None, 'dtype': numpy.float32, 'xp': numpy,
+    }, DropNanProps(threshold=0.5))
+
+    kept = numpy.ones(16, dtype=numpy.bool_)
+    kept[[0, 11]] = False
+
+    assert data.patterns.shape == (14, 2, 2)
+    assert numpy.array_equal(state.scan.data, flat[kept])
+    assert numpy.array_equal(state.scan.initial, flat[kept] + 100.)
+    assert numpy.array_equal(state.scan.tilt, flat[kept] * 0.01)
+
+
+def _make_raster(shape: tuple[int, int], step: float = 1.) -> ScanState:
+    return raster_scan(
+        {'seed': None, 'dtype': numpy.float64, 'xp': numpy},
+        RasterScanProps(shape=shape, step_size=(step, step)),
+    )
+
+
+def _grid_indices(shape: tuple[int, int]) -> NDArray[numpy.int64]:
+    """Row and column index of every point of a `shape` grid, as a (..., 2) array."""
+    return numpy.stack(numpy.indices(shape), axis=-1)
+
+
+def _raster_grid(meta: t.Mapping[str, t.Any]) -> NDArray[numpy.int64]:
+    """`raster_rows` and `raster_cols` stacked into a (..., 2) index array."""
+    return numpy.stack((
+        numpy.array(meta['raster_rows']), numpy.array(meta['raster_cols']),
+    ), axis=-1)
+
+
+def _raster_positions(grid: NDArray[numpy.int64], shape: tuple[int, int], step: float) -> NDArray[numpy.float64]:
+    """Positions `make_raster_scan` assigns to the given grid indices."""
+    return (grid - numpy.array(shape) / 2.) * step
+
+
+def test_raster_scan_meta():
+    # non-square, so a transposed grid can't pass
+    shape = (3, 5)
+    scan = _make_raster(shape)
+
+    assert scan.meta['type'] == 'raster'
+
+    grid = _raster_grid(scan.meta)
+    assert grid.shape == (*shape, 2)
+    assert_array_equal(grid, _grid_indices(shape))
+
+    # each index labels the grid position of the matching scan point
+    assert_allclose(scan.data, _raster_positions(grid, shape, 1.))
+
+    # `meta` is a static pytree field, so it must be hashable
+    hash(scan.meta)
+
+
+@pytest.mark.parametrize('flat_meta', (True, False))
+def test_normalize_scan_shape_keeps_raster_meta(flat_meta: bool):
+    shape = (4, 8)
+    scan = _make_raster(shape)
+
+    if flat_meta:
+        # a state resumed after `drop_nan_patterns`, which flattens scan and metadata alike
+        flat = {k: numpy.array(v).ravel() if k.startswith('raster_') else v for (k, v) in scan.meta.items()}
+        scan = ScanState(scan.data.reshape(-1, 2), scan.initial.reshape(-1, 2), meta=freeze(flat))
+
+    patterns = Patterns(
+        numpy.zeros((*shape, 4, 4), dtype=numpy.float32),
+        numpy.ones((4, 4), dtype=numpy.float32),
+    )
+
+    (_, state) = _normalize_scan_shape(patterns, make_recons_state(scan))
+
+    grid = _raster_grid(state.scan.meta)
+    assert grid.shape == (*shape, 2)
+    assert_array_equal(grid, _grid_indices(shape))
+    assert_allclose(state.scan.data, _raster_positions(grid, shape, 1.))
+
+    hash(state.scan.meta)
+
+
+def test_drop_nan_patterns_keeps_raster_meta():
+    shape = (4, 4)
+    scan = _make_raster(shape)
+
+    patterns = numpy.zeros((*shape, 2, 2), dtype=numpy.float32)
+    patterns[0, 0] = numpy.nan
+    patterns[2, 3] = numpy.nan
+
+    (_, state) = drop_nan_patterns({
+        'data': Patterns(patterns, numpy.ones((2, 2), dtype=numpy.float32)),
+        'state': make_recons_state(scan),
+        'seed': None, 'dtype': numpy.float32, 'xp': numpy,
+    }, DropNanProps(threshold=0.5))
+
+    kept = numpy.ones(16, dtype=numpy.bool_)
+    kept[[0, 11]] = False
+
+    # dropped positions take their indices with them, the rest keep theirs
+    grid = _raster_grid(state.scan.meta)
+    assert grid.shape == (14, 2)
+    assert_array_equal(grid, _grid_indices(shape).reshape(-1, 2)[kept])
+    assert_allclose(state.scan.data, _raster_positions(grid, shape, 1.))
+
+    hash(state.scan.meta)
diff --git a/tests/test_load.py b/tests/test_load.py
index 325cfbc..3b45270 100644
--- a/tests/test_load.py
+++ b/tests/test_load.py
@@ -2,15 +2,15 @@
 import typing as t
 
 import numpy
-from numpy.typing import NDArray
-from numpy.testing import assert_allclose
-import pane
 import pytest
-from phaser.utils.num import get_backend_module
 import tifffile
+from numpy.testing import assert_allclose
+from numpy.typing import NDArray
+
+import pane
+from phaser.hooks import RawData, RawDataHook
+from phaser.utils.num import Sampling, get_backend_module
 
-from phaser.utils.num import Sampling
-from phaser.hooks import RawDataHook, RawData
 from .utils import EXPECTED_PATH, INPUT_FILES_PATH
 
 
diff --git a/tests/test_state.py b/tests/test_state.py
new file mode 100644
index 0000000..2a8cbab
--- /dev/null
+++ b/tests/test_state.py
@@ -0,0 +1,255 @@
+import re
+from pathlib import Path
+
+import h5py
+import numpy
+import pytest
+from frozendict import frozendict
+from numpy.testing import assert_allclose
+
+from phaser.state import (
+    IterState,
+    ObjectState,
+    PartialReconsState,
+    ProbeState,
+    ProgressState,
+    ReconsState,
+    ScanState,
+)
+from phaser.utils.io import hdf5_read_state
+from phaser.utils.num import Sampling, get_backend_module
+from phaser.utils.object import ObjectSampling
+
+from .utils import INPUT_FILES_PATH
+
+# scan and tilt of the checked-in v0.1 fixtures
+V0_1_SCAN = numpy.stack(numpy.meshgrid(
+    numpy.arange(4.) * 0.6, numpy.arange(4.) * 0.6, indexing='ij',
+), axis=-1)
+V0_1_TILT = numpy.linspace(-1., 1., 32).reshape(4, 4, 2)
+
+
+def make_state(*, region: bool = True, meta: bool = True) -> ReconsState:
+    rng = numpy.random.default_rng(0x5CA4)
+
+    probe = ProbeState(
+        Sampling((8, 8), sampling=(0.5, 0.5)),
+        (rng.random((2, 8, 8)) + 1.j * rng.random((2, 8, 8))).astype(numpy.complex64),
+        meta=frozendict({'source': 'test'}) if meta else frozendict(),
+    )
+    obj = ObjectState(
+        ObjectSampling(
+            (16, 16), (0.5, 0.5), (-4., -4.),
+            (-2., -2.) if region else None, (2., 2.) if region else None,
+        ),
+        (rng.random((2, 16, 16)) + 1.j * rng.random((2, 16, 16))).astype(numpy.complex64),
+        numpy.array([10., 10.]),
+        meta=frozendict({'slices': 2}) if meta else frozendict(),
+    )
+    # data, initial and tilt are all distinct
+    scan = ScanState(
+        rng.random((4, 4, 2)) * 2.,
+        rng.random((4, 4, 2)) * 2.,
+        rng.random((4, 4, 2)) * 0.1,
+        meta=frozendict({'step_size': [0.6, 0.6]}) if meta else frozendict(),
+    )
+
+    return ReconsState(
+        iter=IterState(2, 5, 15), wavelength=0.0251,
+        probe=probe, object=obj, scan=scan,
+        progress={'total_loss': ProgressState([1, 2], [3., 2.])},
+    )
+
+
+@pytest.mark.parametrize('region', (True, False))
+def test_state_hdf5_roundtrip(tmp_path: Path, region: bool):
+    state = make_state(region=region)
+    path = tmp_path / 'state.h5'
+    state.write_hdf5(path)
+
+    with h5py.File(path) as f:
+        assert f['probe']['type'][()] == b'pixelated'
+        assert f['object']['type'][()] == b'pixelated'
+
+    read = ReconsState.read_hdf5(path)
+
+    assert read.wavelength == state.wavelength
+    assert read.probe.ty == read.object.ty == 'pixelated'
+    assert (read.iter.engine_num, read.iter.engine_iter, read.iter.total_iter) == (2, 5, 15)
+    assert read.progress['total_loss'].iters == [1, 2]
+    assert read.progress['total_loss'].values == [3., 2.]
+
+    assert read.probe.sampling == state.probe.sampling
+    assert_allclose(read.probe.data, state.probe.data)
+
+    assert read.object.sampling == state.object.sampling
+    assert_allclose(read.object.data, state.object.data)
+    assert_allclose(read.object.thicknesses, state.object.thicknesses)
+
+    if region:
+        assert_allclose(read.object.sampling.region_min, [-2., -2.])
+        assert_allclose(read.object.sampling.region_max, [2., 2.])
+    else:
+        assert read.object.sampling.region_min is None
+        assert read.object.sampling.region_max is None
+
+    assert_allclose(read.scan.data, state.scan.data)
+    assert_allclose(read.scan.initial, state.scan.initial)
+    assert_allclose(read.scan.tilt, state.scan.tilt)
+
+    # `initial` is stored, not synthesized from `data`
+    assert not numpy.allclose(read.scan.data, read.scan.initial)
+    assert not numpy.shares_memory(read.scan.data, read.scan.initial)
+
+
+def test_state_hdf5_roundtrip_no_tilt(tmp_path: Path):
+    state = make_state()
+    state.scan.tilt = None
+    path = tmp_path / 'state.h5'
+    state.write_hdf5(path)
+
+    with h5py.File(path) as f:
+        assert 'tilt' not in f['scan']
+
+    assert ReconsState.read_hdf5(path).scan.tilt is None
+
+
+def test_state_hdf5_meta_roundtrip(tmp_path: Path):
+    meta = frozendict({
+        'str': 'value', 'int': 3, 'float': 1.5, 'bool': True, 'null': None,
+        'list': [1, [2, 3]], 'dict': {'inner': [4]},
+    })
+    # JSON arrays are frozen to tuples, making `meta` hashable
+    expected = frozendict({
+        'str': 'value', 'int': 3, 'float': 1.5, 'bool': True, 'null': None,
+        'list': (1, (2, 3)), 'dict': frozendict({'inner': (4,)}),
+    })
+
+    state = make_state()
+    state.scan.meta = meta
+    path = tmp_path / 'state.h5'
+    state.write_hdf5(path)
+    read = ReconsState.read_hdf5(path)
+
+    assert read.scan.meta == expected
+    hash(read.scan.meta)
+
+    assert read.probe.meta == state.probe.meta
+    assert read.object.meta == state.object.meta
+
+
+def test_state_hdf5_meta_empty(tmp_path: Path):
+    state = make_state(meta=False)
+    path = tmp_path / 'state.h5'
+    state.write_hdf5(path)
+
+    # empty metadata leaves no dataset behind
+    with h5py.File(path) as f:
+        assert 'meta' not in f['scan']
+        assert 'meta' not in f['probe']
+        assert 'meta' not in f['object']
+
+    assert ReconsState.read_hdf5(path).scan.meta == frozendict()
+
+
+@pytest.mark.jax
+def test_state_pytree_leaves_are_numeric():
+    # `tree` dispatches to jax directly, so the backend must be loaded before
+    # the state classes are registered as pytree nodes
+    try:
+        get_backend_module('jax')
+    except ValueError as e:
+        pytest.skip(str(e))
+
+    from phaser.utils import tree
+
+    # every non-static field is an array; discriminators and metadata are static
+    for (path, leaf) in tree.leaves_with_path(make_state()):
+        assert numpy.issubdtype(numpy.asarray(leaf).dtype, numpy.number), \
+            f"non-numeric pytree leaf at '{''.join(map(str, path))}'"
+
+
+def test_state_meta_survives_backend_roundtrip():
+    state = make_state()
+    round_tripped = state.to_xp(numpy).to_numpy()
+
+    assert round_tripped.probe.meta == state.probe.meta
+    assert round_tripped.object.meta == state.object.meta
+    assert round_tripped.scan.meta == state.scan.meta
+
+    assert state.probe.resample(state.probe.sampling).meta == state.probe.meta
+
+
+@pytest.mark.parametrize(('name', 'has_tilt'), (
+    ('state_v0.1.h5', True),
+    ('state_v0.1_no_tilt.h5', False),
+))
+def test_read_state_v0_1(name: str, has_tilt: bool):
+    """Fixtures written by the pre-`ScanState` writer, at commit 55d531e."""
+    state = hdf5_read_state(INPUT_FILES_PATH / name)
+    assert state.scan is not None
+
+    assert_allclose(state.scan.data, V0_1_SCAN)
+    # v0.1 has no `initial`, so the stored scan stands in for it
+    assert_allclose(state.scan.initial, V0_1_SCAN)
+    assert not numpy.shares_memory(state.scan.data, state.scan.initial)
+
+    if has_tilt:
+        assert_allclose(state.scan.tilt, V0_1_TILT)
+    else:
+        assert state.scan.tilt is None
+
+    assert state.scan.meta == frozendict()
+    assert state.probe is not None and state.probe.meta == frozendict()
+    # v0.1 has no 'type', which reads back as pixelated
+    assert state.probe.ty == 'pixelated'
+
+    # v0.1 wrote absent regions as empty datasets
+    assert state.object is not None
+    assert state.object.ty == 'pixelated'
+    assert state.object.sampling.region_min is None
+    assert state.object.sampling.region_max is None
+
+
+@pytest.mark.parametrize(('group', 'ty'), (
+    ('probe', 'fake_type'),
+    ('object', 'fake_type'),
+))
+def test_read_state_unsupported_type(tmp_path: Path, group: str, ty: str):
+    path = tmp_path / 'state.h5'
+    make_state().write_hdf5(path)
+
+    with h5py.File(path, 'r+') as f:
+        del f[group]['type']
+        f[group].create_dataset('type', data=ty, dtype=h5py.string_dtype('utf-8'))
+
+    with pytest.raises(ValueError, match=re.escape(f"Unsupported {group} type '{ty}'")):
+        hdf5_read_state(path)
+
+
+def test_read_state_unsupported_version(tmp_path: Path):
+    path = tmp_path / 'state.h5'
+    with h5py.File(path, 'w') as f:
+        f.create_dataset('type', (), h5py.string_dtype(), "phaser_state")
+        f.create_dataset('version', (), h5py.string_dtype(), "0.3")
+
+    with pytest.raises(ValueError, match=re.escape("Unsupported file version '0.3'")):
+        hdf5_read_state(path)
+
+
+def test_partial_state_hdf5_roundtrip(tmp_path: Path):
+    state = make_state()
+    path = tmp_path / 'state.h5'
+    PartialReconsState(wavelength=state.wavelength, probe=state.probe).write_hdf5(path)
+    read = PartialReconsState.read_hdf5(path)
+
+    assert read.scan is None
+    assert read.object is None
+    assert read.progress is None
+    assert read.iter == IterState.empty()
+
+    assert read.probe is not None
+    assert_allclose(read.probe.data, state.probe.data)
+
+    with pytest.raises(ValueError, match=re.escape("ReconsState missing 'object', 'scan'")):
+        read.to_complete()
diff --git a/tests/utils.py b/tests/utils.py
index 2cd7bb8..0638b65 100644
--- a/tests/utils.py
+++ b/tests/utils.py
@@ -1,14 +1,17 @@
 import builtins
-from contextlib import contextmanager
-from itertools import chain
 import inspect
-from pathlib import Path
 import sys
 import typing as t
+from contextlib import contextmanager
+from itertools import chain
+from pathlib import Path
 
 import numpy
 import pytest
 
+if t.TYPE_CHECKING:
+    from phaser.state import ReconsState, ScanState
+
 CallableT = t.TypeVar('CallableT', bound=t.Callable)
 P = t.ParamSpec('P')
 T = t.TypeVar('T')
@@ -20,7 +23,7 @@
 
 
 def _wrap_pytest(wrapper: CallableT, wrapped: t.Callable,
-                 mod_params: t.Optional[t.Callable[[t.Sequence[inspect.Parameter]], t.Sequence[inspect.Parameter]]] = None
+                 mod_params: t.Callable[[t.Sequence[inspect.Parameter]], t.Sequence[inspect.Parameter]] | None = None
 ) -> CallableT:
     # hacks to allow pytest to find fixtures in wrapped functions
     old_sig = inspect.signature(wrapped)
@@ -39,10 +42,10 @@ def _wrap_pytest(wrapper: CallableT, wrapped: t.Callable,
 
 
 def with_backends(
-    *backends: t.Union[str, t.Iterable[str]]
+    *backends: str | t.Iterable[str]
 ) -> t.Callable[[t.Callable[P, T]], t.Callable[P, T]]:
     """Run a test on the specified compute backends"""
-    backends = t.cast(t.Tuple[str, ...],
+    backends = t.cast(tuple[str, ...],
         tuple(chain.from_iterable((b,) if isinstance(b, str) else b for b in backends))
     )
 
@@ -103,7 +106,7 @@ def write_array(path: Path, arr: numpy.ndarray):
         raise RuntimeError(f"Unable to save file '{path.name}'") from e
 
 
-def check_array_equals_file(name: str, *, out_name: t.Optional[str] = None, decimal: int = 6) -> t.Callable[[t.Callable[..., numpy.ndarray]], t.Callable[..., None]]:
+def check_array_equals_file(name: str, *, out_name: str | None = None, decimal: int = 6) -> t.Callable[[t.Callable[..., numpy.ndarray]], t.Callable[..., None]]:
     def decorator(f: t.Callable[..., numpy.ndarray]):
         @pytest.mark.expected_filename(name)
         def wrapper(*args, file_contents_array: numpy.ndarray, **kwargs):
@@ -128,8 +131,8 @@ def wrapper(*args, file_contents_array: numpy.ndarray, **kwargs):
                 try:
                     print(f"Saving actual result to '{out_path}'")
                     write_array(out_path, actual)
-                except Exception:
-                    print("Failed to save result.")
+                except Exception as e:  # noqa: BLE001
+                    print(f"Failed to save result: {e}")
                 raise
 
         return _wrap_pytest(wrapper, f, lambda params: [*params, inspect.Parameter('file_contents_array', inspect.Parameter.KEYWORD_ONLY)])
@@ -177,4 +180,25 @@ def mock_importerror(*modulenames: t.Iterable[str]):
 
     # if we're the last ones, change __import__ back
     if not len(builtins.__import__.prevented):
-        builtins.__import__ = _import
\ No newline at end of file
+        builtins.__import__ = _import
+
+
+def make_recons_state(scan: 'ScanState') -> 'ReconsState':
+    """Minimal `ReconsState` wrapping `scan`, for tests which only exercise the scan."""
+    from phaser.state import IterState, ObjectState, ProbeState, ReconsState
+    from phaser.utils.num import Sampling
+    from phaser.utils.object import ObjectSampling
+
+    return ReconsState(
+        iter=IterState.empty(), wavelength=1.0,
+        probe=ProbeState(
+            Sampling((4, 4), sampling=(1., 1.)),
+            numpy.zeros((1, 4, 4), dtype=numpy.complex64),
+        ),
+        object=ObjectState(
+            ObjectSampling((4, 4), (1., 1.)),
+            numpy.ones((1, 4, 4), dtype=numpy.complex64),
+            numpy.array([]),
+        ),
+        scan=scan,
+    )
diff --git a/tests/web/test_web_views.py b/tests/web/test_web_views.py
index 62bf305..121b252 100644
--- a/tests/web/test_web_views.py
+++ b/tests/web/test_web_views.py
@@ -217,17 +217,18 @@ def counting_decode(obj):
 
 def test_progress_probes_are_raw_passthrough():
     cache = Cache()
-    scan = encode_obj(numpy.zeros((4, 5, 2), dtype=numpy.float32))
+    positions = encode_obj(numpy.zeros((4, 5, 2), dtype=numpy.float32))
     cache.update_raw({
         'progress': {'total_loss': {'iters': [1], 'values': [0.5]}},
-        'probe': {'sampling': {}, 'data': 'x'}, 'scan': scan,
+        'probe': {'sampling': {}, 'data': 'x'},
+        'scan': {'data': positions, 'initial': 'i', 'tilt': 'ti', 'meta': {}},
     })
 
     assert VIEWS['progress'].compute(cache, {}) == {'total_loss': {'iters': [1], 'values': [0.5]}}
     # the bulk array itself, not the wrapper -- and the same object, not a copy
     assert VIEWS['probes'].compute(cache, {}) is cache.raw['probe']['data']
-    # positions pass through in whatever (..., 2) shape the worker sent; the client flattens
-    assert VIEWS['positions'].compute(cache, {}) is scan
+    # `ScanState.data` only, in whatever (..., 2) shape the worker sent; the client flattens
+    assert VIEWS['positions'].compute(cache, {}) is positions
     assert VIEWS['positions'].deps == frozenset({'scan'})