diff --git a/CHANGELOG.md b/CHANGELOG.md
index d338f76f..61d66709 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,28 @@
 
 ## Unreleased
 
+### Added
+
+- PyTorch cluster-tile selective calls can append caller-owned tile state with
+  `return_state=True`, without changing the default neighbor-list return arity.
+
+### Fixed
+
+- Unbatched JAX naive dual-cutoff PBC neighbor lists now populate both cutoff
+  outputs when using the default `wrap_positions=True`. Previously this path
+  wrapped positions but skipped the fill kernel, leaving zero counts and padded
+  matrices.
+- Batched PyTorch cluster-tile segmented COO validates fixed topology, offsets,
+  counts, and tile-state capacities before launching Warp kernels.
+- Single-system Torch and JAX segmented cluster-tile COO now require one exact
+  physical interval, bound writes by output capacity, fail closed for malformed
+  offsets, and cap compiled/JIT active counts to writable capacity. Batched
+  per-system physical subsegments remain supported.
+- Compiled unified PyTorch cluster-tile dispatch now rejects tensor-valued PBC
+  rather than treating it as fully periodic. Eagerly validate PBC and compile the
+  direct single-system fixed-state route instead.
+- JAX cluster-tile empty selective rebuilds now preserve false-flag state and
+  clear true-flag pair and tile counts while retaining fixed-capacity storage.
 ### Changed (neighbors)
 
 - Improved JAX neighbor-list import performance by deferring dtype-specific
diff --git a/docs/userguide/components/neighborlist.md b/docs/userguide/components/neighborlist.md
index 4df751f9..b1353ba9 100644
--- a/docs/userguide/components/neighborlist.md
+++ b/docs/userguide/components/neighborlist.md
@@ -1245,6 +1245,140 @@ segmented-COO outputs in both the PyTorch and JAX `batch_naive` / `batch_cell_li
 paths (single-system paths take a whole-system flag of shape `(1,)`). It is not
 combined with differentiable per-pair geometry.
 
+Cluster-tile selective rebuilds also require the tile-list buffers. JAX returns that
+state as part of every selective result. Zero-atom calls keep the same arity as
+nonempty calls: selective matrix output contains the matrix triple plus tile state
+(and both matrix triples for dual cutoff), while segmented COO remains the normal
+seven-array single-system or ten-array batched tuple.
+
+PyTorch keeps the existing return arity by default. To receive the tile state, select
+the cluster-tile method explicitly and pass `return_state=True` with
+`rebuild_flags`:
+
+- Single-system matrix output appends `(num_tiles, tile_row_group,
+  tile_col_group)`, yielding six tensors for one cutoff and nine for two.
+- Single-system segmented COO returns `(neighbor_list, pair_offsets,
+  pair_counts, neighbor_list_shifts)` and appends `(num_tiles,
+  tile_row_group, tile_col_group)` with `return_state=True`, yielding seven
+  tensors. The fixed buffers have shapes `(2, max_pairs)`, `(2,)`, `(1,)`, and
+  `(max_pairs, 3)`; `pair_offsets` is `[0, max_pairs]`, and
+  `pair_counts[0]` gives the active prefix. Values after that prefix are
+  inactive.
+- Batched matrix output appends `(tile_offsets, tile_counts, num_tiles,
+  tile_row_group, tile_col_group, tile_system)`, yielding nine tensors for one
+  cutoff and twelve for two.
+- Batched segmented COO appends the same six state tensors to its four topology
+  outputs, yielding ten tensors.
+
+```python
+from nvalchemiops.torch.neighbors import neighbor_list
+
+(
+    neighbor_matrix,
+    num_neighbors,
+    shifts,
+    tile_offsets,
+    tile_counts,
+    num_tiles,
+    tile_row_group,
+    tile_col_group,
+    tile_system,
+) = neighbor_list(
+    positions,
+    cutoff,
+    cell=cells,
+    pbc=pbc,
+    batch_ptr=batch_ptr,
+    method="batch_cluster_tile",
+    rebuild_flags=rebuild_flags,
+    return_state=True,
+)
+```
+
+For a single-system fixed-capacity COO workflow, retain the public topology
+buffers and tile state from the previous step:
+
+```python
+(
+    neighbor_list_coo,
+    pair_offsets,
+    pair_counts,
+    shifts_coo,
+    num_tiles,
+    tile_row_group,
+    tile_col_group,
+) = neighbor_list(
+    positions,
+    cutoff,
+    cell=cell,
+    pbc=pbc,
+    method="cluster_tile",
+    return_neighbor_list=True,
+    rebuild_flags=rebuild_flags,  # shape (1,)
+    neighbor_list=neighbor_list_coo,
+    pair_offsets=pair_offsets,  # tensor([0, max_pairs], dtype=torch.int32)
+    pair_counts=pair_counts,
+    neighbor_list_shifts=shifts_coo,
+    num_tiles=num_tiles,
+    tile_row_group=tile_row_group,
+    tile_col_group=tile_col_group,
+    return_state=True,
+)
+```
+
+Pass the returned topology buffers and state back on the next call. `return_state`
+is intentionally rejected for auto-selected and non-cluster-tile methods, and it
+requires `rebuild_flags`. Caller-supplied state tensors are returned by identity:
+the suffix aliases the same buffers, which the next selective call updates in place.
+When only the public state suffix is supplied, temporary sorting scratch is allocated
+internally without replacing those state buffers.
+
+Use two phases for PyTorch selective COO. First make an eager bootstrap call with
+every flag true; it may allocate omitted topology and tile buffers and returns the
+reusable state. Every later eager call containing a false flag must supply all fixed
+topology and tile-state buffers, otherwise it raises `ValueError` before a kernel
+launch. Single-system Torch and JAX COO state has one segment:
+`pair_offsets` must stay `[0, physical_capacity]`, and `pair_counts` has one value.
+
+```python
+# Eager bootstrap: omit every persistent buffer.
+state = cluster_tile_neighbor_list(
+    positions, cutoff, cell, format="coo", return_state=True,
+    rebuild_flags=torch.ones(1, dtype=torch.bool, device=positions.device),
+)
+```
+
+The direct single-system `cluster_tile_neighbor_list` route also supports
+`torch.compile(fullgraph=True)` for this steady-state phase. Bootstrap and validate
+capacities eagerly, then compile a function that accepts fixed-shape buffers and
+`rebuild_flags`; compiled calls preserve false-flag data and rebuild true-flag data.
+False flags retain the fixed launch graph, so they do not promise zero preprocessing
+or kernel launches. The unified `neighbor_list(..., method="cluster_tile")` dispatch
+is unsupported under compilation because it cannot validate tensor-valued `pbc`
+without host synchronization. Compiled calls also omit
+host-synchronized offset and overflow diagnostics, so retain the eagerly validated
+capacities. The Warp COO query enforces physical output-buffer bounds as defense in
+depth, but mutated or malformed metadata remains unsupported. Compiled Torch and
+JIT-compiled JAX cannot synchronize to raise a data-dependent metadata error. If
+caller-provided single-system offsets no longer equal `[0, physical_capacity]`, a
+true rebuild is suppressed before it writes and the returned active count is zero.
+For valid offsets, an overflowed count is capped at writable capacity. A false
+rebuild retains valid prior buffers and counts; malformed metadata still returns a
+zero active count without changing those buffers. This prevents a returned count
+from naming unwritten entries; it does not make mutated metadata supported. Batched
+cluster-tile fullgraph support is not provided.
+
+```python
+@torch.compile(fullgraph=True)
+def reuse(flags, neighbor_list, pair_offsets, pair_counts, shifts, num_tiles, row, col):
+    return cluster_tile_neighbor_list(
+        positions, cutoff, cell, format="coo", return_state=True,
+        rebuild_flags=flags, neighbor_list=neighbor_list, pair_offsets=pair_offsets,
+        pair_counts=pair_counts, neighbor_list_shifts=shifts, num_tiles=num_tiles,
+        tile_row_group=row, tile_col_group=col,
+    )
+```
+
 ### Dual Cutoff
 
 Compute two neighbor lists with different cutoffs simultaneously:
diff --git a/nvalchemiops/jax/neighbors/batch_cluster_tile.py b/nvalchemiops/jax/neighbors/batch_cluster_tile.py
index b5fbcd4f..ed7c71f1 100644
--- a/nvalchemiops/jax/neighbors/batch_cluster_tile.py
+++ b/nvalchemiops/jax/neighbors/batch_cluster_tile.py
@@ -1335,7 +1335,9 @@ def batch_query_cluster_tile(
         with pair outputs.
     rebuild_flags : jax.Array, shape (S,), dtype=bool, optional
         Per-system selective rebuild flags. Requires ``tile_offsets``,
-        ``tile_counts``, and ``batch_idx``.
+        ``tile_counts``, and ``batch_idx``. For an empty batch, false-flag
+        segments retain their active pair and tile counts while true-flag
+        segments clear them; fixed-capacity arrays remain unchanged.
     tile_offsets : jax.Array, shape (S + 1,), dtype=int32, optional
         Per-system tile segment offsets. Required with ``rebuild_flags``.
     tile_counts : jax.Array, shape (S,), dtype=int32, optional
@@ -1954,7 +1956,10 @@ def batch_cluster_tile_neighbor_list(
     For ``format == "matrix"``:
         ``(neighbor_matrix, num_neighbors, neighbor_matrix_shifts)``, with
         optional ``(*, distances)`` and/or ``(*, vectors)`` appended when
-        ``return_distances`` / ``return_vectors`` is True.
+        ``return_distances`` / ``return_vectors`` is True. With
+        ``rebuild_flags``, append ``(tile_offsets, tile_counts, num_tiles,
+        tile_row_group, tile_col_group, tile_system)`` after one matrix triple,
+        or after both triples when ``cutoff2`` is provided.
     For ``format == "coo"``:
         ``(neighbor_list, neighbor_ptr, neighbor_list_shifts)`` in compact
         mode, or ``(neighbor_list, pair_offsets, pair_counts,
@@ -2087,6 +2092,29 @@ def batch_cluster_tile_neighbor_list(
         nn0 = jnp.empty(0, dtype=jnp.int32)
         ns0 = jnp.empty((0, int(max_neighbors), 3), dtype=jnp.int32)
         if format == "coo":
+            if selective:
+                empty_pair_counts = jnp.where(
+                    rebuild_flags,
+                    jnp.zeros_like(previous_pair_counts),
+                    previous_pair_counts,
+                )
+                empty_tile_counts = jnp.where(
+                    rebuild_flags,
+                    jnp.zeros_like(previous_tile_counts),
+                    previous_tile_counts,
+                )
+                return (
+                    previous_neighbor_list,
+                    pair_offsets,
+                    empty_pair_counts,
+                    previous_neighbor_list_shifts,
+                    tile_offsets,
+                    empty_tile_counts,
+                    previous_num_tiles,
+                    previous_tile_row_group,
+                    previous_tile_col_group,
+                    previous_tile_system,
+                )
             coo_base = (
                 jnp.empty((2, 0), dtype=jnp.int32),
                 jnp.zeros(1, dtype=jnp.int32),
@@ -2156,7 +2184,11 @@ def batch_cluster_tile_neighbor_list(
             return (
                 *matrix_out,
                 tile_offsets,
-                previous_tile_counts,
+                jnp.where(
+                    rebuild_flags,
+                    jnp.zeros_like(previous_tile_counts),
+                    previous_tile_counts,
+                ),
                 previous_num_tiles,
                 previous_tile_row_group,
                 previous_tile_col_group,
diff --git a/nvalchemiops/jax/neighbors/cluster_tile.py b/nvalchemiops/jax/neighbors/cluster_tile.py
index deea4ac3..6627c7ff 100644
--- a/nvalchemiops/jax/neighbors/cluster_tile.py
+++ b/nvalchemiops/jax/neighbors/cluster_tile.py
@@ -1079,8 +1079,10 @@ def query_cluster_tile(
     neighbor_matrix_shifts2 : jax.Array, shape (natom, max_neighbors, 3), dtype=int32, optional
         Second shift buffer for dual-cutoff mode.
     rebuild_flags : jax.Array, shape (1,), dtype=bool, optional
-        When set, skips atoms flagged as unchanged; requires all
-        ``previous_*`` buffers to be passed.
+        When set, requires all ``previous_*`` buffers. For an empty system a
+        false flag preserves the previous active counts and tile state, while
+        a true flag clears the active pair and tile counts without rewriting
+        fixed-capacity topology storage.
     return_vectors : bool, optional
         If True, also fill and return per-pair displacement vectors.
     return_distances : bool, optional
@@ -1396,6 +1398,74 @@ def query_cluster_tile(
     return neighbor_matrix, num_neighbors, neighbor_matrix_shifts
 
 
+def _validate_single_segment_coo_state(
+    *,
+    pair_offsets: jax.Array,
+    pair_counts: jax.Array,
+    rebuild_flags: jax.Array | None,
+    neighbor_list: jax.Array | None,
+    neighbor_list_shifts: jax.Array | None,
+    max_pairs: int,
+) -> int:
+    """Validate static single-system segmented-COO metadata and buffers."""
+    if pair_offsets.dtype != jnp.int32 or pair_offsets.shape != (2,):
+        raise ValueError("pair_offsets must have shape (2,) and dtype int32")
+    if pair_counts.dtype != jnp.int32 or pair_counts.shape != (1,):
+        raise ValueError("pair_counts must have shape (1,) and dtype int32")
+    if rebuild_flags is not None and (
+        rebuild_flags.dtype != jnp.bool_ or rebuild_flags.shape != (1,)
+    ):
+        raise ValueError("rebuild_flags must have shape (1,) and dtype bool")
+
+    physical_capacity: int | None = None
+    if neighbor_list is not None:
+        if (
+            neighbor_list.dtype != jnp.int32
+            or neighbor_list.ndim != 2
+            or neighbor_list.shape[0] != 2
+        ):
+            raise ValueError(
+                "neighbor_list must have shape (2, capacity) and dtype int32"
+            )
+        physical_capacity = int(neighbor_list.shape[1])
+    if neighbor_list_shifts is not None:
+        if (
+            neighbor_list_shifts.dtype != jnp.int32
+            or neighbor_list_shifts.ndim != 2
+            or neighbor_list_shifts.shape[1] != 3
+        ):
+            raise ValueError(
+                "neighbor_list_shifts must have shape (capacity, 3) and dtype int32"
+            )
+        shift_capacity = int(neighbor_list_shifts.shape[0])
+        if physical_capacity is not None and shift_capacity != physical_capacity:
+            raise ValueError(
+                "neighbor_list_shifts must match neighbor_list physical capacity"
+            )
+        physical_capacity = shift_capacity
+    if physical_capacity is None:
+        physical_capacity = int(max_pairs)
+    if physical_capacity < 0:
+        raise ValueError("max_pairs must be nonnegative")
+    return physical_capacity
+
+
+def _normalize_single_segment_coo_count(
+    *,
+    pair_offsets: jax.Array,
+    pair_counts: jax.Array,
+    physical_capacity: int,
+) -> jax.Array:
+    """Fail closed for malformed single-system COO metadata without host sync."""
+    offsets_valid = (pair_offsets[0] == 0) & (pair_offsets[1] == physical_capacity)
+    clamped_count = jnp.clip(pair_counts, 0, physical_capacity)
+    return jnp.where(
+        offsets_valid,
+        clamped_count,
+        jnp.zeros_like(pair_counts),
+    )
+
+
 def query_cluster_tile_coo(
     sorted_atom_index: jax.Array,
     sorted_pos_x: jax.Array,
@@ -1417,9 +1487,9 @@ def query_cluster_tile_coo(
 ) -> tuple[jax.Array, ...]:
     """Convert the tile pair list to flat COO form.
 
-    In compact mode the pair count is data-dependent (eager-only).  In
-    segmented mode the output arrays have fixed shapes dictated by
-    ``pair_offsets`` so the call is ``jit``-compatible.
+    In compact mode the pair count is data-dependent (eager-only). In
+    segmented mode caller-provided fixed buffers define output shapes, so the
+    call is ``jit``-compatible.
 
     Parameters
     ----------
@@ -1444,16 +1514,18 @@ def query_cluster_tile_coo(
     natom : int
         Real atom count (excluding padding).
     max_pairs : int
-        Upper bound on the number of pair entries in compact mode.
-        Ignored in segmented mode (derived from ``pair_offsets[-1]``).
+        Upper bound on the number of pair entries in compact mode. In
+        segmented mode it sizes missing fixed buffers; supplied buffer shapes
+        define physical capacity.
     rebuild_flags : jax.Array, shape (1,), dtype=bool, optional
         Selective-rebuild flag.  Requires ``pair_offsets`` and
         ``pair_counts`` (segmented mode only).
-    pair_offsets : jax.Array, shape (ngroup + 1,), dtype=int32, optional
-        CSR-style offsets into the fixed segmented pair buffer.
+    pair_offsets : jax.Array, shape (2,), dtype=int32, optional
+        Exact ``[0, physical_capacity]`` interval for the fixed segmented
+        pair buffer.
         Pass together with ``pair_counts`` to activate segmented mode.
-    pair_counts : jax.Array, shape (ngroup,), dtype=int32, optional
-        Per-group pair counts for the fixed segmented buffer.
+    pair_counts : jax.Array, shape (1,), dtype=int32, optional
+        Filled pair count for the fixed segmented buffer.
     neighbor_list : jax.Array, shape (2, max_pairs), dtype=int32, optional
         Pre-allocated COO list buffer (segmented mode).
     neighbor_list_shifts : jax.Array, shape (max_pairs, 3), dtype=int32, optional
@@ -1484,6 +1556,15 @@ def query_cluster_tile_coo(
         raise ValueError("Pass both 'pair_offsets' and 'pair_counts', or neither.")
     if rebuild_flags is not None and not segmented:
         raise ValueError("rebuild_flags requires pair_offsets and pair_counts")
+    if segmented:
+        physical_capacity = _validate_single_segment_coo_state(
+            pair_offsets=pair_offsets,
+            pair_counts=pair_counts,
+            rebuild_flags=rebuild_flags,
+            neighbor_list=neighbor_list,
+            neighbor_list_shifts=neighbor_list_shifts,
+            max_pairs=max_pairs,
+        )
     coo_registration = _CLUSTER_TILE_QUERIES["coo_segmented" if segmented else "coo"]
     coo_registration.preload(device_source=sorted_pos_x)
 
@@ -1492,11 +1573,13 @@ def query_cluster_tile_coo(
 
     pair_counter = jnp.zeros(1, dtype=jnp.int32)
     if segmented:
-        max_pairs = int(pair_offsets[-1])
         if neighbor_list is None:
-            neighbor_list = jnp.zeros((2, max_pairs), dtype=jnp.int32)
+            neighbor_list = jnp.zeros((2, physical_capacity), dtype=jnp.int32)
         if neighbor_list_shifts is None:
-            neighbor_list_shifts = jnp.zeros((max_pairs, 3), dtype=jnp.int32)
+            neighbor_list_shifts = jnp.zeros(
+                (physical_capacity, 3),
+                dtype=jnp.int32,
+            )
         coo_list = neighbor_list.T.copy()
         coo_shifts = neighbor_list_shifts
         if rebuild_flags is None:
@@ -1521,7 +1604,12 @@ def query_cluster_tile_coo(
             coo_shifts,
             float(cutoff),
             int(natom),
-            int(max_pairs),
+            physical_capacity,
+        )
+        pair_counts = _normalize_single_segment_coo_count(
+            pair_offsets=pair_offsets,
+            pair_counts=pair_counts,
+            physical_capacity=physical_capacity,
         )
         del pair_counter
         return coo_list.T, pair_offsets, pair_counts, coo_shifts
@@ -1708,9 +1796,11 @@ def cluster_tile_neighbor_list(
     max_pairs : int, optional
         Upper bound for compact COO output; defaults to ``N * max_neighbors``.
     pair_offsets, previous_pair_counts, previous_neighbor_list, previous_neighbor_list_shifts : jax.Array, optional
-        Fixed segmented-COO buffers used with ``rebuild_flags`` and
-        ``format="coo"``. The return tuple preserves the fixed buffer
-        shapes and reports the updated ``pair_counts``.
+        Single fixed segmented-COO state used with ``rebuild_flags`` and
+        ``format="coo"``. ``pair_offsets`` must be
+        ``[0, previous_neighbor_list.shape[1]]`` and
+        ``previous_pair_counts`` must have shape ``(1,)``. The return tuple
+        preserves fixed buffer shapes and reports the updated count.
     return_vectors, return_distances : bool, default False
         If True, append per-pair displacement vectors / scalar distances
         to the matrix-format return tuple. Matrix format only.
@@ -1724,7 +1814,10 @@ def cluster_tile_neighbor_list(
     For ``format == "matrix"``:
         ``(neighbor_matrix, num_neighbors, neighbor_matrix_shifts)``, with
         optional ``(*, distances)`` and/or ``(*, vectors)`` appended when
-        ``return_distances`` / ``return_vectors`` is True.
+        ``return_distances`` / ``return_vectors`` is True. With
+        ``rebuild_flags``, append ``(num_tiles, tile_row_group,
+        tile_col_group)`` after one matrix triple, or after both triples when
+        ``cutoff2`` is provided.
     For ``format == "coo"``:
         ``(neighbor_list, neighbor_ptr, neighbor_list_shifts)`` in compact
         mode, or ``(neighbor_list, pair_offsets, pair_counts,
@@ -1832,6 +1925,16 @@ def cluster_tile_neighbor_list(
         max_neighbors = estimate_max_neighbors(
             cutoff2 if cutoff2 is not None else cutoff
         )
+    single_segment_capacity: int | None = None
+    if selective and format == "coo":
+        single_segment_capacity = _validate_single_segment_coo_state(
+            pair_offsets=pair_offsets,
+            pair_counts=previous_pair_counts,
+            rebuild_flags=rebuild_flags,
+            neighbor_list=previous_neighbor_list,
+            neighbor_list_shifts=previous_neighbor_list_shifts,
+            max_pairs=max_pairs if max_pairs is not None else N * max_neighbors,
+        )
 
     if N == 0:
         # Public docstring promises ``N >= 0``.  Returning zero-shaped
@@ -1842,6 +1945,31 @@ def cluster_tile_neighbor_list(
         nn0 = jnp.empty(0, dtype=jnp.int32)
         ns0 = jnp.empty((0, int(max_neighbors), 3), dtype=jnp.int32)
         if format == "coo":
+            if selective:
+                empty_pair_counts = jnp.where(
+                    rebuild_flags,
+                    jnp.zeros_like(previous_pair_counts),
+                    previous_pair_counts,
+                )
+                empty_pair_counts = _normalize_single_segment_coo_count(
+                    pair_offsets=pair_offsets,
+                    pair_counts=empty_pair_counts,
+                    physical_capacity=single_segment_capacity,
+                )
+                empty_num_tiles = jnp.where(
+                    rebuild_flags,
+                    jnp.zeros_like(previous_num_tiles),
+                    previous_num_tiles,
+                )
+                return (
+                    previous_neighbor_list,
+                    pair_offsets,
+                    empty_pair_counts,
+                    previous_neighbor_list_shifts,
+                    empty_num_tiles,
+                    previous_tile_row_group,
+                    previous_tile_col_group,
+                )
             coo_base = (
                 jnp.empty((2, 0), dtype=jnp.int32),
                 jnp.zeros(1, dtype=jnp.int32),
@@ -1903,7 +2031,11 @@ def cluster_tile_neighbor_list(
         if selective:
             return (
                 *matrix_out,
-                previous_num_tiles,
+                jnp.where(
+                    rebuild_flags,
+                    jnp.zeros_like(previous_num_tiles),
+                    previous_num_tiles,
+                ),
                 previous_tile_row_group,
                 previous_tile_col_group,
             )
@@ -2034,7 +2166,7 @@ def cluster_tile_neighbor_list(
                 cell_topology,
                 cutoff,
                 N,
-                int(pair_offsets[-1]),
+                single_segment_capacity,
                 rebuild_flags=rebuild_flags,
                 pair_offsets=pair_offsets,
                 pair_counts=previous_pair_counts,
diff --git a/nvalchemiops/jax/neighbors/naive_dual_cutoff.py b/nvalchemiops/jax/neighbors/naive_dual_cutoff.py
index 4a3fe905..4f2e1b53 100644
--- a/nvalchemiops/jax/neighbors/naive_dual_cutoff.py
+++ b/nvalchemiops/jax/neighbors/naive_dual_cutoff.py
@@ -440,6 +440,83 @@ def naive_neighbor_list_dual_cutoff(
                 per_atom_cell_offsets,
                 launch_dims=(total_atoms,),
             )
+            if rebuild_flags is not None:
+                rf = rebuild_flags.flatten()[:1].astype(jnp.bool_)
+                num_neighbors1 = jnp.where(
+                    rf[0], jnp.zeros_like(num_neighbors1), num_neighbors1
+                )
+                num_neighbors2 = jnp.where(
+                    rf[0], jnp.zeros_like(num_neighbors2), num_neighbors2
+                )
+                (
+                    neighbor_matrix1,
+                    neighbor_matrix_shifts1,
+                    num_neighbors1,
+                    neighbor_matrix2,
+                    neighbor_matrix_shifts2,
+                    num_neighbors2,
+                ) = _DIRECT_NAIVE_DUAL_KERNELS[("wrap_on_entry", True)][
+                    positions.dtype
+                ](
+                    positions_wrapped,
+                    per_atom_cell_offsets,
+                    float(cutoff1 * cutoff1),
+                    float(cutoff2 * cutoff2),
+                    cell,
+                    shift_range_per_dimension,
+                    empty_num_shifts,
+                    empty_batch_idx,
+                    empty_batch_ptr,
+                    empty_target_indices,
+                    neighbor_matrix1,
+                    neighbor_matrix_shifts1,
+                    num_neighbors1,
+                    neighbor_matrix2,
+                    neighbor_matrix_shifts2,
+                    num_neighbors2,
+                    empty_vectors,
+                    empty_distances,
+                    empty_pair_params,
+                    empty_energies,
+                    empty_forces,
+                    rf,
+                    launch_dims=(1, max_shifts_per_system, total_atoms),
+                )
+            else:
+                (
+                    neighbor_matrix1,
+                    neighbor_matrix_shifts1,
+                    num_neighbors1,
+                    neighbor_matrix2,
+                    neighbor_matrix_shifts2,
+                    num_neighbors2,
+                ) = _DIRECT_NAIVE_DUAL_KERNELS[("wrap_on_entry", False)][
+                    positions.dtype
+                ](
+                    positions_wrapped,
+                    per_atom_cell_offsets,
+                    float(cutoff1 * cutoff1),
+                    float(cutoff2 * cutoff2),
+                    cell,
+                    shift_range_per_dimension,
+                    empty_num_shifts,
+                    empty_batch_idx,
+                    empty_batch_ptr,
+                    empty_target_indices,
+                    neighbor_matrix1,
+                    neighbor_matrix_shifts1,
+                    num_neighbors1,
+                    neighbor_matrix2,
+                    neighbor_matrix_shifts2,
+                    num_neighbors2,
+                    empty_vectors,
+                    empty_distances,
+                    empty_pair_params,
+                    empty_energies,
+                    empty_forces,
+                    empty_rebuild_flags,
+                    launch_dims=(1, max_shifts_per_system, total_atoms),
+                )
         else:
             if rebuild_flags is not None:
                 rf = rebuild_flags.flatten()[:1].astype(jnp.bool_)
diff --git a/nvalchemiops/neighbors/cluster_tile/kernels.py b/nvalchemiops/neighbors/cluster_tile/kernels.py
index 5c2dd5ea..0e930a56 100644
--- a/nvalchemiops/neighbors/cluster_tile/kernels.py
+++ b/nvalchemiops/neighbors/cluster_tile/kernels.py
@@ -541,6 +541,18 @@ def _bbox_valid(
     return 0
 
 
+@wp.func
+def _coo_segment_is_valid(
+    start: wp.int32,
+    stop: wp.int32,
+    max_pairs: wp.int32,
+) -> wp.int32:
+    """Return whether a segmented COO interval lies within physical storage."""
+    if start >= 0 and stop >= start and stop <= max_pairs:
+        return wp.int32(1)
+    return wp.int32(0)
+
+
 @lru_cache(maxsize=None)
 def _get_reset_cluster_tile_counts_kernel(*, selective: bool) -> wp.Kernel:
     """Build a counter-reset kernel for cluster-tile segmented state."""
@@ -1245,7 +1257,8 @@ def _kernel(
         natom : wp.int32
             Number of real atoms before padding.
         max_pairs : wp.int32
-            Compact COO capacity. Sentinel in segmented COO specializations.
+            Physical COO output capacity for both compact and segmented modes.
+            Segmented offsets must describe intervals within this capacity.
         tile_row_group, tile_col_group, tile_system : wp.array, shape (num_tiles,), dtype=wp.int32
             Cluster-tile records from the build launcher.
         num_tiles : wp.array, shape (1,), dtype=wp.int32
@@ -1286,10 +1299,11 @@ def _kernel(
         -----
         - Thread launch: Tiled launch with one block per allocated tile slot and block dimension ``TILE_GROUP_SIZE``.
         - Modifies: ``pair_counter`` or ``pair_counts``, ``coo_list``, ``coo_shifts``, and enabled pair-output buffers.
-        Segmented COO writes use ``local = atomic_add(pair_counts, isys,
-        count)`` and write only if ``local < pair_offsets[isys + 1] -
-        pair_offsets[isys]``; overflow is reported by the wrapper from the
-        final count.
+        Segmented COO writes reserve from ``pair_counts`` only for valid output
+        intervals. Non-batched queries require the complete physical interval
+        ``[0, max_pairs]``; batched queries accept physically bounded
+        per-system subsegments. Valid segments retain attempted counts for
+        wrapper overflow reporting; invalid segments retain zero active pairs.
 
         See Also
         --------
@@ -1388,14 +1402,35 @@ def _kernel(
                 base = wp.int32(0)
                 capacity = max_pairs
                 offset = wp.int32(0)
+                segment_valid = wp.int32(1)
                 # Full-fill: reserve two slots per valid pair (forward (i, j)
                 # and reverse (j, i)), so the COO list carries each ordered
                 # pair once -- matching cell_list / naive (half_fill=False).
                 if COO_SEGMENTED:
-                    capacity = pair_offsets[system_idx + 1] - pair_offsets[system_idx]
-                    offset = pair_offsets[system_idx]
-                    if lane == 0:
-                        base = wp.atomic_add(pair_counts, system_idx, 2 * i_count)
+                    start = pair_offsets[system_idx]
+                    stop = pair_offsets[system_idx + 1]
+                    if BATCHED:
+                        segment_valid = _coo_segment_is_valid(
+                            start,
+                            stop,
+                            max_pairs,
+                        )
+                    else:
+                        if start == 0 and stop == max_pairs:
+                            segment_valid = wp.int32(1)
+                        else:
+                            segment_valid = wp.int32(0)
+                    if segment_valid == 1:
+                        capacity = stop - start
+                        offset = start
+                        if lane == 0:
+                            base = wp.atomic_add(
+                                pair_counts,
+                                system_idx,
+                                2 * i_count,
+                            )
+                    else:
+                        capacity = wp.int32(0)
                 else:
                     if lane == 0:
                         base = wp.atomic_add(pair_counter, 0, 2 * i_count)
@@ -1412,7 +1447,13 @@ def _kernel(
                 write_pos = offset + local_pos
                 write_rev = offset + local_rev
 
-                if valid == 1 and local_pos < capacity:
+                if (
+                    valid == 1
+                    and segment_valid == 1
+                    and local_pos < capacity
+                    and write_pos >= 0
+                    and write_pos < max_pairs
+                ):
                     coo_list[write_pos, 0] = i_orig
                     coo_list[write_pos, 1] = j_orig_t
                     coo_shifts[write_pos, 0] = s_shift[0]
@@ -1435,7 +1476,13 @@ def _kernel(
                             pair_energies[write_pos] = pair_energy
                             pair_forces[write_pos] = pair_force
 
-                if valid == 1 and local_rev < capacity:
+                if (
+                    valid == 1
+                    and segment_valid == 1
+                    and local_rev < capacity
+                    and write_rev >= 0
+                    and write_rev < max_pairs
+                ):
                     coo_list[write_rev, 0] = j_orig_t
                     coo_list[write_rev, 1] = i_orig
                     coo_shifts[write_rev, 0] = -s_shift[0]
diff --git a/nvalchemiops/torch/neighbors/__init__.py b/nvalchemiops/torch/neighbors/__init__.py
index 0dce2576..9fe2a8ba 100644
--- a/nvalchemiops/torch/neighbors/__init__.py
+++ b/nvalchemiops/torch/neighbors/__init__.py
@@ -246,6 +246,15 @@ def neighbor_list(
             Boolean flags selecting which systems to re-enumerate; systems whose
             flag is ``False`` keep their previous output (per-system skip for the
             batched methods, whole-list flag for single-system methods).
+        pair_offsets, pair_counts : torch.Tensor, optional
+            Fixed segmented COO metadata for explicit
+            ``method="cluster_tile"`` with ``return_neighbor_list=True`` and
+            ``rebuild_flags``. Single-system tensors have shapes ``(2,)`` and
+            ``(1,)`` int32 and are returned with the fixed-capacity COO buffers.
+        return_state : bool, default=False
+            With ``rebuild_flags`` and an explicit ``method="cluster_tile"`` or
+            ``method="batch_cluster_tile"``, append reusable tile state to the
+            result. See the selected method's return contract for exact tensors.
         pair_fn : warp.Function or CompiledPairFn, optional
             Inline Warp pair potential evaluated as neighbors are enumerated;
             requires ``pair_params`` and fills ``pair_energies`` / ``pair_forces``.
@@ -305,6 +314,22 @@ def neighbor_list(
 
         When ``cutoff2`` is provided, the pattern repeats for the second cutoff with interleaved
         components (neighbor_data2, num_neighbor_data2, neighbor_shift_data2) appended to the tuple.
+        Explicit cluster-tile methods append their documented tile-state suffix
+        when ``return_state=True``.
+
+        Single-system selective COO calls to explicit ``method="cluster_tile"``
+        return ``(neighbor_list, pair_offsets, pair_counts,
+        neighbor_list_shifts)`` instead of the compact COO pointer tuple. With
+        ``return_state=True``, ``(num_tiles, tile_row_group, tile_col_group)``
+        is appended.
+
+        Batched selective calls to explicit ``method="batch_cluster_tile"`` append
+        ``(tile_offsets, tile_counts, num_tiles, tile_row_group, tile_col_group,
+        tile_system)`` when ``return_state=True``. The resulting matrix tuple has
+        nine tensors (or twelve with ``cutoff2``); segmented COO returns
+        ``(neighbor_list, pair_offsets, pair_counts, neighbor_list_shifts,
+        tile_offsets, tile_counts, num_tiles, tile_row_group, tile_col_group,
+        tile_system)``.
 
     Examples
     --------
@@ -352,6 +377,12 @@ def neighbor_list(
         or kwargs.get("pair_forces") is not None
     )
     rebuild_flags = kwargs.get("rebuild_flags")
+    return_state = bool(kwargs.get("return_state", False))
+    if return_state and method is None:
+        raise ValueError(
+            "return_state=True requires an explicit cluster_tile method "
+            "(method='cluster_tile' or method='batch_cluster_tile')"
+        )
     selected_naive_strategy = "auto"
     selected_cell_strategy = "auto"
 
@@ -426,6 +457,10 @@ def _apply_auto_suboptions(
                 method
             )
             _apply_auto_suboptions(fg_native, fg_cell, fg_path)
+    if return_state and method not in ("cluster_tile", "batch_cluster_tile"):
+        raise ValueError(
+            "return_state=True is supported only by explicit cluster_tile methods"
+        )
     match method:
         case "naive":
             return naive_neighbor_list(
diff --git a/nvalchemiops/torch/neighbors/_dispatch.py b/nvalchemiops/torch/neighbors/_dispatch.py
index 1f28a5e7..2c123320 100644
--- a/nvalchemiops/torch/neighbors/_dispatch.py
+++ b/nvalchemiops/torch/neighbors/_dispatch.py
@@ -382,6 +382,13 @@ def _reject_unsupported_cluster_tile_combo(
     half_fill: bool,
 ) -> None:
     """Raise ``NotImplementedError`` for unsupported explicit cluster-tile combos."""
+    if torch.compiler.is_compiling():
+        raise NotImplementedError(
+            "compiled unified cluster-tile dispatch cannot validate tensor-valued "
+            "pbc without synchronization; validate fully periodic pbc eagerly and "
+            "compile the direct single-system cluster_tile_neighbor_list path; "
+            "batched fullgraph is unsupported"
+        )
     if pbc is None:
         raise NotImplementedError(
             "method='cluster_tile' / 'batch_cluster_tile' is "
@@ -389,10 +396,7 @@ def _reject_unsupported_cluster_tile_combo(
             "are not supported.  Use method='naive' or 'cell_list', "
             "or pass a cell with fully periodic pbc."
         )
-    try:
-        all_periodic = bool(pbc.all().item())
-    except RuntimeError:
-        all_periodic = True
+    all_periodic = bool(pbc.all().item())
     if not all_periodic:
         raise NotImplementedError(
             "method='cluster_tile' / 'batch_cluster_tile' is "
diff --git a/nvalchemiops/torch/neighbors/batch_cluster_tile.py b/nvalchemiops/torch/neighbors/batch_cluster_tile.py
index e4e196a0..91e7a045 100644
--- a/nvalchemiops/torch/neighbors/batch_cluster_tile.py
+++ b/nvalchemiops/torch/neighbors/batch_cluster_tile.py
@@ -60,6 +60,7 @@
     selective_zero_num_neighbors as wp_selective_zero_num_neighbors,
 )
 from nvalchemiops.neighbors.output_args import _has_partial_or_pair_outputs
+from nvalchemiops.torch.neighbors.neighbor_utils import _validate_segmented_coo_state
 from nvalchemiops.torch.types import get_wp_dtype
 
 if TYPE_CHECKING:
@@ -1624,6 +1625,25 @@ def batch_query_cluster_tile_coo(
         Output pair forces; modified in-place when ``pair_fn`` is set.
     """
 
+    if pair_offsets is not None:
+        if pair_counts is None:
+            raise ValueError("pair_counts is required with pair_offsets")
+        _validate_segmented_coo_state(
+            device=coo_list.device,
+            num_systems=int(pair_offsets.shape[0]) - 1,
+            neighbor_list=coo_list.transpose(0, 1),
+            neighbor_list_shifts=coo_shifts,
+            pair_offsets=pair_offsets,
+            pair_counts=pair_counts,
+            rebuild_flags=rebuild_flags,
+            tile_offsets=tile_offsets,
+            tile_counts=tile_counts,
+            num_tiles=num_tiles,
+            tile_row_group=tile_row_group,
+            tile_col_group=tile_col_group,
+            tile_system=tile_system,
+        )
+
     if inv_cell_batch is None:
         inv_cell_batch = torch.linalg.inv(cell_batch).contiguous()
     pair_counter.zero_()
@@ -1851,6 +1871,75 @@ def _batch_cluster_tile_pair_outputs_forward(
     )
 
 
+def _validate_batch_matrix_state(
+    *,
+    device: torch.device,
+    natom: int,
+    num_systems: int,
+    max_neighbors: int,
+    neighbor_matrix: torch.Tensor,
+    num_neighbors: torch.Tensor,
+    neighbor_matrix_shifts: torch.Tensor,
+    num_tiles: torch.Tensor,
+    tile_offsets: torch.Tensor,
+    tile_counts: torch.Tensor,
+    tile_row_group: torch.Tensor,
+    tile_col_group: torch.Tensor,
+    tile_system: torch.Tensor,
+    neighbor_matrix2: torch.Tensor | None = None,
+    num_neighbors2: torch.Tensor | None = None,
+    neighbor_matrix_shifts2: torch.Tensor | None = None,
+) -> None:
+    """Validate caller-owned selective matrix and tile state."""
+    state = {
+        "neighbor_matrix": neighbor_matrix,
+        "num_neighbors": num_neighbors,
+        "neighbor_matrix_shifts": neighbor_matrix_shifts,
+        "num_tiles": num_tiles,
+        "tile_offsets": tile_offsets,
+        "tile_counts": tile_counts,
+        "tile_row_group": tile_row_group,
+        "tile_col_group": tile_col_group,
+        "tile_system": tile_system,
+    }
+    for name, value in state.items():
+        if value.device != device or value.dtype != torch.int32:
+            raise ValueError(f"{name} must be an int32 tensor on positions.device")
+    if neighbor_matrix.shape != (natom, max_neighbors):
+        raise ValueError("neighbor_matrix must have shape (N, max_neighbors)")
+    if num_neighbors.shape != (natom,):
+        raise ValueError("num_neighbors must have shape (N,)")
+    if neighbor_matrix_shifts.shape != (natom, max_neighbors, 3):
+        raise ValueError("neighbor_matrix_shifts must have shape (N, max_neighbors, 3)")
+    if (
+        tile_offsets.shape != (num_systems + 1,)
+        or tile_counts.shape != (num_systems,)
+        or num_tiles.shape != (1,)
+        or tile_row_group.ndim != 1
+        or tile_col_group.shape != tile_row_group.shape
+        or tile_system.shape != tile_row_group.shape
+    ):
+        raise ValueError("batched tile state has an invalid shape")
+    for name, matrix, counts, shifts in (
+        ("secondary", neighbor_matrix2, num_neighbors2, neighbor_matrix_shifts2),
+    ):
+        if matrix is not None:
+            if (
+                counts is None
+                or shifts is None
+                or matrix.device != device
+                or matrix.dtype != torch.int32
+                or counts.device != device
+                or counts.dtype != torch.int32
+                or shifts.device != device
+                or shifts.dtype != torch.int32
+                or matrix.shape != (natom, max_neighbors)
+                or counts.shape != (natom,)
+                or shifts.shape != (natom, max_neighbors, 3)
+            ):
+                raise ValueError(f"{name} matrix state has an invalid shape or dtype")
+
+
 def batch_cluster_tile_neighbor_list(
     positions: torch.Tensor,
     cutoff: float,
@@ -1862,6 +1951,7 @@ def batch_cluster_tile_neighbor_list(
     max_pairs: int | None = None,
     cutoff2: float | None = None,
     rebuild_flags: torch.Tensor | None = None,
+    return_state: bool = False,
     # matrix-format outputs
     neighbor_matrix: torch.Tensor | None = None,
     neighbor_matrix_shifts: torch.Tensor | None = None,
@@ -1949,6 +2039,12 @@ def batch_cluster_tile_neighbor_list(
     rebuild_flags : torch.Tensor, shape (num_systems,), dtype=torch.bool, optional
         Per-system selective rebuild flags. Supported for matrix output and
         segmented COO output.
+    return_state : bool, default=False
+        When used with ``rebuild_flags``, append the reusable tile state
+        ``(tile_offsets, tile_counts, num_tiles, tile_row_group,
+        tile_col_group, tile_system)`` to matrix or segmented-COO output.
+        These are the same tensor objects as caller-supplied state buffers.
+        Requires ``rebuild_flags``.
     tile_offsets, tile_counts : torch.Tensor, optional
         Fixed per-system tile offsets and OUTPUT tile counters for segmented
         tile-list state. Use ``estimate_batch_cluster_tile_segments`` to size
@@ -1956,6 +2052,10 @@ def batch_cluster_tile_neighbor_list(
     pair_offsets, pair_counts : torch.Tensor, optional
         Fixed per-system COO offsets and OUTPUT pair counters for segmented
         COO output. Compact COO cannot be combined with ``rebuild_flags``.
+        Before launch, offsets, counts, topology buffers, and tile state are
+        validated against their physical capacities. Eager calls containing a
+        false flag require caller-owned topology and tile state; batched
+        ``torch.compile(fullgraph=True)`` selective COO is not supported.
     neighbor_matrix, num_neighbors, neighbor_matrix_shifts : torch.Tensor, optional
         Pre-allocated matrix-format outputs. All-or-nothing across the trio.
     neighbor_list, neighbor_list_shifts, pair_counter : torch.Tensor, optional
@@ -1994,13 +2094,18 @@ def batch_cluster_tile_neighbor_list(
           neighbor_matrix_shifts)``, with optional ``(*, distances)``
           and/or ``(*, vectors)`` appended when ``return_distances`` /
           ``return_vectors`` is True, and optional ``(*, pair_energies,
-          pair_forces)`` when ``pair_fn`` is set.
+          pair_forces)`` when ``pair_fn`` is set. Selective calls with
+          ``return_state=True`` append ``(tile_offsets, tile_counts,
+          num_tiles, tile_row_group, tile_col_group, tile_system)``, yielding
+          nine tensors for one cutoff or twelve for two cutoffs.
         - ``"coo"``: ``(neighbor_list, neighbor_ptr, neighbor_list_shifts)``
           via the direct ``batch_query_cluster_tile_coo`` path (no
           matrix intermediate). ``neighbor_ptr`` is reconstructed from
           ``bincount(neighbor_list[0])``. With segmented COO, returns
           ``(neighbor_list, pair_offsets, pair_counts, neighbor_list_shifts)``
-          without trimming the caller-owned fixed segments.
+          without trimming the caller-owned fixed segments; selective calls
+          with ``return_state=True`` append the six state tensors above,
+          yielding ten tensors.
         - ``"tile"``: 11-tuple ``(num_tiles, tile_row_group,
           tile_col_group, tile_system, sorted_atom_index, sorted_pos_x,
           sorted_pos_y, sorted_pos_z, batch_idx_sorted, batch_ptr_padded,
@@ -2043,6 +2148,8 @@ def batch_cluster_tile_neighbor_list(
         or pair_energies is not None
         or pair_forces is not None
     )
+    if return_state and rebuild_flags is None:
+        raise ValueError("return_state=True requires rebuild_flags")
     if has_pair_outputs and format == "tile":
         raise NotImplementedError(
             "Pair outputs (return_vectors / return_distances / pair_fn) "
@@ -2066,6 +2173,11 @@ def batch_cluster_tile_neighbor_list(
         rebuild_flags is not None
         and format == "coo"
         and (pair_offsets is None or pair_counts is None)
+        and (
+            torch.compiler.is_compiling()
+            or not bool(rebuild_flags.all().item())
+            or not return_state
+        )
     ):
         raise ValueError(
             "cluster_tile selective COO requires pair_offsets and pair_counts"
@@ -2078,6 +2190,85 @@ def batch_cluster_tile_neighbor_list(
         raise ValueError("batch_ptr must have length at least 2")
     device = positions.device
     N = int(batch_ptr[-1].item())
+    num_systems = int(batch_ptr.shape[0]) - 1
+    is_compiling = torch.compiler.is_compiling()
+    if rebuild_flags is not None and format == "coo":
+        if neighbor_list is not None and pair_offsets is None and pair_counts is None:
+            raise ValueError(
+                "selective segmented COO bootstrap requires complete caller-owned "
+                "state or no persistent state"
+            )
+    eager_has_false = (
+        rebuild_flags is not None
+        and not is_compiling
+        and not bool(rebuild_flags.all().item())
+    )
+
+    if eager_has_false:
+        required_state: dict[str, torch.Tensor | None] = {
+            "num_tiles": num_tiles,
+            "tile_row_group": tile_row_group,
+            "tile_col_group": tile_col_group,
+            "tile_system": tile_system,
+            "tile_offsets": tile_offsets,
+            "tile_counts": tile_counts,
+        }
+        if format == "coo":
+            required_state.update(
+                {
+                    "neighbor_list": neighbor_list,
+                    "neighbor_list_shifts": neighbor_list_shifts,
+                    "pair_offsets": pair_offsets,
+                    "pair_counts": pair_counts,
+                }
+            )
+        else:
+            required_state.update(
+                {
+                    "neighbor_matrix": neighbor_matrix,
+                    "num_neighbors": num_neighbors,
+                    "neighbor_matrix_shifts": neighbor_matrix_shifts,
+                }
+            )
+            if cutoff2 is not None:
+                required_state.update(
+                    {
+                        "neighbor_matrix2": neighbor_matrix2,
+                        "num_neighbors2": num_neighbors2,
+                        "neighbor_matrix_shifts2": neighbor_matrix_shifts2,
+                    }
+                )
+        missing_state = [
+            name for name, value in required_state.items() if value is None
+        ]
+        if missing_state:
+            raise ValueError(
+                "selective segmented COO with rebuild_flags=False requires "
+                "caller-owned topology and tile state: " + ", ".join(missing_state)
+            )
+    if is_compiling and rebuild_flags is not None and format == "matrix":
+        compiled_state = (
+            neighbor_matrix,
+            num_neighbors,
+            neighbor_matrix_shifts,
+            num_tiles,
+            tile_offsets,
+            tile_counts,
+            tile_row_group,
+            tile_col_group,
+            tile_system,
+        )
+        if any(value is None for value in compiled_state):
+            raise ValueError(
+                "selective compiled matrix requires complete caller-owned state"
+            )
+        if cutoff2 is not None and any(
+            value is None
+            for value in (neighbor_matrix2, num_neighbors2, neighbor_matrix_shifts2)
+        ):
+            raise ValueError(
+                "selective compiled dual-cutoff matrix requires both output triples"
+            )
 
     if max_neighbors is None:
         max_neighbors = max(
@@ -2086,20 +2277,64 @@ def batch_cluster_tile_neighbor_list(
         )
     if fill_value is None:
         fill_value = N
+    build_cutoff = cutoff if cutoff2 is None else max(float(cutoff), float(cutoff2))
+    if max_tiles_per_group is None:
+        max_tiles_per_group = estimate_batch_max_tiles_per_group(
+            batch_ptr, build_cutoff, cell_batch
+        )
+    if (
+        rebuild_flags is not None
+        and format == "matrix"
+        and neighbor_matrix is not None
+        and num_neighbors is not None
+        and neighbor_matrix_shifts is not None
+        and num_tiles is not None
+        and tile_offsets is not None
+        and tile_counts is not None
+        and tile_row_group is not None
+        and tile_col_group is not None
+        and tile_system is not None
+    ):
+        _validate_batch_matrix_state(
+            device=device,
+            natom=N,
+            num_systems=num_systems,
+            max_neighbors=int(max_neighbors),
+            neighbor_matrix=neighbor_matrix,
+            num_neighbors=num_neighbors,
+            neighbor_matrix_shifts=neighbor_matrix_shifts,
+            num_tiles=num_tiles,
+            tile_offsets=tile_offsets,
+            tile_counts=tile_counts,
+            tile_row_group=tile_row_group,
+            tile_col_group=tile_col_group,
+            tile_system=tile_system,
+            neighbor_matrix2=neighbor_matrix2,
+            num_neighbors2=num_neighbors2,
+            neighbor_matrix_shifts2=neighbor_matrix_shifts2,
+        )
 
     needs_segmented_tiles = rebuild_flags is not None or pair_offsets is not None
     if needs_segmented_tiles and tile_offsets is None:
         _tile_caps, tile_offsets, _pair_caps, default_pair_offsets = (
-            estimate_batch_cluster_tile_segments(batch_ptr, int(max_neighbors))
+            estimate_batch_cluster_tile_segments(
+                batch_ptr,
+                int(max_neighbors),
+                max_tiles_per_group=int(max_tiles_per_group),
+            )
         )
-        tile_counts = torch.empty(
+        tile_counts = torch.zeros(
             int(tile_offsets.shape[0]) - 1, dtype=torch.int32, device=device
         )
         if format == "coo" and pair_offsets is None:
             pair_offsets = default_pair_offsets
-            pair_counts = torch.empty(
+            pair_counts = torch.zeros(
                 int(pair_offsets.shape[0]) - 1, dtype=torch.int32, device=device
             )
+            if max_pairs is not None and int(max_pairs) != int(pair_offsets[-1].item()):
+                raise ValueError(
+                    "max_pairs must equal the estimated segmented COO capacity"
+                )
     elif needs_segmented_tiles and tile_counts is None:
         raise ValueError("tile_counts is required when tile_offsets is provided")
 
@@ -2141,16 +2376,13 @@ def batch_cluster_tile_neighbor_list(
             return (*base, distances_out)
         return (*base, vectors_out)
 
-    # Tile candidates must cover the larger radius so the cutoff2 matrix cannot
-    # miss pairs in the (cutoff, cutoff2] shell; the query filters each matrix
-    # by its own cutoff.
-    build_cutoff = cutoff if cutoff2 is None else max(float(cutoff), float(cutoff2))
-
     if sorted_atom_index is None:
-        if max_tiles_per_group is None:
-            max_tiles_per_group = estimate_batch_max_tiles_per_group(
-                batch_ptr, build_cutoff, cell_batch
-            )
+        previous_tile_state = (
+            num_tiles,
+            tile_row_group,
+            tile_col_group,
+            tile_system,
+        )
         (
             sorted_atom_index,
             sort_inv,
@@ -2177,6 +2409,52 @@ def batch_cluster_tile_neighbor_list(
             dtype=positions.dtype,
             max_tiles_per_group=max_tiles_per_group,
         )
+        (
+            previous_num_tiles,
+            previous_tile_row_group,
+            previous_tile_col_group,
+            previous_tile_system,
+        ) = previous_tile_state
+        if previous_num_tiles is not None:
+            num_tiles = previous_num_tiles
+        if previous_tile_row_group is not None:
+            tile_row_group = previous_tile_row_group
+        if previous_tile_col_group is not None:
+            tile_col_group = previous_tile_col_group
+        if previous_tile_system is not None:
+            tile_system = previous_tile_system
+
+    segmented_pair_capacity: int | None = None
+    if format == "coo" and pair_offsets is not None:
+        if neighbor_list is None:
+            capacity = int(pair_offsets[-1].item())
+            neighbor_list = torch.empty(
+                (2, capacity),
+                dtype=torch.int32,
+                device=device,
+            )
+        if neighbor_list_shifts is None:
+            neighbor_list_shifts = torch.empty(
+                (neighbor_list.shape[1], 3),
+                dtype=torch.int32,
+                device=device,
+            )
+        segmented_pair_capacity = _validate_segmented_coo_state(
+            device=device,
+            num_systems=num_systems,
+            neighbor_list=neighbor_list,
+            neighbor_list_shifts=neighbor_list_shifts,
+            pair_offsets=pair_offsets,
+            pair_counts=pair_counts,
+            rebuild_flags=rebuild_flags,
+            tile_offsets=tile_offsets,
+            tile_counts=tile_counts,
+            num_tiles=num_tiles,
+            tile_row_group=tile_row_group,
+            tile_col_group=tile_col_group,
+            tile_system=tile_system,
+        )
+
     batch_build_cluster_tile_list(
         positions,
         build_cutoff,
@@ -2231,7 +2509,9 @@ def batch_cluster_tile_neighbor_list(
     if format == "coo":
         segmented_coo = pair_offsets is not None
         if segmented_coo:
-            max_pairs = int(pair_offsets[-1].item())
+            if segmented_pair_capacity is None:
+                raise RuntimeError("segmented COO state was not validated")
+            max_pairs = segmented_pair_capacity
         elif max_pairs is None:
             max_pairs = N * max_neighbors
         if neighbor_list is None:
@@ -2295,12 +2575,23 @@ def batch_cluster_tile_neighbor_list(
                     int(pair_counts[isys].item()),
                     system_index=isys,
                 )
-            return (
-                coo_buf.transpose(0, 1).contiguous(),
+            outputs = (
+                neighbor_list,
                 pair_offsets,
                 pair_counts,
                 neighbor_list_shifts,
             )
+            if return_state:
+                return (
+                    *outputs,
+                    tile_offsets,
+                    tile_counts,
+                    num_tiles,
+                    tile_row_group,
+                    tile_col_group,
+                    tile_system,
+                )
+            return outputs
 
         npairs = int(pair_counter.item())
         if npairs > int(max_pairs):
@@ -2441,7 +2732,7 @@ def batch_cluster_tile_neighbor_list(
             )
 
     if cutoff2 is not None:
-        return (
+        outputs = (
             neighbor_matrix,
             num_neighbors,
             neighbor_matrix_shifts,
@@ -2449,4 +2740,16 @@ def batch_cluster_tile_neighbor_list(
             num_neighbors2,
             neighbor_matrix_shifts2,
         )
-    return neighbor_matrix, num_neighbors, neighbor_matrix_shifts
+    else:
+        outputs = (neighbor_matrix, num_neighbors, neighbor_matrix_shifts)
+    if return_state:
+        return (
+            *outputs,
+            tile_offsets,
+            tile_counts,
+            num_tiles,
+            tile_row_group,
+            tile_col_group,
+            tile_system,
+        )
+    return outputs
diff --git a/nvalchemiops/torch/neighbors/cluster_tile.py b/nvalchemiops/torch/neighbors/cluster_tile.py
index bc91af63..fb50e3f6 100644
--- a/nvalchemiops/torch/neighbors/cluster_tile.py
+++ b/nvalchemiops/torch/neighbors/cluster_tile.py
@@ -56,6 +56,10 @@
     selective_zero_num_neighbors_single as wp_selective_zero_num_neighbors_single,
 )
 from nvalchemiops.neighbors.output_args import _has_partial_or_pair_outputs
+from nvalchemiops.torch.neighbors.neighbor_utils import (
+    _normalize_compiled_single_segment_coo_count,
+    _validate_segmented_coo_state,
+)
 from nvalchemiops.torch.types import get_wp_dtype
 
 if TYPE_CHECKING:
@@ -382,6 +386,8 @@ def _build_cluster_tile_list_op(
     num_tiles: torch.Tensor,
     tile_row_group: torch.Tensor,
     tile_col_group: torch.Tensor,
+    rebuild_flags: torch.Tensor,
+    use_rebuild_flags: bool,
 ) -> None:
     """Compute Morton codes + argsort + SoA gather in torch, then run
     bbox reduction + tile enumeration on the warp side.
@@ -401,7 +407,8 @@ def _build_cluster_tile_list_op(
     wp_device = str(device)
     wp_dtype = get_wp_dtype(positions.dtype)
 
-    num_tiles.zero_()
+    if not use_rebuild_flags:
+        num_tiles.zero_()
 
     wp_cell = _mat33f_from_torch(cell)
     wp_inv_cell = _mat33f_from_torch(inv_cell)
@@ -487,6 +494,11 @@ def _spread(x: torch.Tensor) -> torch.Tensor:
         dtype=wp.int32,
         return_ctype=True,
     )
+    wp_rebuild_flags = wp.from_torch(
+        rebuild_flags,
+        dtype=wp.bool,
+        return_ctype=True,
+    )
     wp_build_cluster_tile_list(
         sorted_pos_x=wp_sorted_pos_x,
         sorted_pos_y=wp_sorted_pos_y,
@@ -505,6 +517,7 @@ def _spread(x: torch.Tensor) -> torch.Tensor:
         group_ext_x_buffer=wp_group_ext_x,
         group_ext_y_buffer=wp_group_ext_y,
         group_ext_z_buffer=wp_group_ext_z,
+        rebuild_flags=wp_rebuild_flags if use_rebuild_flags else None,
     )
 
 
@@ -528,6 +541,8 @@ def _(
     num_tiles: torch.Tensor,
     tile_row_group: torch.Tensor,
     tile_col_group: torch.Tensor,
+    rebuild_flags: torch.Tensor,
+    use_rebuild_flags: bool,
 ) -> None:
     return None
 
@@ -550,6 +565,8 @@ def build_cluster_tile_list(
     num_tiles: torch.Tensor,
     tile_row_group: torch.Tensor,
     tile_col_group: torch.Tensor,
+    *,
+    rebuild_flags: torch.Tensor | None = None,
 ) -> None:
     """Build cluster-tile neighbor list state into pre-allocated tensors.
 
@@ -610,6 +627,7 @@ def build_cluster_tile_list(
     cell_mat, inv_cell_mat = _cell_invcell_from_cell(cell)
     cell_mat = cell_mat.to(positions.dtype)
     inv_cell_mat = inv_cell_mat.to(positions.dtype)
+    dummy_rebuild_flags = torch.empty(1, dtype=torch.bool, device=positions.device)
     _build_cluster_tile_list_op(
         positions,
         cutoff,
@@ -629,6 +647,8 @@ def build_cluster_tile_list(
         num_tiles,
         tile_row_group,
         tile_col_group,
+        rebuild_flags if rebuild_flags is not None else dummy_rebuild_flags,
+        rebuild_flags is not None,
     )
 
 
@@ -1265,6 +1285,87 @@ def _(
     return None
 
 
+@torch.library.custom_op(
+    "nvalchemiops::_query_cluster_tile_coo_segmented",
+    mutates_args=("pair_counter", "pair_counts", "coo_list", "coo_shifts"),
+)
+def _query_cluster_tile_coo_segmented_op(
+    cutoff: float,
+    natom: int,
+    max_pairs: int,
+    cell: torch.Tensor,
+    inv_cell: torch.Tensor,
+    sorted_atom_index: torch.Tensor,
+    sorted_pos_x: torch.Tensor,
+    sorted_pos_y: torch.Tensor,
+    sorted_pos_z: torch.Tensor,
+    num_tiles: torch.Tensor,
+    tile_row_group: torch.Tensor,
+    tile_col_group: torch.Tensor,
+    rebuild_flags: torch.Tensor,
+    pair_counter: torch.Tensor,
+    pair_offsets: torch.Tensor,
+    pair_counts: torch.Tensor,
+    coo_list: torch.Tensor,
+    coo_shifts: torch.Tensor,
+    n_tiles: int,
+) -> None:
+    """Run fixed-capacity selective COO conversion in place."""
+    device = sorted_pos_x.device
+    wp_device = str(device)
+    wp_dtype = get_wp_dtype(sorted_pos_x.dtype)
+    wp_query_cluster_tile_coo(
+        sorted_atom_index=wp.from_torch(
+            sorted_atom_index, dtype=wp.int32, return_ctype=True
+        ),
+        sorted_pos_x=wp.from_torch(sorted_pos_x, dtype=wp_dtype, return_ctype=True),
+        sorted_pos_y=wp.from_torch(sorted_pos_y, dtype=wp_dtype, return_ctype=True),
+        sorted_pos_z=wp.from_torch(sorted_pos_z, dtype=wp_dtype, return_ctype=True),
+        num_tiles=wp.from_torch(num_tiles, dtype=wp.int32, return_ctype=True),
+        tile_row_group=wp.from_torch(tile_row_group, dtype=wp.int32, return_ctype=True),
+        tile_col_group=wp.from_torch(tile_col_group, dtype=wp.int32, return_ctype=True),
+        cell=_mat33f_from_torch(cell),
+        inv_cell=_mat33f_from_torch(inv_cell),
+        cutoff=float(cutoff),
+        natom=int(natom),
+        max_pairs=int(max_pairs),
+        pair_counter=wp.from_torch(pair_counter, dtype=wp.int32, return_ctype=True),
+        coo_list=wp.from_torch(coo_list, dtype=wp.int32, return_ctype=True),
+        coo_shifts=wp.from_torch(coo_shifts, dtype=wp.int32, return_ctype=True),
+        wp_dtype=wp_dtype,
+        device=wp_device,
+        n_tiles=int(n_tiles),
+        rebuild_flags=wp.from_torch(rebuild_flags, dtype=wp.bool, return_ctype=True),
+        pair_offsets=wp.from_torch(pair_offsets, dtype=wp.int32, return_ctype=True),
+        pair_counts=wp.from_torch(pair_counts, dtype=wp.int32, return_ctype=True),
+    )
+
+
+@_query_cluster_tile_coo_segmented_op.register_fake
+def _(
+    cutoff: float,
+    natom: int,
+    max_pairs: int,
+    cell: torch.Tensor,
+    inv_cell: torch.Tensor,
+    sorted_atom_index: torch.Tensor,
+    sorted_pos_x: torch.Tensor,
+    sorted_pos_y: torch.Tensor,
+    sorted_pos_z: torch.Tensor,
+    num_tiles: torch.Tensor,
+    tile_row_group: torch.Tensor,
+    tile_col_group: torch.Tensor,
+    rebuild_flags: torch.Tensor,
+    pair_counter: torch.Tensor,
+    pair_offsets: torch.Tensor,
+    pair_counts: torch.Tensor,
+    coo_list: torch.Tensor,
+    coo_shifts: torch.Tensor,
+    n_tiles: int,
+) -> None:
+    return None
+
+
 @torch.library.custom_op(
     "nvalchemiops::_query_cluster_tile_coo_optional_no_pair_fn",
     mutates_args=(
@@ -1464,6 +1565,9 @@ def query_cluster_tile_coo(
     coo_list: torch.Tensor,
     coo_shifts: torch.Tensor,
     *,
+    rebuild_flags: torch.Tensor | None = None,
+    pair_offsets: torch.Tensor | None = None,
+    pair_counts: torch.Tensor | None = None,
     return_vectors: bool = False,
     return_distances: bool = False,
     pair_fn: wp.Function | None = None,
@@ -1476,8 +1580,10 @@ def query_cluster_tile_coo(
     """Convert the tile pair list to flat COO format in place.
 
     Cluster-tile does not support partial neighbor lists; there is no
-    ``target_indices`` kwarg.  Optional pair outputs use flat COO buffers of
+    ``target_indices`` kwarg. Optional pair outputs use flat COO buffers of
     length ``max_pairs`` and are written in the same order as ``coo_list``.
+    Selective segmented COO is topology-only and cannot combine with pair
+    outputs.
 
     Parameters
     ----------
@@ -1512,6 +1618,15 @@ def query_cluster_tile_coo(
         Output flat pair list; each row is ``(i, j)``. Modified in-place.
     coo_shifts : torch.Tensor, shape (max_pairs, 3), dtype=int32
         Output periodic image shift vectors for each pair. Modified in-place.
+    rebuild_flags : torch.Tensor, shape (1,), dtype=bool, optional
+        Enables selective fixed-capacity COO. Requires ``pair_offsets`` and
+        ``pair_counts``. A false flag preserves the supplied topology buffers;
+        a true flag rebuilds their active prefix.
+    pair_offsets : torch.Tensor, shape (2,), dtype=int32, optional
+        Fixed COO segment boundaries ``[0, max_pairs]``. Must be supplied with
+        ``rebuild_flags`` and ``pair_counts``.
+    pair_counts : torch.Tensor, shape (1,), dtype=int32, optional
+        Number of active pairs in the segment. Modified in-place.
     return_vectors : bool, optional
         Write per-pair Cartesian displacement vectors to ``neighbor_vectors``.
         Default is ``False``.
@@ -1542,6 +1657,42 @@ def query_cluster_tile_coo(
         Row-padded matrix-format alternative.
     """
 
+    segmented_inputs = (rebuild_flags, pair_offsets, pair_counts)
+    segmented = rebuild_flags is not None
+    if any(value is not None for value in segmented_inputs) and not all(
+        value is not None for value in segmented_inputs
+    ):
+        raise ValueError(
+            "rebuild_flags, pair_offsets, and pair_counts must be supplied together"
+        )
+    if segmented:
+        if _has_partial_or_pair_outputs(
+            return_vectors=return_vectors,
+            return_distances=return_distances,
+            pair_fn=pair_fn,
+            pair_params=pair_params,
+            neighbor_vectors=neighbor_vectors,
+            neighbor_distances=neighbor_distances,
+            pair_energies=pair_energies,
+            pair_forces=pair_forces,
+        ):
+            raise ValueError(
+                "cluster_tile selective rebuild cannot be combined with pair outputs"
+            )
+        if coo_list.shape != (max_pairs, 2):
+            raise ValueError("coo_list must have shape (max_pairs, 2)")
+        physical_capacity = _validate_segmented_coo_state(
+            device=coo_list.device,
+            num_systems=1,
+            neighbor_list=coo_list.transpose(0, 1),
+            neighbor_list_shifts=coo_shifts,
+            pair_offsets=pair_offsets,
+            pair_counts=pair_counts,
+            rebuild_flags=rebuild_flags,
+        )
+        if physical_capacity != max_pairs:
+            raise ValueError("max_pairs must equal coo_list.shape[0]")
+
     cell_mat, inv_cell_mat = _cell_invcell_from_cell(cell)
     cell_mat = cell_mat.to(sorted_pos_x.dtype)
     inv_cell_mat = inv_cell_mat.to(sorted_pos_x.dtype)
@@ -1557,6 +1708,36 @@ def query_cluster_tile_coo(
             raise NeighborOverflowError(tile_capacity, n_tiles)
     pair_counter.zero_()
 
+    if segmented:
+        _query_cluster_tile_coo_segmented_op(
+            cutoff,
+            natom,
+            max_pairs,
+            cell_mat,
+            inv_cell_mat,
+            sorted_atom_index,
+            sorted_pos_x,
+            sorted_pos_y,
+            sorted_pos_z,
+            num_tiles,
+            tile_row_group,
+            tile_col_group,
+            rebuild_flags,
+            pair_counter,
+            pair_offsets,
+            pair_counts,
+            coo_list,
+            coo_shifts,
+            n_tiles,
+        )
+        if torch.compiler.is_compiling():
+            _normalize_compiled_single_segment_coo_count(
+                pair_offsets=pair_offsets,
+                pair_counts=pair_counts,
+                physical_capacity=physical_capacity,
+            )
+        return
+
     if _has_partial_or_pair_outputs(
         return_vectors=return_vectors,
         return_distances=return_distances,
@@ -1765,6 +1946,7 @@ def cluster_tile_neighbor_list(
     max_pairs: int | None = None,
     cutoff2: float | None = None,
     rebuild_flags: torch.Tensor | None = None,
+    return_state: bool = False,
     # matrix-format outputs
     neighbor_matrix: torch.Tensor | None = None,
     neighbor_matrix_shifts: torch.Tensor | None = None,
@@ -1775,6 +1957,8 @@ def cluster_tile_neighbor_list(
     # coo-format outputs
     neighbor_list: torch.Tensor | None = None,
     neighbor_list_shifts: torch.Tensor | None = None,
+    pair_offsets: torch.Tensor | None = None,
+    pair_counts: torch.Tensor | None = None,
     pair_counter: torch.Tensor | None = None,
     # scratch buffers
     sorted_atom_index: torch.Tensor | None = None,
@@ -1838,13 +2022,12 @@ def cluster_tile_neighbor_list(
           ``(neighbor_matrix, num_neighbors, neighbor_matrix_shifts)`` —
           the dense ``(N, max_neighbors)`` row-padded form used by
           ``cell_list`` and ``naive``.
-        - ``"coo"``: returns
-          ``(neighbor_list, neighbor_ptr, neighbor_list_shifts)`` — flat
-          pair list emitted directly by ``query_cluster_tile_coo`` (no matrix
-          intermediate).  ``neighbor_ptr`` is reconstructed from
-          ``bincount(neighbor_list[0])`` (cheap; requires a single CPU
-          sync on ``pair_counter[0]`` that's needed for the trim
-          anyway).
+        - ``"coo"``: compact calls return
+          ``(neighbor_list, neighbor_ptr, neighbor_list_shifts)`` — a trimmed
+          flat pair list emitted directly by ``query_cluster_tile_coo``. With
+          ``rebuild_flags``, fixed-capacity segmented COO returns
+          ``(neighbor_list, pair_offsets, pair_counts,
+          neighbor_list_shifts)`` instead.
         - ``"tile"``: returns the native cluster-pair tile state as a
           7-tuple
           ``(num_tiles, tile_row_group, tile_col_group,
@@ -1857,16 +2040,26 @@ def cluster_tile_neighbor_list(
     max_pairs : int, optional
         Upper bound for COO output; defaults to ``N * max_neighbors``.
     rebuild_flags : torch.Tensor, shape (1,), dtype=bool, optional
-        Matrix-format selective rebuild flag. Requires previous tile state
-        and previous matrix outputs. When ``rebuild_flags[0]`` is False,
-        the previous outputs are returned unchanged.
+        Selective rebuild flag. An eager all-true call may bootstrap omitted
+        state. When false, caller-owned fixed topology and tile state are
+        required before launch. COO calls additionally require
+        ``pair_offsets`` and ``pair_counts``. The explicit single-system COO
+        route supports ``torch.compile(fullgraph=True)`` only after eager
+        bootstrap with fixed-shape validated buffers.
+    return_state : bool, default=False
+        With ``rebuild_flags``, append ``(num_tiles, tile_row_group,
+        tile_col_group)`` after topology outputs. These are the same tensor
+        objects as caller-supplied state buffers.
     neighbor_matrix, num_neighbors, neighbor_matrix_shifts : optional
         Pre-allocated matrix-format outputs.  All-or-nothing only across
         the trio; supply all three or none.
     neighbor_list, neighbor_list_shifts, pair_counter : optional
-        Pre-allocated COO-format outputs.  Shapes
-        ``(2, max_pairs)``, ``(max_pairs, 3)``, ``(1,)`` int32.
-        Same all-or-nothing semantics.
+        Pre-allocated compact COO outputs with shapes ``(2, max_pairs)``,
+        ``(max_pairs, 3)``, ``(1,)`` int32.
+    pair_offsets, pair_counts : optional
+        Fixed selective-COO segment metadata with shapes ``(2,)`` and
+        ``(1,)`` int32. Requires ``rebuild_flags``, caller-owned
+        ``neighbor_list``, and ``neighbor_list_shifts``.
     sorted_atom_index, morton_codes, sorted_pos_x, sorted_pos_y, sorted_pos_z, group_ctr_x, group_ctr_y, group_ctr_z, group_ext_x, group_ext_y, group_ext_z, num_tiles, tile_row_group, tile_col_group : torch.Tensor, optional
         Pre-allocated scratch buffers (shapes as returned by
         ``allocate_cluster_tile_list``).  All-or-nothing: either
@@ -1906,8 +2099,15 @@ def cluster_tile_neighbor_list(
           ``(*, vectors)`` appended when ``return_distances`` /
           ``return_vectors`` is True, and optional ``(*, pair_energies,
           pair_forces)`` when ``pair_fn`` is set. With ``cutoff2``, returns
-          the primary group followed by the secondary cutoff group.
-        - ``"coo"``: ``(neighbor_list, neighbor_ptr, neighbor_list_shifts)``.
+          the primary group followed by the secondary cutoff group. Selective
+          calls with ``return_state=True`` append ``(num_tiles,
+          tile_row_group, tile_col_group)``, yielding six tensors for one
+          cutoff or nine for two cutoffs.
+        - ``"coo"``: compact calls return ``(neighbor_list, neighbor_ptr,
+          neighbor_list_shifts)``. Selective calls return fixed-capacity
+          ``(neighbor_list, pair_offsets, pair_counts,
+          neighbor_list_shifts)`` and append tile state when
+          ``return_state=True``.
         - ``"tile"``: ``(num_tiles, tile_row_group, tile_col_group,
           sorted_atom_index, sorted_pos_x, sorted_pos_y, sorted_pos_z)``.
 
@@ -1949,6 +2149,14 @@ def cluster_tile_neighbor_list(
     )
     dual_cutoff = cutoff2 is not None
     selective = rebuild_flags is not None
+    is_compiling = torch.compiler.is_compiling()
+    eager_all_true = (
+        selective and not is_compiling and bool(rebuild_flags.flatten()[0].item())
+    )
+    if return_state and not selective:
+        raise ValueError("return_state=True requires rebuild_flags")
+    if not selective and (pair_offsets is not None or pair_counts is not None):
+        raise ValueError("pair_offsets and pair_counts require rebuild_flags")
     if has_pair_outputs and format == "tile":
         raise NotImplementedError(
             "Pair outputs (return_vectors / return_distances / pair_fn) "
@@ -1965,9 +2173,10 @@ def cluster_tile_neighbor_list(
                 "cluster_tile cutoff2 cannot be combined with pair outputs"
             )
     if selective:
-        if format != "matrix":
+        if format not in ("matrix", "coo"):
             raise ValueError(
-                "cluster_tile selective rebuild is supported only with format='matrix'"
+                "cluster_tile selective rebuild is supported only with "
+                "format='matrix' or format='coo'"
             )
         if has_pair_outputs:
             raise ValueError(
@@ -1977,10 +2186,24 @@ def cluster_tile_neighbor_list(
             "num_tiles": num_tiles,
             "tile_row_group": tile_row_group,
             "tile_col_group": tile_col_group,
-            "neighbor_matrix": neighbor_matrix,
-            "num_neighbors": num_neighbors,
-            "neighbor_matrix_shifts": neighbor_matrix_shifts,
         }
+        if format == "matrix":
+            required.update(
+                {
+                    "neighbor_matrix": neighbor_matrix,
+                    "num_neighbors": num_neighbors,
+                    "neighbor_matrix_shifts": neighbor_matrix_shifts,
+                }
+            )
+        else:
+            required.update(
+                {
+                    "neighbor_list": neighbor_list,
+                    "pair_offsets": pair_offsets,
+                    "pair_counts": pair_counts,
+                    "neighbor_list_shifts": neighbor_list_shifts,
+                }
+            )
         if dual_cutoff:
             required.update(
                 {
@@ -1990,7 +2213,8 @@ def cluster_tile_neighbor_list(
                 }
             )
         missing = [name for name, value in required.items() if value is None]
-        if missing:
+        bootstrap_coo = format == "coo" and eager_all_true and bool(return_state)
+        if missing and not bootstrap_coo:
             raise ValueError(
                 "rebuild_flags requires previous cluster_tile state: "
                 + ", ".join(missing)
@@ -1998,16 +2222,113 @@ def cluster_tile_neighbor_list(
     N = positions.shape[0]
     device = positions.device
     if max_neighbors is None:
-        max_neighbors = max(
-            estimate_max_neighbors(cutoff2 if cutoff2 is not None else cutoff),
-            TILE_GROUP_SIZE,
+        max_neighbors = (
+            int(neighbor_matrix.shape[1])
+            if format == "matrix" and neighbor_matrix is not None
+            else max(
+                estimate_max_neighbors(cutoff2 if cutoff2 is not None else cutoff),
+                TILE_GROUP_SIZE,
+            )
         )
+    elif (
+        format == "matrix"
+        and neighbor_matrix is not None
+        and (neighbor_matrix.shape != (N, int(max_neighbors)))
+    ):
+        raise ValueError("neighbor_matrix must have shape (N, max_neighbors)")
     if fill_value is None:
         fill_value = N
+    if selective and format == "matrix":
+        matrix_state = {
+            "neighbor_matrix": neighbor_matrix,
+            "num_neighbors": num_neighbors,
+            "neighbor_matrix_shifts": neighbor_matrix_shifts,
+            "num_tiles": num_tiles,
+            "tile_row_group": tile_row_group,
+            "tile_col_group": tile_col_group,
+        }
+        for name, value in matrix_state.items():
+            if value is None or value.device != device or value.dtype != torch.int32:
+                raise ValueError(f"{name} must be an int32 tensor on positions.device")
+        if (
+            num_neighbors.shape != (N,)
+            or neighbor_matrix_shifts.shape != (N, int(max_neighbors), 3)
+            or num_tiles.shape != (1,)
+            or tile_row_group.ndim != 1
+            or tile_col_group.shape != tile_row_group.shape
+        ):
+            raise ValueError("selective matrix state has an invalid shape")
+        if dual_cutoff and (
+            neighbor_matrix2 is None
+            or num_neighbors2 is None
+            or neighbor_matrix_shifts2 is None
+            or neighbor_matrix2.dtype != torch.int32
+            or num_neighbors2.dtype != torch.int32
+            or neighbor_matrix_shifts2.dtype != torch.int32
+            or neighbor_matrix2.device != device
+            or num_neighbors2.device != device
+            or neighbor_matrix_shifts2.device != device
+            or neighbor_matrix2.shape != (N, int(max_neighbors))
+            or num_neighbors2.shape != (N,)
+            or neighbor_matrix_shifts2.shape != (N, int(max_neighbors), 3)
+        ):
+            raise ValueError("selective secondary matrix state has an invalid shape")
+    if selective and format == "coo" and eager_all_true:
+        coo_state = (
+            neighbor_list,
+            pair_offsets,
+            pair_counts,
+            neighbor_list_shifts,
+            num_tiles,
+            tile_row_group,
+            tile_col_group,
+        )
+        if any(value is not None for value in coo_state) and any(
+            value is None for value in coo_state
+        ):
+            raise ValueError(
+                "selective COO bootstrap requires complete caller-owned state "
+                "or no persistent state"
+            )
+        if neighbor_list is None:
+            capacity = max_pairs if max_pairs is not None else N * max_neighbors
+            neighbor_list = torch.empty(
+                (2, capacity), dtype=torch.int32, device=positions.device
+            )
+            neighbor_list_shifts = torch.empty(
+                (capacity, 3), dtype=torch.int32, device=positions.device
+            )
+            pair_offsets = torch.tensor(
+                [0, capacity], dtype=torch.int32, device=positions.device
+            )
+            pair_counts = torch.zeros(1, dtype=torch.int32, device=positions.device)
+    if selective and format == "coo":
+        max_pairs_from_buffers = _validate_segmented_coo_state(
+            device=positions.device,
+            num_systems=1,
+            neighbor_list=neighbor_list,
+            neighbor_list_shifts=neighbor_list_shifts,
+            pair_offsets=pair_offsets,
+            pair_counts=pair_counts,
+            rebuild_flags=rebuild_flags,
+        )
+        if max_pairs is not None and int(max_pairs) != max_pairs_from_buffers:
+            raise ValueError("max_pairs must equal neighbor_list.shape[1]")
 
-    if selective and not bool(rebuild_flags.flatten()[0].item()):
-        if dual_cutoff:
-            return (
+    if (
+        selective
+        and not torch.compiler.is_compiling()
+        and not bool(rebuild_flags.flatten()[0].item())
+    ):
+        if format == "coo":
+            outputs = (
+                neighbor_list,
+                pair_offsets,
+                pair_counts,
+                neighbor_list_shifts,
+            )
+        elif dual_cutoff:
+            outputs = (
                 neighbor_matrix,
                 num_neighbors,
                 neighbor_matrix_shifts,
@@ -2015,7 +2336,11 @@ def cluster_tile_neighbor_list(
                 num_neighbors2,
                 neighbor_matrix_shifts2,
             )
-        return neighbor_matrix, num_neighbors, neighbor_matrix_shifts
+        else:
+            outputs = (neighbor_matrix, num_neighbors, neighbor_matrix_shifts)
+        if return_state:
+            return (*outputs, num_tiles, tile_row_group, tile_col_group)
+        return outputs
 
     # Autograd routing: when only distances/vectors are requested (no
     # pair_fn / energies / forces), route through the autograd primitive
@@ -2058,8 +2383,19 @@ def cluster_tile_neighbor_list(
     # Allocate scratch if caller didn't supply.  ``sorted_atom_index`` is the
     # all-or-nothing sentinel.
     if sorted_atom_index is None:
-        cell_volume = float(_cell_volume(cell))
-        max_tiles_per_group = estimate_max_tiles_per_group(N, build_cutoff, cell_volume)
+        previous_tile_state = (num_tiles, tile_row_group, tile_col_group)
+        if selective and torch.compiler.is_compiling():
+            if tile_row_group is None:
+                raise ValueError("selective compiled COO requires tile_row_group")
+            groups = max(1, (N + TILE_GROUP_SIZE - 1) // TILE_GROUP_SIZE)
+            max_tiles_per_group = max(1, int(tile_row_group.shape[0]) // groups)
+        else:
+            cell_volume = float(_cell_volume(cell))
+            max_tiles_per_group = estimate_max_tiles_per_group(
+                N,
+                build_cutoff,
+                cell_volume,
+            )
         (
             sorted_atom_index,
             morton_codes,
@@ -2081,6 +2417,14 @@ def cluster_tile_neighbor_list(
             dtype=positions.dtype,
             max_tiles_per_group=max_tiles_per_group,
         )
+        if selective:
+            previous_num_tiles, previous_row, previous_col = previous_tile_state
+            if previous_num_tiles is not None:
+                num_tiles = previous_num_tiles
+            if previous_row is not None:
+                tile_row_group = previous_row
+            if previous_col is not None:
+                tile_col_group = previous_col
     build_cluster_tile_list(
         positions,
         build_cutoff,
@@ -2099,6 +2443,7 @@ def cluster_tile_neighbor_list(
         num_tiles,
         tile_row_group,
         tile_col_group,
+        rebuild_flags=rebuild_flags if selective else None,
     )
 
     if format == "tile":
@@ -2118,7 +2463,9 @@ def cluster_tile_neighbor_list(
         )
 
     if format == "coo":
-        if max_pairs is None:
+        if selective:
+            max_pairs = int(neighbor_list.shape[1])
+        elif max_pairs is None:
             max_pairs = N * max_neighbors
         # ``query_cluster_tile_coo`` writes row-major (max_pairs, 2); we transpose to
         # package-canonical (2, num_pairs) on the way out.  Pre-allocation
@@ -2157,6 +2504,9 @@ def cluster_tile_neighbor_list(
             pair_counter,
             coo_buf,
             neighbor_list_shifts,
+            rebuild_flags=rebuild_flags,
+            pair_offsets=pair_offsets,
+            pair_counts=pair_counts,
             return_vectors=return_vectors,
             return_distances=return_distances,
             pair_fn=pair_fn,
@@ -2166,6 +2516,23 @@ def cluster_tile_neighbor_list(
             pair_energies=pair_energies,
             pair_forces=pair_forces,
         )
+        if selective:
+            if not torch.compiler.is_compiling():
+                pair_capacity = int(pair_offsets[1].item()) - int(
+                    pair_offsets[0].item()
+                )
+                pair_count = int(pair_counts[0].item())
+                if pair_count > pair_capacity:
+                    raise NeighborOverflowError(pair_capacity, pair_count)
+            outputs = (
+                neighbor_list,
+                pair_offsets,
+                pair_counts,
+                neighbor_list_shifts,
+            )
+            if return_state:
+                return (*outputs, num_tiles, tile_row_group, tile_col_group)
+            return outputs
         # Trim to actual pair count and rebuild CSR neighbor_ptr.
         # The ``.item()`` is the only sync; it's needed for the slice
         # anyway, so the bincount is on a CPU-known-length tensor.
@@ -2319,7 +2686,7 @@ def cluster_tile_neighbor_list(
             )
 
     if dual_cutoff:
-        return (
+        outputs = (
             neighbor_matrix,
             num_neighbors,
             neighbor_matrix_shifts,
@@ -2327,4 +2694,8 @@ def cluster_tile_neighbor_list(
             num_neighbors2,
             neighbor_matrix_shifts2,
         )
-    return neighbor_matrix, num_neighbors, neighbor_matrix_shifts
+    else:
+        outputs = (neighbor_matrix, num_neighbors, neighbor_matrix_shifts)
+    if return_state:
+        return (*outputs, num_tiles, tile_row_group, tile_col_group)
+    return outputs
diff --git a/nvalchemiops/torch/neighbors/neighbor_utils.py b/nvalchemiops/torch/neighbors/neighbor_utils.py
index a7d779d2..0b6551b4 100644
--- a/nvalchemiops/torch/neighbors/neighbor_utils.py
+++ b/nvalchemiops/torch/neighbors/neighbor_utils.py
@@ -62,6 +62,216 @@ def _validate_pair_params_present(
         raise ValueError("pair_params is required when pair_fn is provided")
 
 
+def _validate_segmented_coo_structure(
+    *,
+    device: torch.device,
+    num_systems: int,
+    neighbor_list: torch.Tensor,
+    neighbor_list_shifts: torch.Tensor,
+    pair_offsets: torch.Tensor,
+    pair_counts: torch.Tensor,
+    rebuild_flags: torch.Tensor | None,
+    tile_offsets: torch.Tensor | None = None,
+    tile_counts: torch.Tensor | None = None,
+    num_tiles: torch.Tensor | None = None,
+    tile_row_group: torch.Tensor | None = None,
+    tile_col_group: torch.Tensor | None = None,
+    tile_system: torch.Tensor | None = None,
+) -> int:
+    """Validate segmented COO shapes, dtypes, devices, and capacities."""
+    tensors = {
+        "neighbor_list": neighbor_list,
+        "neighbor_list_shifts": neighbor_list_shifts,
+        "pair_offsets": pair_offsets,
+        "pair_counts": pair_counts,
+    }
+    if rebuild_flags is not None:
+        tensors["rebuild_flags"] = rebuild_flags
+    optional_tensors = {
+        "tile_offsets": tile_offsets,
+        "tile_counts": tile_counts,
+        "num_tiles": num_tiles,
+        "tile_row_group": tile_row_group,
+        "tile_col_group": tile_col_group,
+        "tile_system": tile_system,
+    }
+    tensors.update(
+        {name: value for name, value in optional_tensors.items() if value is not None}
+    )
+    for name, value in tensors.items():
+        if value.device != device:
+            raise ValueError(f"{name} must match positions.device")
+
+    if (
+        neighbor_list.dtype != torch.int32
+        or neighbor_list.ndim != 2
+        or neighbor_list.shape[0] != 2
+    ):
+        raise ValueError("neighbor_list must have shape (2, capacity) and dtype int32")
+    capacity = int(neighbor_list.shape[1])
+    if neighbor_list_shifts.dtype != torch.int32 or neighbor_list_shifts.shape != (
+        capacity,
+        3,
+    ):
+        raise ValueError(
+            "neighbor_list_shifts must have shape (neighbor_list capacity, 3) "
+            "and dtype int32"
+        )
+    if pair_offsets.dtype != torch.int32 or pair_offsets.shape != (num_systems + 1,):
+        raise ValueError(
+            "pair_offsets must have shape (num_systems + 1,) and dtype int32"
+        )
+    if pair_counts.dtype != torch.int32 or pair_counts.shape != (num_systems,):
+        raise ValueError("pair_counts must have shape (num_systems,) and dtype int32")
+    if rebuild_flags is not None and (
+        rebuild_flags.dtype != torch.bool or rebuild_flags.shape != (num_systems,)
+    ):
+        raise ValueError("rebuild_flags must have shape (num_systems,) and dtype bool")
+
+    tile_values = (tile_offsets, tile_counts, num_tiles, tile_row_group, tile_col_group)
+    if any(value is not None for value in tile_values):
+        if any(value is None for value in tile_values):
+            raise ValueError("segmented COO tile state must be supplied completely")
+        if (
+            tile_offsets.dtype != torch.int32
+            or tile_offsets.shape != (num_systems + 1,)
+            or tile_counts.dtype != torch.int32
+            or tile_counts.shape != (num_systems,)
+            or num_tiles.dtype != torch.int32
+            or num_tiles.shape != (1,)
+        ):
+            raise ValueError(
+                "segmented COO tile metadata has an invalid shape or dtype"
+            )
+        if (
+            tile_row_group.dtype != torch.int32
+            or tile_col_group.dtype != torch.int32
+            or tile_row_group.ndim != 1
+            or tile_col_group.ndim != 1
+            or tile_row_group.shape != tile_col_group.shape
+        ):
+            raise ValueError(
+                "tile row and column buffers must be matching 1D int32 arrays"
+            )
+        if tile_system is not None and (
+            tile_system.dtype != torch.int32
+            or tile_system.ndim != 1
+            or tile_system.shape != tile_row_group.shape
+        ):
+            raise ValueError("tile_system must match tile row and column buffer shapes")
+
+    return capacity
+
+
+def _validate_segmented_coo_values(
+    *,
+    neighbor_list: torch.Tensor,
+    pair_offsets: torch.Tensor,
+    pair_counts: torch.Tensor,
+    tile_offsets: torch.Tensor | None = None,
+    tile_counts: torch.Tensor | None = None,
+    tile_row_group: torch.Tensor | None = None,
+) -> None:
+    """Eagerly validate segmented COO metadata values before mutation."""
+    capacity = int(neighbor_list.shape[1])
+    pair_offsets_host = pair_offsets.cpu()
+    pair_counts_host = pair_counts.cpu()
+    if int(pair_offsets_host[0]) != 0:
+        raise ValueError("pair_offsets must start at zero")
+    if bool((pair_offsets_host[1:] < pair_offsets_host[:-1]).any()):
+        raise ValueError("pair_offsets must be nondecreasing")
+    if int(pair_offsets_host[-1]) != capacity:
+        raise ValueError("pair_offsets final value must equal neighbor_list capacity")
+    pair_capacities = pair_offsets_host[1:] - pair_offsets_host[:-1]
+    if bool((pair_counts_host < 0).any()) or bool(
+        (pair_counts_host > pair_capacities).any()
+    ):
+        raise ValueError("pair_counts must lie within their pair_offsets segments")
+
+    if tile_offsets is not None:
+        tile_offsets_host = tile_offsets.cpu()
+        tile_counts_host = tile_counts.cpu()
+        if int(tile_offsets_host[0]) != 0:
+            raise ValueError("tile_offsets must start at zero")
+        if bool((tile_offsets_host[1:] < tile_offsets_host[:-1]).any()):
+            raise ValueError("tile_offsets must be nondecreasing")
+        if int(tile_offsets_host[-1]) > int(tile_row_group.shape[0]):
+            raise ValueError("tile_offsets exceed the physical tile-buffer capacity")
+        tile_capacities = tile_offsets_host[1:] - tile_offsets_host[:-1]
+        if bool((tile_counts_host < 0).any()) or bool(
+            (tile_counts_host > tile_capacities).any()
+        ):
+            raise ValueError("tile_counts must lie within their tile_offsets segments")
+
+
+def _validate_segmented_coo_state(
+    *,
+    device: torch.device,
+    num_systems: int,
+    neighbor_list: torch.Tensor,
+    neighbor_list_shifts: torch.Tensor,
+    pair_offsets: torch.Tensor,
+    pair_counts: torch.Tensor,
+    rebuild_flags: torch.Tensor | None,
+    tile_offsets: torch.Tensor | None = None,
+    tile_counts: torch.Tensor | None = None,
+    num_tiles: torch.Tensor | None = None,
+    tile_row_group: torch.Tensor | None = None,
+    tile_col_group: torch.Tensor | None = None,
+    tile_system: torch.Tensor | None = None,
+) -> int:
+    """Validate fixed-capacity segmented COO state before a kernel launch."""
+    capacity = _validate_segmented_coo_structure(
+        device=device,
+        num_systems=num_systems,
+        neighbor_list=neighbor_list,
+        neighbor_list_shifts=neighbor_list_shifts,
+        pair_offsets=pair_offsets,
+        pair_counts=pair_counts,
+        rebuild_flags=rebuild_flags,
+        tile_offsets=tile_offsets,
+        tile_counts=tile_counts,
+        num_tiles=num_tiles,
+        tile_row_group=tile_row_group,
+        tile_col_group=tile_col_group,
+        tile_system=tile_system,
+    )
+    if not torch.compiler.is_compiling():
+        _validate_segmented_coo_values(
+            neighbor_list=neighbor_list,
+            pair_offsets=pair_offsets,
+            pair_counts=pair_counts,
+            tile_offsets=tile_offsets,
+            tile_counts=tile_counts,
+            tile_row_group=tile_row_group,
+        )
+    return capacity
+
+
+def _normalize_compiled_single_segment_coo_count(
+    *,
+    pair_offsets: torch.Tensor,
+    pair_counts: torch.Tensor,
+    physical_capacity: int,
+) -> None:
+    """Fail closed for malformed compiled single-segment COO metadata.
+
+    This compiled-path helper uses only device-side int32 operations. A false
+    rebuild flag returns before the Warp query validates metadata, while an
+    overflowed attempted count is not final until every query block completes.
+    It therefore validates the exact fixed interval and clamps the final count
+    here. This is deliberately not a generic batched normalizer.
+    """
+    offsets_valid = (pair_offsets[0] == 0) & (pair_offsets[1] == physical_capacity)
+    clamped_count = torch.clamp(pair_counts, min=0, max=physical_capacity)
+    normalized_counts = torch.where(
+        offsets_valid,
+        clamped_count,
+        torch.zeros_like(pair_counts),
+    )
+    pair_counts.copy_(normalized_counts)
+
+
 def compute_naive_num_shifts(
     cell: torch.Tensor,
     cutoff: float,
diff --git a/test/neighbors/bindings/jax/test_batch_cluster_tile.py b/test/neighbors/bindings/jax/test_batch_cluster_tile.py
index 763de87a..2a21e6f9 100644
--- a/test/neighbors/bindings/jax/test_batch_cluster_tile.py
+++ b/test/neighbors/bindings/jax/test_batch_cluster_tile.py
@@ -1103,3 +1103,235 @@ def test_rebuild_flags_require_previous_state(self):
                 max_neighbors=16,
                 rebuild_flags=jnp.array([True], dtype=jnp.bool_),
             )
+
+
+class TestJaxBatchClusterTileEmptySelective:
+    """Selective zero-atom returns keep their fixed tuple layout."""
+
+    @staticmethod
+    def _tile_state():
+        """Return distinctive batched tile state buffers."""
+        return {
+            "tile_offsets": jnp.array([0, 4], dtype=jnp.int32),
+            "previous_tile_counts": jnp.array([3], dtype=jnp.int32),
+            "previous_num_tiles": jnp.array([3], dtype=jnp.int32),
+            "previous_tile_row_group": jnp.arange(4, dtype=jnp.int32),
+            "previous_tile_col_group": jnp.arange(4, dtype=jnp.int32),
+            "previous_tile_system": jnp.zeros(4, dtype=jnp.int32),
+        }
+
+    @pytest.mark.parametrize("rebuild_flag", [False, True])
+    def test_empty_selective_coo_returns_segmented_state(self, rebuild_flag):
+        """Zero atoms retain the ten-array batched COO contract."""
+        batch_ptr = jnp.array([0, 0], dtype=jnp.int32)
+        pair_offsets = jnp.array([0, 5], dtype=jnp.int32)
+        previous_pair_counts = jnp.array([4], dtype=jnp.int32)
+        previous_neighbor_list = jnp.full((2, 5), 7, dtype=jnp.int32)
+        previous_neighbor_list_shifts = jnp.full((5, 3), 7, dtype=jnp.int32)
+        tile_state = self._tile_state()
+
+        out = batch_cluster_tile_neighbor_list(
+            jnp.empty((0, 3), dtype=jnp.float32),
+            1.0,
+            jnp.eye(3, dtype=jnp.float32)[jnp.newaxis, :, :],
+            batch_ptr,
+            max_neighbors=8,
+            format="coo",
+            rebuild_flags=jnp.array([rebuild_flag], dtype=jnp.bool_),
+            pair_offsets=pair_offsets,
+            previous_pair_counts=previous_pair_counts,
+            previous_neighbor_list=previous_neighbor_list,
+            previous_neighbor_list_shifts=previous_neighbor_list_shifts,
+            **tile_state,
+        )
+
+        (
+            neighbor_list,
+            offsets,
+            pair_counts,
+            shifts,
+            returned_tile_offsets,
+            tile_counts,
+            num_tiles,
+            row,
+            col,
+            system,
+        ) = out
+        assert neighbor_list.shape == (2, 5)
+        assert shifts.shape == (5, 3)
+        np.testing.assert_array_equal(np.asarray(offsets), np.asarray(pair_offsets))
+        np.testing.assert_array_equal(
+            np.asarray(pair_counts),
+            np.array([0 if rebuild_flag else 4], dtype=np.int32),
+        )
+        np.testing.assert_array_equal(
+            np.asarray(neighbor_list), np.asarray(previous_neighbor_list)
+        )
+        np.testing.assert_array_equal(
+            np.asarray(shifts),
+            np.asarray(previous_neighbor_list_shifts),
+        )
+        np.testing.assert_array_equal(
+            np.asarray(returned_tile_offsets), np.asarray(tile_state["tile_offsets"])
+        )
+        np.testing.assert_array_equal(
+            np.asarray(tile_counts),
+            np.array([0 if rebuild_flag else 3], dtype=np.int32),
+        )
+        np.testing.assert_array_equal(
+            np.asarray(num_tiles), np.asarray(tile_state["previous_num_tiles"])
+        )
+        np.testing.assert_array_equal(
+            np.asarray(row), np.asarray(tile_state["previous_tile_row_group"])
+        )
+        np.testing.assert_array_equal(
+            np.asarray(col), np.asarray(tile_state["previous_tile_col_group"])
+        )
+        np.testing.assert_array_equal(
+            np.asarray(system), np.asarray(tile_state["previous_tile_system"])
+        )
+
+    @pytest.mark.parametrize("rebuild_flag", [False, True])
+    def test_empty_selective_matrix_returns_state(self, rebuild_flag):
+        """Zero atoms retain the nine-array selective matrix contract."""
+        tile_state = self._tile_state()
+        out = batch_cluster_tile_neighbor_list(
+            jnp.empty((0, 3), dtype=jnp.float32),
+            1.0,
+            jnp.eye(3, dtype=jnp.float32)[jnp.newaxis, :, :],
+            jnp.array([0, 0], dtype=jnp.int32),
+            max_neighbors=8,
+            rebuild_flags=jnp.array([rebuild_flag], dtype=jnp.bool_),
+            previous_neighbor_matrix=jnp.empty((0, 8), dtype=jnp.int32),
+            previous_num_neighbors=jnp.empty(0, dtype=jnp.int32),
+            previous_neighbor_matrix_shifts=jnp.empty((0, 8, 3), dtype=jnp.int32),
+            **tile_state,
+        )
+
+        assert len(out) == 9
+        assert tuple(value.shape for value in out[:3]) == (
+            (0, 8),
+            (0,),
+            (0, 8, 3),
+        )
+        for returned, name in zip(out[3:], tile_state):
+            expected = (
+                np.zeros(1, dtype=np.int32)
+                if rebuild_flag and name == "previous_tile_counts"
+                else np.asarray(tile_state[name])
+            )
+            np.testing.assert_array_equal(
+                np.asarray(returned),
+                expected,
+            )
+
+    def test_empty_two_system_selective_matrix_returns_state(self):
+        """Two empty systems retain per-system offsets, counts, and flags."""
+        tile_state = {
+            "tile_offsets": jnp.array([0, 2, 4], dtype=jnp.int32),
+            "previous_tile_counts": jnp.array([1, 2], dtype=jnp.int32),
+            "previous_num_tiles": jnp.array([3], dtype=jnp.int32),
+            "previous_tile_row_group": jnp.arange(4, dtype=jnp.int32),
+            "previous_tile_col_group": jnp.arange(4, dtype=jnp.int32),
+            "previous_tile_system": jnp.array([0, 0, 1, 1], dtype=jnp.int32),
+        }
+        out = batch_cluster_tile_neighbor_list(
+            jnp.empty((0, 3), dtype=jnp.float32),
+            1.0,
+            jnp.repeat(jnp.eye(3, dtype=jnp.float32)[jnp.newaxis, :, :], 2, axis=0),
+            jnp.array([0, 0, 0], dtype=jnp.int32),
+            max_neighbors=8,
+            rebuild_flags=jnp.array([False, True], dtype=jnp.bool_),
+            previous_neighbor_matrix=jnp.empty((0, 8), dtype=jnp.int32),
+            previous_num_neighbors=jnp.empty(0, dtype=jnp.int32),
+            previous_neighbor_matrix_shifts=jnp.empty((0, 8, 3), dtype=jnp.int32),
+            **tile_state,
+        )
+
+        assert len(out) == 9
+        expected_state = (
+            tile_state["tile_offsets"],
+            jnp.array([1, 0], dtype=jnp.int32),
+            tile_state["previous_num_tiles"],
+            tile_state["previous_tile_row_group"],
+            tile_state["previous_tile_col_group"],
+            tile_state["previous_tile_system"],
+        )
+        for returned, expected in zip(out[3:], expected_state):
+            np.testing.assert_array_equal(
+                np.asarray(returned),
+                np.asarray(expected),
+            )
+
+    def test_empty_selective_dual_cutoff_matrix_returns_both_outputs(self):
+        """Zero atoms retain both matrix triples plus batched selective state."""
+        tile_state = self._tile_state()
+        matrix = jnp.empty((0, 8), dtype=jnp.int32)
+        counts = jnp.empty(0, dtype=jnp.int32)
+        shifts = jnp.empty((0, 8, 3), dtype=jnp.int32)
+
+        out = batch_cluster_tile_neighbor_list(
+            jnp.empty((0, 3), dtype=jnp.float32),
+            1.0,
+            jnp.eye(3, dtype=jnp.float32)[jnp.newaxis, :, :],
+            jnp.array([0, 0], dtype=jnp.int32),
+            cutoff2=2.0,
+            max_neighbors=8,
+            rebuild_flags=jnp.array([True], dtype=jnp.bool_),
+            previous_neighbor_matrix=matrix,
+            previous_num_neighbors=counts,
+            previous_neighbor_matrix_shifts=shifts,
+            previous_neighbor_matrix2=matrix,
+            previous_num_neighbors2=counts,
+            previous_neighbor_matrix_shifts2=shifts,
+            **tile_state,
+        )
+
+        assert len(out) == 12
+        assert tuple(value.shape for value in out[:6]) == (
+            (0, 8),
+            (0,),
+            (0, 8, 3),
+            (0, 8),
+            (0,),
+            (0, 8, 3),
+        )
+        np.testing.assert_array_equal(
+            np.asarray(out[6]), np.asarray(tile_state["tile_offsets"])
+        )
+
+    @pytest.mark.parametrize("rebuild_flag, expected_count", [(False, 4), (True, 0)])
+    def test_jit_empty_selective_coo_keeps_fixed_arity(
+        self,
+        rebuild_flag,
+        expected_count,
+    ):
+        """JIT preserves the ten-array zero-atom selective COO contract."""
+        tile_state = self._tile_state()
+
+        @jax.jit
+        def build(positions, rebuild_flags):
+            return batch_cluster_tile_neighbor_list(
+                positions,
+                1.0,
+                jnp.eye(3, dtype=jnp.float32)[jnp.newaxis, :, :],
+                jnp.array([0, 0], dtype=jnp.int32),
+                max_neighbors=8,
+                format="coo",
+                rebuild_flags=rebuild_flags,
+                pair_offsets=jnp.array([0, 5], dtype=jnp.int32),
+                previous_pair_counts=jnp.array([4], dtype=jnp.int32),
+                previous_neighbor_list=jnp.full((2, 5), 7, dtype=jnp.int32),
+                previous_neighbor_list_shifts=jnp.full((5, 3), 7, dtype=jnp.int32),
+                **tile_state,
+            )
+
+        out = build(
+            jnp.empty((0, 3), dtype=jnp.float32),
+            jnp.array([rebuild_flag], dtype=jnp.bool_),
+        )
+
+        assert len(out) == 10
+        assert out[0].shape == (2, 5)
+        assert out[2].shape == (1,)
+        assert int(out[2][0]) == expected_count
diff --git a/test/neighbors/bindings/jax/test_cluster_tile.py b/test/neighbors/bindings/jax/test_cluster_tile.py
index c33ef9f0..6b30cb52 100644
--- a/test/neighbors/bindings/jax/test_cluster_tile.py
+++ b/test/neighbors/bindings/jax/test_cluster_tile.py
@@ -31,6 +31,7 @@
     build_cluster_tile_list,
     cluster_tile_neighbor_list,
     estimate_cluster_tile_list_sizes,
+    query_cluster_tile_coo,
 )
 
 from .conftest import requires_gpu
@@ -866,6 +867,44 @@ def test_rebuild_flags_coo_true_writes_segment_counts(self):
         assert int(pair_counts[0]) <= int(pair_offsets[1] - pair_offsets[0])
         assert int(nt2[0]) > 0
 
+    def test_rebuild_flags_coo_clamps_valid_overflow_to_physical_capacity(self):
+        """Valid fixed metadata cannot report more pairs than it can store."""
+        positions = jnp.zeros((64, 3), dtype=jnp.float32)
+        cell = _orthorhombic_cell(6.0)
+        num_tiles, tile_row_group, tile_col_group, *_ = cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell,
+            format="tile",
+        )
+        capacity = 64
+        out = cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell,
+            max_neighbors=64,
+            format="coo",
+            rebuild_flags=jnp.array([True], dtype=jnp.bool_),
+            previous_num_tiles=num_tiles,
+            previous_tile_row_group=tile_row_group,
+            previous_tile_col_group=tile_col_group,
+            pair_offsets=jnp.array([0, capacity], dtype=jnp.int32),
+            previous_pair_counts=jnp.zeros(1, dtype=jnp.int32),
+            previous_neighbor_list=jnp.full((2, capacity), -77, dtype=jnp.int32),
+            previous_neighbor_list_shifts=jnp.full(
+                (capacity, 3),
+                -77,
+                dtype=jnp.int32,
+            ),
+        )
+
+        np.testing.assert_array_equal(
+            np.asarray(out[2]),
+            np.array([capacity], dtype=np.int32),
+        )
+        assert np.all(np.asarray(out[0]) != -77)
+        assert np.all(np.asarray(out[3]) != -77)
+
     def test_rebuild_flags_coo_requires_segment_buffers(self):
         positions = jnp.zeros((32, 3), dtype=jnp.float32)
         cell = _orthorhombic_cell(4.0)
@@ -888,6 +927,358 @@ def test_rebuild_flags_coo_requires_segment_buffers(self):
                 previous_tile_col_group=tile_col_group,
             )
 
+    def test_rebuild_flags_coo_requires_one_single_system_segment(self):
+        """Single-system segmented COO rejects unused extra metadata segments."""
+        positions = jnp.zeros((32, 3), dtype=jnp.float32)
+        cell = _orthorhombic_cell(6.0)
+        num_tiles, tile_row_group, tile_col_group, *_ = cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell,
+            format="tile",
+        )
+
+        with pytest.raises(ValueError, match=r"pair_offsets must have shape \(2,\)"):
+            cluster_tile_neighbor_list(
+                positions,
+                2.0,
+                cell,
+                max_neighbors=64,
+                format="coo",
+                rebuild_flags=jnp.array([False], dtype=jnp.bool_),
+                previous_num_tiles=num_tiles,
+                previous_tile_row_group=tile_row_group,
+                previous_tile_col_group=tile_col_group,
+                pair_offsets=jnp.array([0, 4, 8], dtype=jnp.int32),
+                previous_pair_counts=jnp.array([0, 0], dtype=jnp.int32),
+                previous_neighbor_list=jnp.full((2, 8), -77, dtype=jnp.int32),
+                previous_neighbor_list_shifts=jnp.full((8, 3), -77, dtype=jnp.int32),
+            )
+
+    @pytest.mark.parametrize(
+        ("pair_offsets_values", "rebuild_flag"),
+        [
+            ((64, 0), True),
+            ((0, 32), True),
+            ((1, 64), True),
+            ((64, 0), False),
+            ((0, 32), False),
+            ((1, 64), False),
+        ],
+        ids=[
+            "reversed-true",
+            "undersized-true",
+            "nonzero-start-true",
+            "reversed-false",
+            "undersized-false",
+            "nonzero-start-false",
+        ],
+    )
+    def test_jit_selective_coo_invalid_offsets_fail_closed(
+        self,
+        pair_offsets_values,
+        rebuild_flag,
+    ):
+        """JIT returns no active entries for malformed single-system metadata."""
+        positions = jnp.zeros((32, 3), dtype=jnp.float32)
+        cell = _orthorhombic_cell(6.0)
+        num_tiles, tile_row_group, tile_col_group, *_ = cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell,
+            format="tile",
+        )
+        neighbor_list = jnp.full((2, 64), -77, dtype=jnp.int32)
+        neighbor_list_shifts = jnp.full((64, 3), -77, dtype=jnp.int32)
+
+        @jax.jit
+        def run(pair_offsets, pair_counts, runtime_rebuild_flags):
+            return cluster_tile_neighbor_list(
+                positions,
+                2.0,
+                cell,
+                max_neighbors=64,
+                format="coo",
+                rebuild_flags=runtime_rebuild_flags,
+                previous_num_tiles=num_tiles,
+                previous_tile_row_group=tile_row_group,
+                previous_tile_col_group=tile_col_group,
+                pair_offsets=pair_offsets,
+                previous_pair_counts=pair_counts,
+                previous_neighbor_list=neighbor_list,
+                previous_neighbor_list_shifts=neighbor_list_shifts,
+            )
+
+        out = run(
+            jnp.array(pair_offsets_values, dtype=jnp.int32),
+            jnp.array([0 if rebuild_flag else 7], dtype=jnp.int32),
+            jnp.array([rebuild_flag], dtype=jnp.bool_),
+        )
+
+        np.testing.assert_array_equal(np.asarray(out[2]), np.zeros(1, dtype=np.int32))
+        np.testing.assert_array_equal(np.asarray(out[0]), np.asarray(neighbor_list))
+        np.testing.assert_array_equal(
+            np.asarray(out[3]),
+            np.asarray(neighbor_list_shifts),
+        )
+
+    def test_jit_query_coo_uses_fixed_buffer_capacity(self):
+        """Fixed buffers keep a traced max_pairs argument out of sizing logic."""
+        positions = jnp.zeros((32, 3), dtype=jnp.float32)
+        cell = _orthorhombic_cell(6.0)
+        state = build_cluster_tile_list(positions, 2.0, cell)
+        neighbor_list = jnp.full((2, 64), -77, dtype=jnp.int32)
+        neighbor_list_shifts = jnp.full((64, 3), -77, dtype=jnp.int32)
+
+        @jax.jit
+        def run(runtime_max_pairs):
+            return query_cluster_tile_coo(
+                state[0],
+                state[2],
+                state[3],
+                state[4],
+                state[11],
+                state[12],
+                state[13],
+                cell,
+                2.0,
+                32,
+                runtime_max_pairs,
+                rebuild_flags=jnp.array([False], dtype=jnp.bool_),
+                pair_offsets=jnp.array([0, 64], dtype=jnp.int32),
+                pair_counts=jnp.array([7], dtype=jnp.int32),
+                neighbor_list=neighbor_list,
+                neighbor_list_shifts=neighbor_list_shifts,
+            )
+
+        out = run(jnp.array(1, dtype=jnp.int32))
+
+        np.testing.assert_array_equal(np.asarray(out[2]), np.array([7], dtype=np.int32))
+        np.testing.assert_array_equal(np.asarray(out[0]), np.asarray(neighbor_list))
+        np.testing.assert_array_equal(
+            np.asarray(out[3]),
+            np.asarray(neighbor_list_shifts),
+        )
+
+
+class TestJaxClusterTileEmptySelective:
+    """Selective zero-atom returns keep their fixed tuple layout."""
+
+    @staticmethod
+    def _tile_state():
+        """Return distinctive single-system tile state buffers."""
+        return {
+            "previous_num_tiles": jnp.array([3], dtype=jnp.int32),
+            "previous_tile_row_group": jnp.arange(4, dtype=jnp.int32),
+            "previous_tile_col_group": jnp.arange(4, dtype=jnp.int32),
+        }
+
+    @pytest.mark.parametrize("rebuild_flag", [False, True])
+    def test_empty_selective_coo_returns_segmented_state(self, rebuild_flag):
+        """Zero atoms retain the seven-array single-system COO contract."""
+        pair_offsets = jnp.array([0, 5], dtype=jnp.int32)
+        previous_pair_counts = jnp.array([4], dtype=jnp.int32)
+        previous_neighbor_list = jnp.full((2, 5), 7, dtype=jnp.int32)
+        previous_neighbor_list_shifts = jnp.full((5, 3), 7, dtype=jnp.int32)
+        tile_state = self._tile_state()
+
+        out = cluster_tile_neighbor_list(
+            jnp.empty((0, 3), dtype=jnp.float32),
+            1.0,
+            _orthorhombic_cell(4.0),
+            max_neighbors=8,
+            format="coo",
+            rebuild_flags=jnp.array([rebuild_flag], dtype=jnp.bool_),
+            pair_offsets=pair_offsets,
+            previous_pair_counts=previous_pair_counts,
+            previous_neighbor_list=previous_neighbor_list,
+            previous_neighbor_list_shifts=previous_neighbor_list_shifts,
+            **tile_state,
+        )
+
+        neighbor_list, offsets, pair_counts, shifts, num_tiles, row, col = out
+        assert neighbor_list.shape == (2, 5)
+        assert shifts.shape == (5, 3)
+        np.testing.assert_array_equal(np.asarray(offsets), np.asarray(pair_offsets))
+        np.testing.assert_array_equal(
+            np.asarray(pair_counts),
+            np.array([0 if rebuild_flag else 4], dtype=np.int32),
+        )
+        np.testing.assert_array_equal(
+            np.asarray(neighbor_list), np.asarray(previous_neighbor_list)
+        )
+        np.testing.assert_array_equal(
+            np.asarray(shifts),
+            np.asarray(previous_neighbor_list_shifts),
+        )
+        np.testing.assert_array_equal(
+            np.asarray(num_tiles),
+            np.array([0 if rebuild_flag else 3], dtype=np.int32),
+        )
+        np.testing.assert_array_equal(
+            np.asarray(row), np.asarray(tile_state["previous_tile_row_group"])
+        )
+        np.testing.assert_array_equal(
+            np.asarray(col), np.asarray(tile_state["previous_tile_col_group"])
+        )
+
+    def test_empty_selective_coo_invalid_offsets_fail_closed(self):
+        """Malformed metadata zeros an inactive empty-system count."""
+        pair_offsets = jnp.array([0, 3], dtype=jnp.int32)
+        previous_pair_counts = jnp.array([4], dtype=jnp.int32)
+        previous_neighbor_list = jnp.full((2, 5), 7, dtype=jnp.int32)
+        previous_neighbor_list_shifts = jnp.full((5, 3), 7, dtype=jnp.int32)
+
+        out = cluster_tile_neighbor_list(
+            jnp.empty((0, 3), dtype=jnp.float32),
+            1.0,
+            _orthorhombic_cell(4.0),
+            max_neighbors=8,
+            format="coo",
+            rebuild_flags=jnp.array([False], dtype=jnp.bool_),
+            pair_offsets=pair_offsets,
+            previous_pair_counts=previous_pair_counts,
+            previous_neighbor_list=previous_neighbor_list,
+            previous_neighbor_list_shifts=previous_neighbor_list_shifts,
+            **self._tile_state(),
+        )
+
+        np.testing.assert_array_equal(np.asarray(out[2]), np.zeros(1, dtype=np.int32))
+        np.testing.assert_array_equal(
+            np.asarray(out[0]), np.asarray(previous_neighbor_list)
+        )
+        np.testing.assert_array_equal(
+            np.asarray(out[3]),
+            np.asarray(previous_neighbor_list_shifts),
+        )
+
+    def test_empty_selective_coo_rejects_mismatched_fixed_buffer_shapes(self):
+        """Fixed COO and shift buffers must describe one shared capacity."""
+        with pytest.raises(ValueError, match="neighbor_list_shifts"):
+            cluster_tile_neighbor_list(
+                jnp.empty((0, 3), dtype=jnp.float32),
+                1.0,
+                _orthorhombic_cell(4.0),
+                max_neighbors=8,
+                format="coo",
+                rebuild_flags=jnp.array([False], dtype=jnp.bool_),
+                pair_offsets=jnp.array([0, 5], dtype=jnp.int32),
+                previous_pair_counts=jnp.array([4], dtype=jnp.int32),
+                previous_neighbor_list=jnp.full((2, 5), 7, dtype=jnp.int32),
+                previous_neighbor_list_shifts=jnp.full((4, 3), 7, dtype=jnp.int32),
+                **self._tile_state(),
+            )
+
+    @pytest.mark.parametrize("rebuild_flag", [False, True])
+    def test_empty_selective_matrix_returns_state(self, rebuild_flag):
+        """Zero atoms retain the six-array selective matrix contract."""
+        tile_state = self._tile_state()
+        out = cluster_tile_neighbor_list(
+            jnp.empty((0, 3), dtype=jnp.float32),
+            1.0,
+            _orthorhombic_cell(4.0),
+            max_neighbors=8,
+            rebuild_flags=jnp.array([rebuild_flag], dtype=jnp.bool_),
+            previous_neighbor_matrix=jnp.empty((0, 8), dtype=jnp.int32),
+            previous_num_neighbors=jnp.empty(0, dtype=jnp.int32),
+            previous_neighbor_matrix_shifts=jnp.empty((0, 8, 3), dtype=jnp.int32),
+            **tile_state,
+        )
+
+        assert len(out) == 6
+        assert out[0].shape == (0, 8)
+        assert out[1].shape == (0,)
+        assert out[2].shape == (0, 8, 3)
+        for returned, name in zip(
+            out[3:],
+            (
+                "previous_num_tiles",
+                "previous_tile_row_group",
+                "previous_tile_col_group",
+            ),
+        ):
+            expected = (
+                np.zeros(1, dtype=np.int32)
+                if rebuild_flag and name == "previous_num_tiles"
+                else np.asarray(tile_state[name])
+            )
+            np.testing.assert_array_equal(
+                np.asarray(returned),
+                expected,
+            )
+
+    def test_empty_selective_dual_cutoff_matrix_returns_both_outputs(self):
+        """Zero atoms retain both matrix triples plus selective tile state."""
+        tile_state = self._tile_state()
+        matrix = jnp.empty((0, 8), dtype=jnp.int32)
+        counts = jnp.empty(0, dtype=jnp.int32)
+        shifts = jnp.empty((0, 8, 3), dtype=jnp.int32)
+
+        out = cluster_tile_neighbor_list(
+            jnp.empty((0, 3), dtype=jnp.float32),
+            1.0,
+            _orthorhombic_cell(4.0),
+            cutoff2=2.0,
+            max_neighbors=8,
+            rebuild_flags=jnp.array([True], dtype=jnp.bool_),
+            previous_neighbor_matrix=matrix,
+            previous_num_neighbors=counts,
+            previous_neighbor_matrix_shifts=shifts,
+            previous_neighbor_matrix2=matrix,
+            previous_num_neighbors2=counts,
+            previous_neighbor_matrix_shifts2=shifts,
+            **tile_state,
+        )
+
+        assert len(out) == 9
+        assert tuple(value.shape for value in out[:6]) == (
+            (0, 8),
+            (0,),
+            (0, 8, 3),
+            (0, 8),
+            (0,),
+            (0, 8, 3),
+        )
+        np.testing.assert_array_equal(
+            np.asarray(out[6]),
+            np.zeros(1, dtype=np.int32),
+        )
+
+    @pytest.mark.parametrize("rebuild_flag, expected_count", [(False, 4), (True, 0)])
+    def test_jit_empty_selective_coo_keeps_fixed_arity(
+        self,
+        rebuild_flag,
+        expected_count,
+    ):
+        """JIT preserves the seven-array zero-atom selective COO contract."""
+        tile_state = self._tile_state()
+
+        @jax.jit
+        def build(positions, rebuild_flags):
+            return cluster_tile_neighbor_list(
+                positions,
+                1.0,
+                _orthorhombic_cell(4.0),
+                max_neighbors=8,
+                format="coo",
+                rebuild_flags=rebuild_flags,
+                pair_offsets=jnp.array([0, 5], dtype=jnp.int32),
+                previous_pair_counts=jnp.array([4], dtype=jnp.int32),
+                previous_neighbor_list=jnp.full((2, 5), 7, dtype=jnp.int32),
+                previous_neighbor_list_shifts=jnp.full((5, 3), 7, dtype=jnp.int32),
+                **tile_state,
+            )
+
+        out = build(
+            jnp.empty((0, 3), dtype=jnp.float32),
+            jnp.array([rebuild_flag], dtype=jnp.bool_),
+        )
+
+        assert len(out) == 7
+        assert out[0].shape == (2, 5)
+        assert out[2].shape == (1,)
+        assert int(out[2][0]) == expected_count
+
 
 class TestJaxClusterTileNeighborListDispatcher:
     """Unified JAX neighbor_list routes cluster-tile kwargs explicitly."""
diff --git a/test/neighbors/bindings/jax/test_naive_dual_cutoff.py b/test/neighbors/bindings/jax/test_naive_dual_cutoff.py
index f3df03fb..98479ad1 100644
--- a/test/neighbors/bindings/jax/test_naive_dual_cutoff.py
+++ b/test/neighbors/bindings/jax/test_naive_dual_cutoff.py
@@ -21,9 +21,13 @@
 
 import jax
 import jax.numpy as jnp
+import numpy as np
 import pytest
 
-from nvalchemiops.jax.neighbors import naive_neighbor_list_dual_cutoff
+from nvalchemiops.jax.neighbors import (
+    naive_neighbor_list,
+    naive_neighbor_list_dual_cutoff,
+)
 
 from .conftest import create_simple_cubic_system_jax, requires_gpu
 
@@ -32,6 +36,24 @@
 dual_module = import_module("nvalchemiops.jax.neighbors.naive_dual_cutoff")
 
 
+def _active_neighbor_shift_rows(
+    neighbor_matrix: jax.Array,
+    shifts: jax.Array,
+    counts: jax.Array,
+    atom_index: int,
+) -> list[tuple[int, int, int, int]]:
+    """Return sorted active ``(neighbor, sx, sy, sz)`` rows for one atom."""
+    count = int(np.asarray(counts[atom_index]))
+    rows = np.concatenate(
+        (
+            np.asarray(neighbor_matrix[atom_index, :count])[:, None],
+            np.asarray(shifts[atom_index, :count]),
+        ),
+        axis=1,
+    )
+    return sorted(tuple(row) for row in rows.tolist())
+
+
 class TestNaiveDualCutoffCorrectness:
     """Test correctness of naive dual cutoff neighbor list."""
 
@@ -123,6 +145,92 @@ def test_matrix_format_with_pbc(self):
         assert jnp.all(num_neighbors2 >= 0)
         assert jnp.all(num_neighbors2 >= num_neighbors1)
 
+    @pytest.mark.parametrize("dtype", [jnp.float32, jnp.float64])
+    @pytest.mark.parametrize(
+        "wrap_kwargs",
+        [
+            pytest.param({}, id="default_wrapped"),
+            pytest.param({"wrap_positions": False}, id="prewrapped"),
+        ],
+    )
+    def test_pbc_wrap_modes_match_single_cutoff(
+        self,
+        dtype,
+        wrap_kwargs,
+    ):
+        """Both PBC wrap modes match independent single-cutoff results."""
+        positions, cell, pbc = create_simple_cubic_system_jax(
+            num_atoms=8,
+            cell_size=2.0,
+            dtype=dtype,
+        )
+        cutoff1 = 1.1
+        cutoff2 = 1.5
+        max_neighbors1 = 15
+        max_neighbors2 = 25
+
+        (
+            neighbor_matrix1,
+            num_neighbors1,
+            neighbor_matrix_shifts1,
+            neighbor_matrix2,
+            num_neighbors2,
+            neighbor_matrix_shifts2,
+        ) = naive_neighbor_list_dual_cutoff(
+            positions,
+            cutoff1=cutoff1,
+            cutoff2=cutoff2,
+            cell=cell,
+            pbc=pbc,
+            max_neighbors1=max_neighbors1,
+            max_neighbors2=max_neighbors2,
+            **wrap_kwargs,
+        )
+        reference1 = naive_neighbor_list(
+            positions,
+            cutoff=cutoff1,
+            cell=cell,
+            pbc=pbc,
+            max_neighbors=max_neighbors1,
+            **wrap_kwargs,
+        )
+        reference2 = naive_neighbor_list(
+            positions,
+            cutoff=cutoff2,
+            cell=cell,
+            pbc=pbc,
+            max_neighbors=max_neighbors2,
+            **wrap_kwargs,
+        )
+
+        assert jnp.any(num_neighbors1 > 0)
+        assert jnp.any(num_neighbors2 > 0)
+        assert jnp.array_equal(num_neighbors1, reference1[1])
+        assert jnp.array_equal(num_neighbors2, reference2[1])
+        for atom_index in range(positions.shape[0]):
+            assert _active_neighbor_shift_rows(
+                neighbor_matrix1,
+                neighbor_matrix_shifts1,
+                num_neighbors1,
+                atom_index,
+            ) == _active_neighbor_shift_rows(
+                reference1[0],
+                reference1[2],
+                reference1[1],
+                atom_index,
+            )
+            assert _active_neighbor_shift_rows(
+                neighbor_matrix2,
+                neighbor_matrix_shifts2,
+                num_neighbors2,
+                atom_index,
+            ) == _active_neighbor_shift_rows(
+                reference2[0],
+                reference2[2],
+                reference2[1],
+                atom_index,
+            )
+
 
 class TestNaiveDualCutoffEdgeCases:
     """Test edge cases for naive dual cutoff neighbor list."""
@@ -326,6 +434,84 @@ def jitted_dual(positions):
 class TestNaiveDualCutoffSelectiveRebuildFlags:
     """Test selective rebuild (rebuild_flags) for naive_neighbor_list_dual_cutoff JAX."""
 
+    def test_wrapped_pbc_false_flag_preserves_both_cutoffs(self, dtype):
+        """Wrapped PBC false flags preserve both matrices, counts, and shifts."""
+        positions, cell, pbc = create_simple_cubic_system_jax(
+            num_atoms=8, cell_size=2.0, dtype=dtype
+        )
+        positions = positions.at[0, 0].add(2.0)
+        kwargs = {
+            "cell": cell,
+            "pbc": pbc,
+            "max_neighbors1": 15,
+            "max_neighbors2": 25,
+        }
+        reference = naive_neighbor_list_dual_cutoff(positions, 1.1, 1.5, **kwargs)
+        preserved = naive_neighbor_list_dual_cutoff(
+            positions,
+            1.1,
+            1.5,
+            neighbor_matrix1=reference[0],
+            num_neighbors1=reference[1],
+            neighbor_matrix_shifts1=reference[2],
+            neighbor_matrix2=reference[3],
+            num_neighbors2=reference[4],
+            neighbor_matrix_shifts2=reference[5],
+            rebuild_flags=jnp.zeros(1, dtype=jnp.bool_),
+            **kwargs,
+        )
+
+        for got, expected in zip(preserved, reference):
+            assert jnp.array_equal(got, expected)
+
+    def test_wrapped_pbc_true_flag_rebuilds_both_cutoffs(self, dtype):
+        """Wrapped PBC true flags rebuild both cutoff matrices and shifts."""
+        positions, cell, pbc = create_simple_cubic_system_jax(
+            num_atoms=8, cell_size=2.0, dtype=dtype
+        )
+        positions = positions.at[0, 0].add(2.0)
+        kwargs = {
+            "cell": cell,
+            "pbc": pbc,
+            "max_neighbors1": 15,
+            "max_neighbors2": 25,
+        }
+        reference = naive_neighbor_list_dual_cutoff(positions, 1.1, 1.5, **kwargs)
+        rebuilt = naive_neighbor_list_dual_cutoff(
+            positions,
+            1.1,
+            1.5,
+            neighbor_matrix1=jnp.zeros_like(reference[0]),
+            num_neighbors1=jnp.zeros_like(reference[1]),
+            neighbor_matrix_shifts1=jnp.zeros_like(reference[2]),
+            neighbor_matrix2=jnp.zeros_like(reference[3]),
+            num_neighbors2=jnp.zeros_like(reference[4]),
+            neighbor_matrix_shifts2=jnp.zeros_like(reference[5]),
+            rebuild_flags=jnp.ones(1, dtype=jnp.bool_),
+            **kwargs,
+        )
+
+        for matrix_index, count_index, shifts_index in ((0, 1, 2), (3, 4, 5)):
+            assert jnp.array_equal(rebuilt[count_index], reference[count_index])
+            for atom_index, count in enumerate(reference[count_index].tolist()):
+                rebuilt_pairs = jnp.concatenate(
+                    [
+                        rebuilt[matrix_index][atom_index, :count, None],
+                        rebuilt[shifts_index][atom_index, :count],
+                    ],
+                    axis=1,
+                )
+                reference_pairs = jnp.concatenate(
+                    [
+                        reference[matrix_index][atom_index, :count, None],
+                        reference[shifts_index][atom_index, :count],
+                    ],
+                    axis=1,
+                )
+                assert {tuple(row) for row in rebuilt_pairs.tolist()} == {
+                    tuple(row) for row in reference_pairs.tolist()
+                }
+
     def test_no_rebuild_preserves_data(self, dtype):
         """Flag=False: neighbor data should remain unchanged."""
         positions, _, _ = create_simple_cubic_system_jax(
diff --git a/test/neighbors/bindings/torch/test_batch_cluster_tile.py b/test/neighbors/bindings/torch/test_batch_cluster_tile.py
index d012d9cd..d6b715b4 100644
--- a/test/neighbors/bindings/torch/test_batch_cluster_tile.py
+++ b/test/neighbors/bindings/torch/test_batch_cluster_tile.py
@@ -33,8 +33,10 @@
     batch_cluster_tile_neighbor_list,
     batch_query_cluster_tile,
     estimate_batch_cluster_tile_list_sizes,
+    estimate_batch_cluster_tile_segments,
     estimate_batch_max_tiles_per_group,
 )
+from nvalchemiops.torch.neighbors.neighbor_utils import _validate_segmented_coo_state
 
 from ...test_utils import (
     assert_neighbor_lists_equal,
@@ -79,9 +81,70 @@ def _make_batch(
     return positions, cell_batch, batch_ptr
 
 
+def _scratch_kwargs(
+    scratch: tuple[torch.Tensor, ...],
+) -> dict[str, torch.Tensor]:
+    """Map allocator outputs to one-shot API keyword arguments."""
+    names = (
+        "sorted_atom_index",
+        "sort_inv",
+        "sorted_pos_x",
+        "sorted_pos_y",
+        "sorted_pos_z",
+        "batch_idx_sorted",
+        "batch_ptr_padded",
+        "group_system",
+        "group_ptr",
+        "group_ctr_x",
+        "group_ctr_y",
+        "group_ctr_z",
+        "group_ext_x",
+        "group_ext_y",
+        "group_ext_z",
+        "num_tiles",
+        "tile_row_group",
+        "tile_col_group",
+        "tile_system",
+    )
+    return dict(zip(names, scratch))
+
+
 class TestBatchClusterTileValidation:
     """Validate public option combinations rejected before kernel launch."""
 
+    def test_selective_coo_bootstrap_rejects_partial_state(self):
+        """All-true COO bootstrap does not mix caller and allocated state."""
+        positions, cell_batch, batch_ptr = _make_batch([32], [8.0], device="cpu")
+
+        with pytest.raises(ValueError, match="complete caller-owned state"):
+            batch_cluster_tile_neighbor_list(
+                positions,
+                2.0,
+                cell_batch,
+                batch_ptr,
+                format="coo",
+                rebuild_flags=torch.ones(1, dtype=torch.bool),
+                return_state=True,
+                neighbor_list=torch.empty((2, 64), dtype=torch.int32),
+            )
+
+    def test_segmented_coo_validation_allows_nonselective_state(self):
+        """Segmented COO metadata does not require selective rebuild flags."""
+        device = torch.device("cpu")
+        capacity = 8
+
+        result = _validate_segmented_coo_state(
+            device=device,
+            num_systems=1,
+            neighbor_list=torch.empty((2, capacity), dtype=torch.int32),
+            neighbor_list_shifts=torch.empty((capacity, 3), dtype=torch.int32),
+            pair_offsets=torch.tensor([0, capacity], dtype=torch.int32),
+            pair_counts=torch.zeros(1, dtype=torch.int32),
+            rebuild_flags=None,
+        )
+
+        assert result == capacity
+
     def test_pair_outputs_reject_tile_format(self):
         """Pair-output buffers are not supported with tile-format output."""
         positions, cell_batch, batch_ptr = _make_batch([32], [8.0], device="cpu")
@@ -157,6 +220,146 @@ def test_segmented_offsets_are_all_or_nothing(self):
                 tile_counts=counts,
             )
 
+    def test_return_state_requires_rebuild_flags(self):
+        """State return is rejected without selective-rebuild inputs."""
+        positions, cell_batch, batch_ptr = _make_batch([32], [8.0], device="cpu")
+
+        with pytest.raises(ValueError, match="requires rebuild_flags"):
+            batch_cluster_tile_neighbor_list(
+                positions,
+                2.0,
+                cell_batch,
+                batch_ptr,
+                return_state=True,
+            )
+
+    def test_segmented_coo_rejects_undersized_topology_before_build(self, monkeypatch):
+        """Segmented COO capacity mismatches fail before a Warp build launch."""
+        positions, cell_batch, batch_ptr = _make_batch([32], [8.0], device="cpu")
+        pair_offsets = torch.tensor([0, 64], dtype=torch.int32)
+        tile_offsets = torch.tensor([0, 4], dtype=torch.int32)
+
+        def fail_build(*args, **kwargs):
+            del args, kwargs
+            raise AssertionError("segmented COO validation must run before build")
+
+        monkeypatch.setattr(
+            batch_cluster_tile_module,
+            "batch_build_cluster_tile_list",
+            fail_build,
+        )
+
+        with pytest.raises(ValueError, match="neighbor_list.*capacity"):
+            batch_cluster_tile_neighbor_list(
+                positions,
+                2.0,
+                cell_batch,
+                batch_ptr,
+                format="coo",
+                rebuild_flags=torch.ones(1, dtype=torch.bool),
+                neighbor_list=torch.empty((2, 1), dtype=torch.int32),
+                neighbor_list_shifts=torch.empty((64, 3), dtype=torch.int32),
+                pair_offsets=pair_offsets,
+                pair_counts=torch.zeros(1, dtype=torch.int32),
+                tile_offsets=tile_offsets,
+                tile_counts=torch.zeros(1, dtype=torch.int32),
+            )
+
+    def test_segmented_coo_false_flag_requires_prior_state(self, monkeypatch):
+        """A skipped system cannot preserve implicitly allocated COO state."""
+        positions, cell_batch, batch_ptr = _make_batch([32], [8.0], device="cpu")
+
+        def fail_build(*args, **kwargs):
+            del args, kwargs
+            raise AssertionError("selective state validation must run before build")
+
+        monkeypatch.setattr(
+            batch_cluster_tile_module,
+            "batch_build_cluster_tile_list",
+            fail_build,
+        )
+
+        with pytest.raises(ValueError, match="caller-owned"):
+            batch_cluster_tile_neighbor_list(
+                positions,
+                2.0,
+                cell_batch,
+                batch_ptr,
+                format="coo",
+                rebuild_flags=torch.zeros(1, dtype=torch.bool),
+                pair_offsets=torch.tensor([0, 64], dtype=torch.int32),
+                pair_counts=torch.zeros(1, dtype=torch.int32),
+                tile_offsets=torch.tensor([0, 4], dtype=torch.int32),
+                tile_counts=torch.zeros(1, dtype=torch.int32),
+            )
+
+    @pytest.mark.parametrize(
+        ("mutate", "message"),
+        [
+            (
+                lambda kwargs: kwargs.update(
+                    pair_offsets=torch.tensor([1, 64], dtype=torch.int32),
+                ),
+                "pair_offsets must start",
+            ),
+            (
+                lambda kwargs: kwargs.update(
+                    pair_counts=torch.tensor([65], dtype=torch.int32),
+                ),
+                "pair_counts must lie",
+            ),
+            (
+                lambda kwargs: kwargs.update(
+                    neighbor_list_shifts=torch.empty((63, 3), dtype=torch.int32),
+                ),
+                "neighbor_list_shifts",
+            ),
+            (
+                lambda kwargs: kwargs.update(
+                    rebuild_flags=torch.ones(1, dtype=torch.int32),
+                ),
+                "rebuild_flags",
+            ),
+        ],
+    )
+    def test_segmented_coo_rejects_malformed_state_before_build(
+        self,
+        monkeypatch,
+        mutate,
+        message,
+    ):
+        """Malformed segmented metadata is rejected before any Warp launch."""
+        positions, cell_batch, batch_ptr = _make_batch([32], [8.0], device="cpu")
+        kwargs = {
+            "format": "coo",
+            "rebuild_flags": torch.ones(1, dtype=torch.bool),
+            "neighbor_list": torch.empty((2, 64), dtype=torch.int32),
+            "neighbor_list_shifts": torch.empty((64, 3), dtype=torch.int32),
+            "pair_offsets": torch.tensor([0, 64], dtype=torch.int32),
+            "pair_counts": torch.zeros(1, dtype=torch.int32),
+            "tile_offsets": torch.tensor([0, 4], dtype=torch.int32),
+            "tile_counts": torch.zeros(1, dtype=torch.int32),
+        }
+        mutate(kwargs)
+
+        def fail_build(*args, **kwargs):
+            del args, kwargs
+            raise AssertionError("segmented COO validation must run before build")
+
+        monkeypatch.setattr(
+            batch_cluster_tile_module,
+            "batch_build_cluster_tile_list",
+            fail_build,
+        )
+        with pytest.raises(ValueError, match=message):
+            batch_cluster_tile_neighbor_list(
+                positions,
+                2.0,
+                cell_batch,
+                batch_ptr,
+                **kwargs,
+            )
+
 
 def test_batch_cluster_tile_max_tiles_per_group_bypasses_sizing(monkeypatch):
     """Explicit max_tiles_per_group skips geometry-aware sizing sync."""
@@ -771,6 +974,343 @@ def test_format_tile_returns_eleven_tuple(self, device, dtype):
         assert group_ptr.shape[0] == len(sizes) + 1
         assert n_padded >= N
 
+    def test_selective_matrix_return_state_appends_batched_state(self, device, dtype):
+        """Selective matrix calls can return reusable batched tile state."""
+        positions, cell_batch, batch_ptr = _make_batch(
+            [64], [8.0], device=device, dtype=dtype, seed=43
+        )
+
+        out = batch_cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell_batch,
+            batch_ptr,
+            max_neighbors=64,
+            rebuild_flags=torch.ones(1, dtype=torch.bool, device=device),
+            return_state=True,
+        )
+
+        assert len(out) == 9
+        neighbor_matrix, num_neighbors, shifts, *state = out
+        assert neighbor_matrix.shape == (64, 64)
+        assert num_neighbors.shape == (64,)
+        assert shifts.shape == (64, 64, 3)
+        tile_offsets, tile_counts, num_tiles, row, col, system = state
+        assert tile_offsets.shape == (2,)
+        assert tile_counts.shape == (1,)
+        assert num_tiles.shape == (1,)
+        assert row.ndim == col.ndim == system.ndim == 1
+
+    def test_selective_matrix_default_arity_is_unchanged(self, device, dtype):
+        """Selective matrix calls keep the three-array default return."""
+        positions, cell_batch, batch_ptr = _make_batch(
+            [64], [8.0], device=device, dtype=dtype, seed=44
+        )
+
+        out = batch_cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell_batch,
+            batch_ptr,
+            max_neighbors=64,
+            rebuild_flags=torch.ones(1, dtype=torch.bool, device=device),
+        )
+
+        assert len(out) == 3
+
+    @pytest.mark.parametrize("rebuild_flag", [False, True])
+    def test_returned_matrix_state_can_be_passed_to_next_call(
+        self, device, dtype, rebuild_flag
+    ):
+        """The public state suffix is sufficient for the next selective call."""
+        positions, cell_batch, batch_ptr = _make_batch(
+            [64], [8.0], device=device, dtype=dtype, seed=48
+        )
+        first = batch_cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell_batch,
+            batch_ptr,
+            max_neighbors=64,
+            rebuild_flags=torch.ones(1, dtype=torch.bool, device=device),
+            return_state=True,
+        )
+        matrix, counts, shifts, *state = first
+        tile_offsets, tile_counts, num_tiles, row, col, system = state
+
+        second = batch_cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell_batch,
+            batch_ptr,
+            max_neighbors=64,
+            rebuild_flags=torch.tensor([rebuild_flag], dtype=torch.bool, device=device),
+            neighbor_matrix=matrix,
+            num_neighbors=counts,
+            neighbor_matrix_shifts=shifts,
+            tile_offsets=tile_offsets,
+            tile_counts=tile_counts,
+            num_tiles=num_tiles,
+            tile_row_group=row,
+            tile_col_group=col,
+            tile_system=system,
+            return_state=True,
+        )
+
+        for returned, supplied in zip(second[-6:], state):
+            assert returned is supplied
+
+    def test_return_state_mixed_flags_preserve_unflagged_system(self, device, dtype):
+        """Mixed flags keep unflagged rows while returning reusable state."""
+        positions, cell_batch, batch_ptr = _make_batch(
+            [64, 64], [8.0, 8.0], device=device, dtype=dtype, seed=49
+        )
+        first = batch_cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell_batch,
+            batch_ptr,
+            max_neighbors=64,
+            rebuild_flags=torch.ones(2, dtype=torch.bool, device=device),
+            return_state=True,
+        )
+        matrix, counts, shifts, *state = first
+        system0_matrix = matrix[:64].clone()
+        system0_counts = counts[:64].clone()
+        system0_shifts = shifts[:64].clone()
+        system0_tile_count = state[1][0].clone()
+        moved = positions.clone()
+        moved[:64] = torch.roll(moved[:64], shifts=1, dims=0)
+        moved[64:, 0] += 0.2
+
+        second = batch_cluster_tile_neighbor_list(
+            moved,
+            2.0,
+            cell_batch,
+            batch_ptr,
+            max_neighbors=64,
+            rebuild_flags=torch.tensor([False, True], dtype=torch.bool, device=device),
+            neighbor_matrix=matrix,
+            num_neighbors=counts,
+            neighbor_matrix_shifts=shifts,
+            tile_offsets=state[0],
+            tile_counts=state[1],
+            num_tiles=state[2],
+            tile_row_group=state[3],
+            tile_col_group=state[4],
+            tile_system=state[5],
+            return_state=True,
+        )
+
+        assert len(second) == 9
+        torch.testing.assert_close(second[0][:64], system0_matrix)
+        torch.testing.assert_close(second[1][:64], system0_counts)
+        torch.testing.assert_close(second[2][:64], system0_shifts)
+        torch.testing.assert_close(second[4][0], system0_tile_count)
+        for returned, supplied in zip(second[-6:], state):
+            assert returned is supplied
+
+    @pytest.mark.parametrize("rebuild_flag", [False, True])
+    def test_selective_matrix_state_aliases_caller_buffers(
+        self, device, dtype, rebuild_flag
+    ):
+        """Returned matrix state aliases explicitly supplied scratch buffers."""
+        positions, cell_batch, batch_ptr = _make_batch(
+            [64], [8.0], device=device, dtype=dtype, seed=45
+        )
+        scratch_kwargs = _scratch_kwargs(
+            allocate_batch_cluster_tile_list(
+                batch_ptr, torch.device(device), dtype=dtype
+            )
+        )
+        _tile_caps, tile_offsets, _pair_caps, _pair_offsets = (
+            estimate_batch_cluster_tile_segments(batch_ptr, max_neighbors=64)
+        )
+        tile_counts = torch.zeros(1, dtype=torch.int32, device=device)
+        matrix = torch.full((64, 64), 64, dtype=torch.int32, device=device)
+        counts = torch.zeros(64, dtype=torch.int32, device=device)
+        shifts = torch.zeros((64, 64, 3), dtype=torch.int32, device=device)
+
+        out = batch_cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell_batch,
+            batch_ptr,
+            max_neighbors=64,
+            rebuild_flags=torch.tensor([rebuild_flag], dtype=torch.bool, device=device),
+            neighbor_matrix=matrix,
+            num_neighbors=counts,
+            neighbor_matrix_shifts=shifts,
+            tile_offsets=tile_offsets,
+            tile_counts=tile_counts,
+            return_state=True,
+            **scratch_kwargs,
+        )
+
+        assert len(out) == 9
+        for returned, supplied in zip(
+            out[-6:],
+            (
+                tile_offsets,
+                tile_counts,
+                scratch_kwargs["num_tiles"],
+                scratch_kwargs["tile_row_group"],
+                scratch_kwargs["tile_col_group"],
+                scratch_kwargs["tile_system"],
+            ),
+        ):
+            assert returned is supplied
+
+    @pytest.mark.parametrize("rebuild_flag", [False, True])
+    def test_selective_segmented_coo_return_state(self, device, dtype, rebuild_flag):
+        """Segmented COO appends six caller-owned batched state buffers."""
+        positions, cell_batch, batch_ptr = _make_batch(
+            [64], [8.0], device=device, dtype=dtype, seed=46
+        )
+        scratch_kwargs = _scratch_kwargs(
+            allocate_batch_cluster_tile_list(
+                batch_ptr, torch.device(device), dtype=dtype
+            )
+        )
+        _tile_caps, tile_offsets, _pair_caps, pair_offsets = (
+            estimate_batch_cluster_tile_segments(batch_ptr, max_neighbors=64)
+        )
+        max_pairs = int(pair_offsets[-1].item())
+        tile_counts = torch.zeros(1, dtype=torch.int32, device=device)
+        pair_counts = torch.zeros(1, dtype=torch.int32, device=device)
+
+        out = batch_cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell_batch,
+            batch_ptr,
+            max_neighbors=64,
+            format="coo",
+            rebuild_flags=torch.tensor([rebuild_flag], dtype=torch.bool, device=device),
+            neighbor_list=torch.empty((2, max_pairs), dtype=torch.int32, device=device),
+            neighbor_list_shifts=torch.empty(
+                (max_pairs, 3), dtype=torch.int32, device=device
+            ),
+            pair_counter=torch.zeros(1, dtype=torch.int32, device=device),
+            tile_offsets=tile_offsets,
+            tile_counts=tile_counts,
+            pair_offsets=pair_offsets,
+            pair_counts=pair_counts,
+            return_state=True,
+            **scratch_kwargs,
+        )
+
+        assert len(out) == 10
+        assert out[1] is pair_offsets
+        assert out[2] is pair_counts
+        for returned, supplied in zip(
+            out[-6:],
+            (
+                tile_offsets,
+                tile_counts,
+                scratch_kwargs["num_tiles"],
+                scratch_kwargs["tile_row_group"],
+                scratch_kwargs["tile_col_group"],
+                scratch_kwargs["tile_system"],
+            ),
+        ):
+            assert returned is supplied
+
+    def test_selective_segmented_coo_all_true_bootstrap_allocates_state(
+        self, device, dtype
+    ):
+        """An all-true batched COO call allocates persistent state."""
+        positions, cell_batch, batch_ptr = _make_batch(
+            [64], [8.0], device=device, dtype=dtype, seed=46
+        )
+
+        result = batch_cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell_batch,
+            batch_ptr,
+            max_neighbors=64,
+            format="coo",
+            rebuild_flags=torch.ones(1, dtype=torch.bool, device=device),
+            return_state=True,
+        )
+
+        assert len(result) == 10
+        neighbor_list, pair_offsets, pair_counts, shifts, *state = result
+        assert neighbor_list.shape[0] == 2
+        assert shifts.shape == (neighbor_list.shape[1], 3)
+        assert pair_offsets.shape == (2,)
+        assert pair_counts.shape == (1,)
+
+        reused = batch_cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell_batch,
+            batch_ptr,
+            max_neighbors=64,
+            format="coo",
+            rebuild_flags=torch.zeros(1, dtype=torch.bool, device=device),
+            return_state=True,
+            neighbor_list=neighbor_list,
+            neighbor_list_shifts=shifts,
+            pair_offsets=pair_offsets,
+            pair_counts=pair_counts,
+            tile_offsets=state[0],
+            tile_counts=state[1],
+            num_tiles=state[2],
+            tile_row_group=state[3],
+            tile_col_group=state[4],
+            tile_system=state[5],
+        )
+
+        assert all(reused[index] is result[index] for index in range(10))
+
+    def test_selective_dual_cutoff_return_state(self, device, dtype):
+        """Dual-cutoff matrix output appends state after both matrix triples."""
+        positions, cell_batch, batch_ptr = _make_batch(
+            [64], [8.0], device=device, dtype=dtype, seed=47
+        )
+        scratch_kwargs = _scratch_kwargs(
+            allocate_batch_cluster_tile_list(
+                batch_ptr, torch.device(device), dtype=dtype
+            )
+        )
+        _tile_caps, tile_offsets, _pair_caps, _pair_offsets = (
+            estimate_batch_cluster_tile_segments(batch_ptr, max_neighbors=64)
+        )
+        tile_counts = torch.zeros(1, dtype=torch.int32, device=device)
+        matrix1 = torch.full((64, 64), 64, dtype=torch.int32, device=device)
+        counts1 = torch.zeros(64, dtype=torch.int32, device=device)
+        shifts1 = torch.zeros((64, 64, 3), dtype=torch.int32, device=device)
+        matrix2 = matrix1.clone()
+        counts2 = counts1.clone()
+        shifts2 = shifts1.clone()
+
+        out = batch_cluster_tile_neighbor_list(
+            positions,
+            1.5,
+            cell_batch,
+            batch_ptr,
+            cutoff2=2.0,
+            max_neighbors=64,
+            rebuild_flags=torch.ones(1, dtype=torch.bool, device=device),
+            neighbor_matrix=matrix1,
+            num_neighbors=counts1,
+            neighbor_matrix_shifts=shifts1,
+            neighbor_matrix2=matrix2,
+            num_neighbors2=counts2,
+            neighbor_matrix_shifts2=shifts2,
+            tile_offsets=tile_offsets,
+            tile_counts=tile_counts,
+            return_state=True,
+            **scratch_kwargs,
+        )
+
+        assert len(out) == 12
+        assert out[6] is tile_offsets
+        assert out[7] is tile_counts
+        assert out[8] is scratch_kwargs["num_tiles"]
+
     def test_format_coo_returns_three_tuple(self, device, dtype):
         sizes = [64, 64]
         positions, cell_batch, batch_ptr = _make_batch(
diff --git a/test/neighbors/bindings/torch/test_cluster_tile.py b/test/neighbors/bindings/torch/test_cluster_tile.py
index b46f7ac3..b141862a 100644
--- a/test/neighbors/bindings/torch/test_cluster_tile.py
+++ b/test/neighbors/bindings/torch/test_cluster_tile.py
@@ -27,6 +27,7 @@
     cluster_tile_neighbor_list,
     estimate_cluster_tile_list_sizes,
     query_cluster_tile,
+    query_cluster_tile_coo,
 )
 from nvalchemiops.torch.neighbors.naive import naive_neighbor_list
 
@@ -586,6 +587,364 @@ def compiled_build_and_convert(positions, cutoff, cell, nm, nn, nms):
             s_got = {int(x.item()) for x in nm[i, :n_i]}
             assert s_ref == s_got, f"Row {i} neighbor set mismatch under torch.compile"
 
+    @pytest.mark.slow
+    def test_segmented_coo_query_compile(self, device, dtype):
+        """The selective segmented COO custom op matches eager execution."""
+        torch.manual_seed(5)
+        natom = 64
+        max_pairs = natom * 64
+        positions = torch.rand(natom, 3, dtype=dtype, device=device) * 6.0
+        cell = _orthorhombic_cell(6.0, device, dtype)
+        state = allocate_cluster_tile_list(natom, torch.device(device), dtype=dtype)
+        build_cluster_tile_list(positions, 2.0, cell, *state)
+        rebuild_flags = torch.tensor([True], dtype=torch.bool, device=device)
+        pair_offsets = torch.tensor([0, max_pairs], dtype=torch.int32, device=device)
+
+        eager_counter = torch.zeros(1, dtype=torch.int32, device=device)
+        eager_counts = torch.zeros(1, dtype=torch.int32, device=device)
+        eager_list = torch.empty((max_pairs, 2), dtype=torch.int32, device=device)
+        eager_shifts = torch.empty((max_pairs, 3), dtype=torch.int32, device=device)
+        query_cluster_tile_coo(
+            state[0],
+            state[2],
+            state[3],
+            state[4],
+            state[11],
+            state[12],
+            state[13],
+            cell,
+            2.0,
+            natom,
+            max_pairs,
+            eager_counter,
+            eager_list,
+            eager_shifts,
+            rebuild_flags=rebuild_flags,
+            pair_offsets=pair_offsets,
+            pair_counts=eager_counts,
+        )
+
+        @torch.compile
+        def compiled_query(pair_counter, pair_counts, coo_list, coo_shifts):
+            query_cluster_tile_coo(
+                state[0],
+                state[2],
+                state[3],
+                state[4],
+                state[11],
+                state[12],
+                state[13],
+                cell,
+                2.0,
+                natom,
+                max_pairs,
+                pair_counter,
+                coo_list,
+                coo_shifts,
+                rebuild_flags=rebuild_flags,
+                pair_offsets=pair_offsets,
+                pair_counts=pair_counts,
+            )
+            return pair_counts, coo_list, coo_shifts
+
+        compiled_counter = torch.zeros(1, dtype=torch.int32, device=device)
+        compiled_counts = torch.zeros(1, dtype=torch.int32, device=device)
+        compiled_list = torch.empty((max_pairs, 2), dtype=torch.int32, device=device)
+        compiled_shifts = torch.empty((max_pairs, 3), dtype=torch.int32, device=device)
+        compiled_counts, compiled_list, compiled_shifts = compiled_query(
+            compiled_counter,
+            compiled_counts,
+            compiled_list,
+            compiled_shifts,
+        )
+
+        pair_count = int(eager_counts.item())
+        assert torch.equal(eager_counts, compiled_counts)
+        eager_pairs = {
+            tuple(entry)
+            for entry in torch.cat(
+                [eager_list[:pair_count], eager_shifts[:pair_count]], dim=1
+            )
+            .cpu()
+            .tolist()
+        }
+        compiled_pairs = {
+            tuple(entry)
+            for entry in torch.cat(
+                [compiled_list[:pair_count], compiled_shifts[:pair_count]], dim=1
+            )
+            .cpu()
+            .tolist()
+        }
+        assert eager_pairs == compiled_pairs
+
+    @pytest.mark.slow
+    def test_selective_coo_wrapper_fullgraph_compile(self, device, dtype):
+        """The public selective COO wrapper supports true and false compiled flags."""
+        torch.manual_seed(23)
+        natom = 64
+        max_pairs = natom * 64
+        positions = torch.rand(natom, 3, dtype=dtype, device=device) * 6.0
+        cell = _orthorhombic_cell(6.0, device, dtype)
+        num_tiles, tile_row_group, tile_col_group, *_ = cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell,
+            format="tile",
+        )
+        pair_offsets = torch.tensor([0, max_pairs], dtype=torch.int32, device=device)
+
+        eager_list = torch.empty((2, max_pairs), dtype=torch.int32, device=device)
+        eager_counts = torch.zeros(1, dtype=torch.int32, device=device)
+        eager_shifts = torch.empty((max_pairs, 3), dtype=torch.int32, device=device)
+        cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell,
+            max_neighbors=64,
+            format="coo",
+            rebuild_flags=torch.ones(1, dtype=torch.bool, device=device),
+            return_state=True,
+            num_tiles=num_tiles.clone(),
+            tile_row_group=tile_row_group.clone(),
+            tile_col_group=tile_col_group.clone(),
+            neighbor_list=eager_list,
+            pair_offsets=pair_offsets,
+            pair_counts=eager_counts,
+            neighbor_list_shifts=eager_shifts,
+        )
+
+        compiled_list = torch.empty((2, max_pairs), dtype=torch.int32, device=device)
+        compiled_counts = torch.zeros(1, dtype=torch.int32, device=device)
+        compiled_shifts = torch.empty(
+            (max_pairs, 3),
+            dtype=torch.int32,
+            device=device,
+        )
+        compiled_num_tiles = num_tiles.clone()
+        compiled_row = tile_row_group.clone()
+        compiled_col = tile_col_group.clone()
+
+        @torch.compile(fullgraph=True)
+        def run(rebuild_flags):
+            return cluster_tile_neighbor_list(
+                positions,
+                2.0,
+                cell,
+                max_neighbors=64,
+                format="coo",
+                rebuild_flags=rebuild_flags,
+                return_state=True,
+                num_tiles=compiled_num_tiles,
+                tile_row_group=compiled_row,
+                tile_col_group=compiled_col,
+                neighbor_list=compiled_list,
+                pair_offsets=pair_offsets,
+                pair_counts=compiled_counts,
+                neighbor_list_shifts=compiled_shifts,
+            )
+
+        result = run(torch.ones(1, dtype=torch.bool, device=device))
+        assert result[0].data_ptr() == compiled_list.data_ptr()
+        assert result[2].data_ptr() == compiled_counts.data_ptr()
+        assert torch.equal(compiled_counts, eager_counts)
+        count = int(compiled_counts.item())
+        compiled_pairs = {
+            tuple(pair)
+            for pair in torch.cat(
+                [compiled_list[:, :count].T, compiled_shifts[:count]],
+                dim=1,
+            )
+            .cpu()
+            .tolist()
+        }
+        eager_pairs = {
+            tuple(pair)
+            for pair in torch.cat(
+                [eager_list[:, :count].T, eager_shifts[:count]],
+                dim=1,
+            )
+            .cpu()
+            .tolist()
+        }
+        assert compiled_pairs == eager_pairs
+
+        preserved_list = compiled_list.clone()
+        preserved_counts = compiled_counts.clone()
+        preserved_shifts = compiled_shifts.clone()
+        run(torch.zeros(1, dtype=torch.bool, device=device))
+        assert torch.equal(compiled_list, preserved_list)
+        assert torch.equal(compiled_counts, preserved_counts)
+        assert torch.equal(compiled_shifts, preserved_shifts)
+
+        capture_flags = torch.ones(1, dtype=torch.bool, device=device)
+        warmup_stream = torch.cuda.Stream(device=device)
+        with torch.cuda.stream(warmup_stream):
+            run(capture_flags)
+        warmup_stream.synchronize()
+        graph = torch.cuda.CUDAGraph()
+        torch.cuda.synchronize(device=device)
+        with torch.cuda.graph(graph):
+            run(capture_flags)
+        captured_list = compiled_list.clone()
+        captured_counts = compiled_counts.clone()
+        captured_shifts = compiled_shifts.clone()
+        capture_flags.zero_()
+        graph.replay()
+        torch.cuda.synchronize(device=device)
+        assert torch.equal(compiled_list, captured_list)
+        assert torch.equal(compiled_counts, captured_counts)
+        assert torch.equal(compiled_shifts, captured_shifts)
+
+    @pytest.mark.slow
+    @pytest.mark.parametrize(
+        ("pair_offsets_values", "rebuild_flag"),
+        [
+            ((64, 0), True),
+            ((0, 32), True),
+            ((1, 64), True),
+            ((64, 0), False),
+            ((0, 32), False),
+            ((1, 64), False),
+        ],
+        ids=[
+            "reversed-true",
+            "undersized-true",
+            "nonzero-start-true",
+            "reversed-false",
+            "undersized-false",
+            "nonzero-start-false",
+        ],
+    )
+    def test_selective_coo_wrapper_fullgraph_invalid_offsets_fail_closed(
+        self,
+        device,
+        dtype,
+        pair_offsets_values,
+        rebuild_flag,
+    ):
+        """Compiled malformed offsets must not name active COO entries."""
+        natom = 64
+        capacity = 64
+        positions = torch.zeros((natom, 3), dtype=dtype, device=device)
+        cell = _orthorhombic_cell(6.0, device, dtype)
+        num_tiles, tile_row_group, tile_col_group, *_ = cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell,
+            format="tile",
+        )
+        neighbor_list = torch.full(
+            (2, capacity),
+            -77,
+            dtype=torch.int32,
+            device=device,
+        )
+        neighbor_list_shifts = torch.full(
+            (capacity, 3),
+            -77,
+            dtype=torch.int32,
+            device=device,
+        )
+        pair_offsets = torch.tensor(
+            pair_offsets_values,
+            dtype=torch.int32,
+            device=device,
+        )
+        pair_counts = torch.tensor(
+            [7 if not rebuild_flag else 0],
+            dtype=torch.int32,
+            device=device,
+        )
+
+        @torch.compile(fullgraph=True)
+        def run(
+            rebuild_flags,
+            runtime_pair_offsets,
+            runtime_pair_counts,
+        ):
+            return cluster_tile_neighbor_list(
+                positions,
+                2.0,
+                cell,
+                max_neighbors=64,
+                format="coo",
+                rebuild_flags=rebuild_flags,
+                return_state=True,
+                num_tiles=num_tiles,
+                tile_row_group=tile_row_group,
+                tile_col_group=tile_col_group,
+                neighbor_list=neighbor_list,
+                pair_offsets=runtime_pair_offsets,
+                pair_counts=runtime_pair_counts,
+                neighbor_list_shifts=neighbor_list_shifts,
+            )
+
+        result = run(
+            torch.tensor([rebuild_flag], dtype=torch.bool, device=device),
+            pair_offsets,
+            pair_counts,
+        )
+
+        assert result[2].data_ptr() == pair_counts.data_ptr()
+        assert int(pair_counts.item()) == 0
+        assert torch.all(neighbor_list == -77)
+        assert torch.all(neighbor_list_shifts == -77)
+
+    @pytest.mark.slow
+    def test_selective_coo_wrapper_fullgraph_clamps_overflow_count(self, device, dtype):
+        """Compiled valid segments report at most their written capacity."""
+        natom = 64
+        capacity = 64
+        positions = torch.zeros((natom, 3), dtype=dtype, device=device)
+        cell = _orthorhombic_cell(6.0, device, dtype)
+        num_tiles, tile_row_group, tile_col_group, *_ = cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell,
+            format="tile",
+        )
+        neighbor_list = torch.full(
+            (2, capacity),
+            -77,
+            dtype=torch.int32,
+            device=device,
+        )
+        neighbor_list_shifts = torch.full(
+            (capacity, 3),
+            -77,
+            dtype=torch.int32,
+            device=device,
+        )
+        pair_offsets = torch.tensor([0, capacity], dtype=torch.int32, device=device)
+        pair_counts = torch.zeros(1, dtype=torch.int32, device=device)
+
+        @torch.compile(fullgraph=True)
+        def run(runtime_pair_counts):
+            return cluster_tile_neighbor_list(
+                positions,
+                2.0,
+                cell,
+                max_neighbors=64,
+                format="coo",
+                rebuild_flags=torch.ones(1, dtype=torch.bool, device=device),
+                return_state=True,
+                num_tiles=num_tiles,
+                tile_row_group=tile_row_group,
+                tile_col_group=tile_col_group,
+                neighbor_list=neighbor_list,
+                pair_offsets=pair_offsets,
+                pair_counts=runtime_pair_counts,
+                neighbor_list_shifts=neighbor_list_shifts,
+            )
+
+        result = run(pair_counts)
+
+        assert result[2].data_ptr() == pair_counts.data_ptr()
+        assert int(pair_counts.item()) == capacity
+        assert torch.all(neighbor_list != -77)
+        assert torch.all(neighbor_list_shifts != -77)
+
 
 # =============================================================================
 # Left-handed cells
@@ -784,7 +1143,11 @@ def test_cutoff2_returns_two_matrix_groups(self, device, dtype):
         assert int(nn2.sum().item()) > 0
         assert nm1.shape == nm2.shape == (3, 8)
 
-    def test_rebuild_flags_false_preserves_previous_outputs(self, device, dtype):
+    @pytest.mark.parametrize("rebuild_flag", [False, True])
+    def test_return_state_preserves_caller_owned_buffers(
+        self, device, dtype, rebuild_flag
+    ):
+        """Selective matrix calls return the caller-owned outputs and tile state."""
         torch.manual_seed(101)
         positions = torch.rand(64, 3, dtype=dtype, device=device) * 6.0
         cell = _orthorhombic_cell(6.0, device, dtype)
@@ -797,22 +1160,421 @@ def test_rebuild_flags_false_preserves_previous_outputs(self, device, dtype):
 
         moved = positions.clone()
         moved[0, 0] += 0.25
-        nm2, nn2, shifts2 = cluster_tile_neighbor_list(
-            moved,
-            cutoff,
+        nm2, nn2, shifts2, returned_num_tiles, returned_row, returned_col = (
+            cluster_tile_neighbor_list(
+                moved,
+                cutoff,
+                cell,
+                max_neighbors=64,
+                rebuild_flags=torch.tensor(
+                    [rebuild_flag], dtype=torch.bool, device=device
+                ),
+                num_tiles=num_tiles,
+                tile_row_group=tile_row_group,
+                tile_col_group=tile_col_group,
+                neighbor_matrix=nm,
+                num_neighbors=nn,
+                neighbor_matrix_shifts=shifts,
+                return_state=True,
+            )
+        )
+        assert nm2 is nm
+        assert nn2 is nn
+        assert shifts2 is shifts
+        assert returned_num_tiles is num_tiles
+        assert returned_row is tile_row_group
+        assert returned_col is tile_col_group
+
+    def test_return_state_false_keeps_matrix_arity(self, device, dtype):
+        """Selective matrix calls keep the three-array default return."""
+        positions = torch.rand(64, 3, dtype=dtype, device=device) * 6.0
+        cell = _orthorhombic_cell(6.0, device, dtype)
+        outputs = cluster_tile_neighbor_list(positions, 2.0, cell, max_neighbors=64)
+        num_tiles, row, col, *_ = cluster_tile_neighbor_list(
+            positions, 2.0, cell, format="tile"
+        )
+
+        result = cluster_tile_neighbor_list(
+            positions,
+            2.0,
             cell,
             max_neighbors=64,
             rebuild_flags=torch.tensor([False], dtype=torch.bool, device=device),
             num_tiles=num_tiles,
+            tile_row_group=row,
+            tile_col_group=col,
+            neighbor_matrix=outputs[0],
+            num_neighbors=outputs[1],
+            neighbor_matrix_shifts=outputs[2],
+        )
+
+        assert len(result) == 3
+
+    def test_dual_cutoff_return_state_appends_tile_state(self, device, dtype):
+        """Dual-cutoff selective matrix calls append state after both triples."""
+        positions = torch.rand(64, 3, dtype=dtype, device=device) * 6.0
+        cell = _orthorhombic_cell(6.0, device, dtype)
+        outputs = cluster_tile_neighbor_list(
+            positions, 1.5, cell, max_neighbors=64, cutoff2=2.0
+        )
+        num_tiles, row, col, *_ = cluster_tile_neighbor_list(
+            positions, 2.0, cell, format="tile"
+        )
+
+        result = cluster_tile_neighbor_list(
+            positions,
+            1.5,
+            cell,
+            max_neighbors=64,
+            cutoff2=2.0,
+            rebuild_flags=torch.tensor([False], dtype=torch.bool, device=device),
+            num_tiles=num_tiles,
+            tile_row_group=row,
+            tile_col_group=col,
+            neighbor_matrix=outputs[0],
+            num_neighbors=outputs[1],
+            neighbor_matrix_shifts=outputs[2],
+            neighbor_matrix2=outputs[3],
+            num_neighbors2=outputs[4],
+            neighbor_matrix_shifts2=outputs[5],
+            return_state=True,
+        )
+
+        assert len(result) == 9
+        assert result[-3] is num_tiles
+        assert result[-2] is row
+        assert result[-1] is col
+
+    def test_return_state_requires_rebuild_flags(self, device, dtype):
+        """State return is rejected without selective-rebuild inputs."""
+        positions = torch.rand(64, 3, dtype=dtype, device=device) * 6.0
+        cell = _orthorhombic_cell(6.0, device, dtype)
+
+        with pytest.raises(ValueError, match="requires rebuild_flags"):
+            cluster_tile_neighbor_list(
+                positions, 2.0, cell, max_neighbors=64, return_state=True
+            )
+
+    def test_selective_coo_false_preserves_buffers(self, device, dtype):
+        """A false selective COO flag returns fixed caller-owned buffers."""
+        positions = torch.rand(64, 3, dtype=dtype, device=device) * 6.0
+        cell = _orthorhombic_cell(6.0, device, dtype)
+        max_pairs = 256
+        num_tiles, tile_row_group, tile_col_group, *_ = cluster_tile_neighbor_list(
+            positions, 2.0, cell, format="tile"
+        )
+        neighbor_list = torch.full((2, max_pairs), 41, dtype=torch.int32, device=device)
+        pair_offsets = torch.tensor([0, max_pairs], dtype=torch.int32, device=device)
+        pair_counts = torch.full((1,), 17, dtype=torch.int32, device=device)
+        neighbor_list_shifts = torch.full(
+            (max_pairs, 3), -9, dtype=torch.int32, device=device
+        )
+
+        result = cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell,
+            max_neighbors=64,
+            max_pairs=max_pairs,
+            format="coo",
+            rebuild_flags=torch.tensor([False], dtype=torch.bool, device=device),
+            num_tiles=num_tiles,
             tile_row_group=tile_row_group,
             tile_col_group=tile_col_group,
-            neighbor_matrix=nm,
-            num_neighbors=nn,
-            neighbor_matrix_shifts=shifts,
+            neighbor_list=neighbor_list,
+            pair_offsets=pair_offsets,
+            pair_counts=pair_counts,
+            neighbor_list_shifts=neighbor_list_shifts,
+        )
+
+        assert len(result) == 4
+        assert result[0] is neighbor_list
+        assert result[1] is pair_offsets
+        assert result[2] is pair_counts
+        assert result[3] is neighbor_list_shifts
+        assert int(pair_counts.item()) == 17
+        assert torch.all(neighbor_list == 41)
+        assert torch.all(neighbor_list_shifts == -9)
+
+    def test_query_cluster_tile_coo_selective_false_preserves_buffers(
+        self, device, dtype
+    ):
+        """A false segmented flag preserves caller-owned COO topology buffers."""
+        torch.manual_seed(31)
+        natom = 64
+        max_pairs = 128
+        positions = torch.rand(natom, 3, dtype=dtype, device=device) * 6.0
+        cell = _orthorhombic_cell(6.0, device, dtype)
+        state = allocate_cluster_tile_list(natom, torch.device(device), dtype=dtype)
+        build_cluster_tile_list(positions, 2.0, cell, *state)
+
+        pair_counter = torch.zeros(1, dtype=torch.int32, device=device)
+        pair_offsets = torch.tensor([0, max_pairs], dtype=torch.int32, device=device)
+        pair_counts = torch.full((1,), 17, dtype=torch.int32, device=device)
+        coo_list = torch.full((max_pairs, 2), 41, dtype=torch.int32, device=device)
+        coo_shifts = torch.full((max_pairs, 3), -9, dtype=torch.int32, device=device)
+
+        query_cluster_tile_coo(
+            state[0],
+            state[2],
+            state[3],
+            state[4],
+            state[11],
+            state[12],
+            state[13],
+            cell,
+            2.0,
+            natom,
+            max_pairs,
+            pair_counter,
+            coo_list,
+            coo_shifts,
+            rebuild_flags=torch.tensor([False], dtype=torch.bool, device=device),
+            pair_offsets=pair_offsets,
+            pair_counts=pair_counts,
         )
-        assert torch.equal(nm2, nm)
-        assert torch.equal(nn2, nn)
-        assert torch.equal(shifts2, shifts)
+
+        assert int(pair_counts.item()) == 17
+        assert torch.all(coo_list == 41)
+        assert torch.all(coo_shifts == -9)
+
+    @pytest.mark.parametrize(
+        "missing",
+        [
+            "neighbor_list",
+            "pair_offsets",
+            "pair_counts",
+            "neighbor_list_shifts",
+        ],
+    )
+    def test_selective_coo_requires_fixed_buffers(self, device, dtype, missing):
+        """Selective COO requires every fixed-capacity topology buffer."""
+        positions = torch.rand(64, 3, dtype=dtype, device=device) * 6.0
+        cell = _orthorhombic_cell(6.0, device, dtype)
+        max_pairs = 256
+        num_tiles, tile_row_group, tile_col_group, *_ = cluster_tile_neighbor_list(
+            positions, 2.0, cell, format="tile"
+        )
+        kwargs = {
+            "neighbor_list": torch.empty(
+                (2, max_pairs), dtype=torch.int32, device=device
+            ),
+            "pair_offsets": torch.tensor(
+                [0, max_pairs], dtype=torch.int32, device=device
+            ),
+            "pair_counts": torch.zeros(1, dtype=torch.int32, device=device),
+            "neighbor_list_shifts": torch.empty(
+                (max_pairs, 3), dtype=torch.int32, device=device
+            ),
+        }
+        del kwargs[missing]
+
+        with pytest.raises(ValueError, match=missing):
+            cluster_tile_neighbor_list(
+                positions,
+                2.0,
+                cell,
+                max_neighbors=64,
+                max_pairs=max_pairs,
+                format="coo",
+                rebuild_flags=torch.tensor([False], dtype=torch.bool, device=device),
+                num_tiles=num_tiles,
+                tile_row_group=tile_row_group,
+                tile_col_group=tile_col_group,
+                **kwargs,
+            )
+
+    def test_selective_coo_rejects_invalid_neighbor_list_shape(self, device, dtype):
+        """Selective COO reports malformed topology buffers as a ValueError."""
+        positions = torch.rand(64, 3, dtype=dtype, device=device) * 6.0
+        cell = _orthorhombic_cell(6.0, device, dtype)
+        max_pairs = 256
+        num_tiles, tile_row_group, tile_col_group, *_ = cluster_tile_neighbor_list(
+            positions, 2.0, cell, format="tile"
+        )
+
+        with pytest.raises(ValueError, match="neighbor_list"):
+            cluster_tile_neighbor_list(
+                positions,
+                2.0,
+                cell,
+                max_neighbors=64,
+                max_pairs=max_pairs,
+                format="coo",
+                rebuild_flags=torch.tensor([False], dtype=torch.bool, device=device),
+                num_tiles=num_tiles,
+                tile_row_group=tile_row_group,
+                tile_col_group=tile_col_group,
+                neighbor_list=torch.empty(max_pairs, dtype=torch.int32, device=device),
+                pair_offsets=torch.tensor(
+                    [0, max_pairs], dtype=torch.int32, device=device
+                ),
+                pair_counts=torch.zeros(1, dtype=torch.int32, device=device),
+                neighbor_list_shifts=torch.empty(
+                    (max_pairs, 3), dtype=torch.int32, device=device
+                ),
+            )
+
+    def test_selective_coo_return_state_lifecycle(self, device, dtype):
+        """Selective COO matches compact pairs and reuses returned state."""
+        torch.manual_seed(33)
+        natom = 64
+        max_pairs = natom * 64
+        positions = torch.rand(natom, 3, dtype=dtype, device=device) * 6.0
+        cell = _orthorhombic_cell(6.0, device, dtype)
+        num_tiles, tile_row_group, tile_col_group, *_ = cluster_tile_neighbor_list(
+            positions, 2.0, cell, format="tile"
+        )
+        compact_list, _compact_ptr, compact_shifts = cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell,
+            max_neighbors=64,
+            max_pairs=max_pairs,
+            format="coo",
+        )
+        neighbor_list = torch.empty((2, max_pairs), dtype=torch.int32, device=device)
+        pair_offsets = torch.tensor([0, max_pairs], dtype=torch.int32, device=device)
+        pair_counts = torch.zeros(1, dtype=torch.int32, device=device)
+        neighbor_list_shifts = torch.empty(
+            (max_pairs, 3), dtype=torch.int32, device=device
+        )
+
+        result = cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell,
+            max_neighbors=64,
+            max_pairs=max_pairs,
+            format="coo",
+            rebuild_flags=torch.tensor([True], dtype=torch.bool, device=device),
+            return_state=True,
+            num_tiles=num_tiles,
+            tile_row_group=tile_row_group,
+            tile_col_group=tile_col_group,
+            neighbor_list=neighbor_list,
+            pair_offsets=pair_offsets,
+            pair_counts=pair_counts,
+            neighbor_list_shifts=neighbor_list_shifts,
+        )
+
+        assert len(result) == 7
+        assert result[0] is neighbor_list
+        assert result[1] is pair_offsets
+        assert result[2] is pair_counts
+        assert result[3] is neighbor_list_shifts
+        assert result[4] is num_tiles
+        assert result[5] is tile_row_group
+        assert result[6] is tile_col_group
+
+        pair_count = int(pair_counts.item())
+        segmented_pairs = {
+            tuple(entry)
+            for entry in torch.cat(
+                [neighbor_list[:, :pair_count].T, neighbor_list_shifts[:pair_count]],
+                dim=1,
+            )
+            .cpu()
+            .tolist()
+        }
+        compact_pairs = {
+            tuple(entry)
+            for entry in torch.cat([compact_list.T, compact_shifts], dim=1)
+            .cpu()
+            .tolist()
+        }
+        assert segmented_pairs == compact_pairs
+
+        reused = cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell,
+            max_neighbors=64,
+            max_pairs=max_pairs,
+            format="coo",
+            rebuild_flags=torch.tensor([False], dtype=torch.bool, device=device),
+            return_state=True,
+            num_tiles=result[4],
+            tile_row_group=result[5],
+            tile_col_group=result[6],
+            neighbor_list=result[0],
+            pair_offsets=result[1],
+            pair_counts=result[2],
+            neighbor_list_shifts=result[3],
+        )
+        assert all(reused[index] is result[index] for index in range(7))
+
+    def test_selective_coo_all_true_bootstrap_allocates_reusable_state(
+        self, device, dtype
+    ):
+        """An eager all-true selective COO call allocates reusable state."""
+        positions = torch.rand(64, 3, dtype=dtype, device=device) * 6.0
+        cell = _orthorhombic_cell(6.0, device, dtype)
+
+        result = cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell,
+            max_neighbors=64,
+            format="coo",
+            rebuild_flags=torch.ones(1, dtype=torch.bool, device=device),
+            return_state=True,
+        )
+
+        assert len(result) == 7
+        assert result[0].shape == (2, 64 * 64)
+        assert torch.equal(
+            result[1],
+            torch.tensor([0, 64 * 64], dtype=torch.int32, device=device),
+        )
+        assert result[2].dtype == torch.int32
+
+        reused = cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell,
+            max_neighbors=64,
+            format="coo",
+            rebuild_flags=torch.zeros(1, dtype=torch.bool, device=device),
+            return_state=True,
+            neighbor_list=result[0],
+            pair_offsets=result[1],
+            pair_counts=result[2],
+            neighbor_list_shifts=result[3],
+            num_tiles=result[4],
+            tile_row_group=result[5],
+            tile_col_group=result[6],
+        )
+
+        assert all(reused[index] is result[index] for index in range(7))
+
+    def test_selective_coo_overflow_raises(self, device, dtype):
+        """Selective COO reports counts beyond its fixed pair segment."""
+        positions = torch.zeros((64, 3), dtype=dtype, device=device)
+        cell = _orthorhombic_cell(8.0, device, dtype)
+        num_tiles, tile_row_group, tile_col_group, *_ = cluster_tile_neighbor_list(
+            positions, 2.0, cell, format="tile"
+        )
+
+        with pytest.raises(NeighborOverflowError):
+            cluster_tile_neighbor_list(
+                positions,
+                2.0,
+                cell,
+                max_neighbors=64,
+                max_pairs=1,
+                format="coo",
+                rebuild_flags=torch.tensor([True], dtype=torch.bool, device=device),
+                num_tiles=num_tiles,
+                tile_row_group=tile_row_group,
+                tile_col_group=tile_col_group,
+                neighbor_list=torch.empty((2, 1), dtype=torch.int32, device=device),
+                pair_offsets=torch.tensor([0, 1], dtype=torch.int32, device=device),
+                pair_counts=torch.zeros(1, dtype=torch.int32, device=device),
+                neighbor_list_shifts=torch.empty(
+                    (1, 3), dtype=torch.int32, device=device
+                ),
+            )
 
     def test_matrix_overflow_raises(self, device, dtype):
         positions = torch.rand(64, 3, dtype=dtype, device=device) * 8.0
diff --git a/test/neighbors/bindings/torch/test_neighborlist.py b/test/neighbors/bindings/torch/test_neighborlist.py
index 4d669c7c..d9d72566 100644
--- a/test/neighbors/bindings/torch/test_neighborlist.py
+++ b/test/neighbors/bindings/torch/test_neighborlist.py
@@ -32,6 +32,7 @@
 from nvalchemiops.torch.neighbors.cell_list import (
     cell_list,
 )
+from nvalchemiops.torch.neighbors.cluster_tile import cluster_tile_neighbor_list
 from nvalchemiops.torch.neighbors.naive import (
     naive_neighbor_list,
 )
@@ -46,6 +47,60 @@
 )
 
 
+class TestNeighborListClusterTileCompile:
+    """Validate compiled public cluster-tile dispatcher routes."""
+
+    @pytest.mark.slow
+    def test_compiled_unified_cluster_tile_is_rejected(self):
+        """Unified compiled cluster-tile cannot inspect tensor-valued PBC."""
+        if not torch.cuda.is_available():
+            pytest.skip("cluster_tile requires CUDA")
+        torch.manual_seed(42)
+        natom = 64
+        max_pairs = natom * 64
+        device = "cuda"
+        positions = torch.rand(natom, 3, dtype=torch.float32, device=device) * 6.0
+        cell = torch.eye(3, dtype=torch.float32, device=device).unsqueeze(0) * 6.0
+        pbc = torch.ones(3, dtype=torch.bool, device=device)
+        num_tiles, row, col, *_ = cluster_tile_neighbor_list(
+            positions,
+            2.0,
+            cell,
+            format="tile",
+        )
+        topology = torch.empty((2, max_pairs), dtype=torch.int32, device=device)
+        pair_offsets = torch.tensor([0, max_pairs], dtype=torch.int32, device=device)
+        pair_counts = torch.zeros(1, dtype=torch.int32, device=device)
+        shifts = torch.empty((max_pairs, 3), dtype=torch.int32, device=device)
+
+        @torch.compile(fullgraph=True)
+        def run(rebuild_flags):
+            return neighbor_list(
+                positions,
+                2.0,
+                cell=cell,
+                pbc=pbc,
+                method="cluster_tile",
+                return_neighbor_list=True,
+                rebuild_flags=rebuild_flags,
+                return_state=True,
+                num_tiles=num_tiles,
+                tile_row_group=row,
+                tile_col_group=col,
+                neighbor_list=topology,
+                pair_offsets=pair_offsets,
+                pair_counts=pair_counts,
+                neighbor_list_shifts=shifts,
+                max_neighbors=64,
+            )
+
+        with pytest.raises(
+            torch._dynamo.exc.Unsupported,
+            match="cannot validate tensor-valued pbc without synchronization",
+        ):
+            run(torch.ones(1, dtype=torch.bool, device=device))
+
+
 class TestNeighborListAutoSelection:
     """Test automatic method selection based on estimated neighbor density."""
 
@@ -573,6 +628,89 @@ def fake_batch_cluster_tile(positions, cutoff, cell, batch_ptr, **kwargs):
         torch.testing.assert_close(captured["cell"][0], cell)
         torch.testing.assert_close(captured["cell"][1], cell)
 
+    def test_explicit_cluster_tile_forwards_return_state(self, monkeypatch):
+        """The dispatcher forwards segmented COO state arguments to cluster_tile."""
+        captured = {}
+        pair_offsets = torch.tensor([0, 32], dtype=torch.int32)
+        pair_counts = torch.zeros(1, dtype=torch.int32)
+
+        def fake_cluster_tile(positions, cutoff, cell, **kwargs):
+            del positions, cutoff, cell
+            captured.update(kwargs)
+            return "cluster_tile"
+
+        monkeypatch.setattr(
+            neighbor_module, "cluster_tile_neighbor_list", fake_cluster_tile
+        )
+
+        result = neighbor_module.neighbor_list(
+            torch.zeros(32, 3, dtype=torch.float32),
+            2.0,
+            cell=torch.eye(3) * 8.0,
+            pbc=torch.ones(3, dtype=torch.bool),
+            method="cluster_tile",
+            return_neighbor_list=True,
+            rebuild_flags=torch.tensor([True], dtype=torch.bool),
+            return_state=True,
+            pair_offsets=pair_offsets,
+            pair_counts=pair_counts,
+        )
+
+        assert result == "cluster_tile"
+        assert captured["format"] == "coo"
+        assert captured["return_state"] is True
+        assert captured["pair_offsets"] is pair_offsets
+        assert captured["pair_counts"] is pair_counts
+
+    def test_explicit_batch_cluster_tile_forwards_return_state(self, monkeypatch):
+        """The dispatcher forwards state-return arguments to batch_cluster_tile."""
+        captured = {}
+
+        def fake_batch_cluster_tile(positions, cutoff, cell, batch_ptr, **kwargs):
+            del positions, cutoff, cell, batch_ptr
+            captured.update(kwargs)
+            return "batch_cluster_tile"
+
+        monkeypatch.setattr(
+            neighbor_module,
+            "batch_cluster_tile_neighbor_list",
+            fake_batch_cluster_tile,
+        )
+
+        result = neighbor_module.neighbor_list(
+            torch.zeros(32, 3, dtype=torch.float32),
+            2.0,
+            cell=torch.eye(3).unsqueeze(0) * 8.0,
+            pbc=torch.ones((1, 3), dtype=torch.bool),
+            batch_ptr=torch.tensor([0, 32], dtype=torch.int32),
+            method="batch_cluster_tile",
+            rebuild_flags=torch.tensor([True], dtype=torch.bool),
+            return_state=True,
+        )
+
+        assert result == "batch_cluster_tile"
+        assert captured["return_state"] is True
+
+    @pytest.mark.parametrize("method", [None, "naive", "batch_cell_list"])
+    def test_return_state_rejects_implicit_or_unsupported_method(self, method):
+        """State return requires an explicit cluster-tile method."""
+        match = (
+            "requires an explicit cluster_tile method"
+            if method is None
+            else "supported only by explicit cluster_tile methods"
+        )
+
+        with pytest.raises(ValueError, match=match):
+            neighbor_module.neighbor_list(
+                torch.zeros(32, 3, dtype=torch.float32),
+                2.0,
+                cell=torch.eye(3) * 8.0,
+                pbc=torch.ones(3, dtype=torch.bool),
+                method=method,
+                rebuild_flags=torch.tensor([True], dtype=torch.bool),
+                return_state=True,
+            )
+
     @pytest.mark.parametrize("dtype", [torch.float32, torch.float64])
     @pytest.mark.parametrize("device", ["cpu", "cuda"])
     def test_auto_select_batch_naive_dual_cutoff(self, dtype, device):
@@ -715,6 +853,45 @@ def test_explicit_batch_cluster_tile_no_cell_raises(self, dtype):
                 return_neighbor_list=True,
             )
 
+    @pytest.mark.parametrize("dtype", [torch.float32])
+    def test_explicit_cluster_tile_rejects_partial_pbc(self, dtype):
+        """Explicit cluster-tile rejects any non-periodic dimension."""
+        positions = torch.zeros((64, 3), dtype=dtype)
+        cell = torch.eye(3, dtype=dtype) * 8.0
+        pbc = torch.tensor([True, True, False], dtype=torch.bool)
+
+        with pytest.raises(NotImplementedError, match="pbc with any False"):
+            neighbor_list(
+                positions,
+                cutoff=2.0,
+                cell=cell,
+                pbc=pbc,
+                method="cluster_tile",
+                return_neighbor_list=True,
+            )
+
+    @pytest.mark.parametrize("dtype", [torch.float32])
+    def test_explicit_batch_cluster_tile_rejects_partial_pbc(self, dtype):
+        """Explicit batch cluster-tile rejects any non-periodic dimension."""
+        positions = torch.zeros((128, 3), dtype=dtype)
+        cell = torch.eye(3, dtype=dtype).expand(2, -1, -1).contiguous() * 8.0
+        pbc = torch.tensor(
+            [[True, True, True], [True, False, True]],
+            dtype=torch.bool,
+        )
+        batch_ptr = torch.tensor([0, 64, 128], dtype=torch.int32)
+
+        with pytest.raises(NotImplementedError, match="pbc with any False"):
+            neighbor_list(
+                positions,
+                cutoff=2.0,
+                cell=cell,
+                pbc=pbc,
+                batch_ptr=batch_ptr,
+                method="batch_cluster_tile",
+                return_neighbor_list=True,
+            )
+
     @pytest.mark.parametrize("dtype", [torch.float32])
     def test_explicit_cluster_tile_accepts_cutoff2(self, dtype):
         """method='cluster_tile' accepts cutoff2 on matrix output."""
diff --git a/test/neighbors/test_batch_cluster_tile_kernels.py b/test/neighbors/test_batch_cluster_tile_kernels.py
index d67b67fd..56f9ab6f 100644
--- a/test/neighbors/test_batch_cluster_tile_kernels.py
+++ b/test/neighbors/test_batch_cluster_tile_kernels.py
@@ -30,6 +30,7 @@
     TILE_GROUP_SIZE,
     batch_build_cluster_tile_list,
     batch_query_cluster_tile,
+    batch_query_cluster_tile_coo,
 )
 
 pytestmark = pytest.mark.gpu
@@ -437,3 +438,148 @@ def _spread(x: torch.Tensor) -> torch.Tensor:
                 if j == N:  # sentinel padding
                     continue
                 assert bi_cpu[i] == bi_cpu[j], f"Cross-system pair ({i}, {j}) emitted"
+
+
+class TestBatchSegmentedCooPhysicalBounds:
+    """Ensure per-system segmented COO offsets respect physical storage."""
+
+    def test_oversized_system_segment_preserves_canary(self, device):
+        """An invalid batched segment cannot write beyond max_pairs."""
+        natom_per_system = TILE_GROUP_SIZE
+        natom = 2 * natom_per_system
+        max_pairs = 16
+        backing_capacity = 2048
+        sorted_positions = torch.zeros(
+            natom,
+            dtype=torch.float32,
+            device=device,
+        )
+        sorted_atom_index = torch.arange(natom, dtype=torch.int32, device=device)
+        cell_batch = (
+            torch.eye(
+                3,
+                dtype=torch.float32,
+                device=device,
+            ).repeat(2, 1, 1)
+            * 8.0
+        )
+        inv_cell_batch = torch.linalg.inv(cell_batch).contiguous()
+        num_tiles = torch.tensor([2], dtype=torch.int32, device=device)
+        tile_row_group = torch.tensor([0, 1], dtype=torch.int32, device=device)
+        tile_col_group = torch.tensor([0, 1], dtype=torch.int32, device=device)
+        tile_system = torch.tensor([0, 1], dtype=torch.int32, device=device)
+        tile_offsets = torch.tensor([0, 1, 2], dtype=torch.int32, device=device)
+        tile_counts = torch.ones(2, dtype=torch.int32, device=device)
+        rebuild_flags = torch.ones(2, dtype=torch.bool, device=device)
+        pair_offsets = torch.tensor(
+            [0, 8, backing_capacity],
+            dtype=torch.int32,
+            device=device,
+        )
+        pair_counts = torch.zeros(2, dtype=torch.int32, device=device)
+        pair_counter = torch.zeros(1, dtype=torch.int32, device=device)
+        coo_list = torch.full(
+            (backing_capacity, 2),
+            -77,
+            dtype=torch.int32,
+            device=device,
+        )
+        coo_shifts = torch.full(
+            (backing_capacity, 3),
+            -77,
+            dtype=torch.int32,
+            device=device,
+        )
+
+        batch_query_cluster_tile_coo(
+            sorted_atom_index=wp.from_torch(
+                sorted_atom_index,
+                dtype=wp.int32,
+                return_ctype=True,
+            ),
+            sorted_pos_x=wp.from_torch(
+                sorted_positions,
+                dtype=wp.float32,
+                return_ctype=True,
+            ),
+            sorted_pos_y=wp.from_torch(
+                sorted_positions,
+                dtype=wp.float32,
+                return_ctype=True,
+            ),
+            sorted_pos_z=wp.from_torch(
+                sorted_positions,
+                dtype=wp.float32,
+                return_ctype=True,
+            ),
+            cell_batch=wp.from_torch(
+                cell_batch,
+                dtype=wp.mat33f,
+                return_ctype=True,
+            ),
+            inv_cell_batch=wp.from_torch(
+                inv_cell_batch,
+                dtype=wp.mat33f,
+                return_ctype=True,
+            ),
+            num_tiles=wp.from_torch(num_tiles, dtype=wp.int32, return_ctype=True),
+            tile_row_group=wp.from_torch(
+                tile_row_group,
+                dtype=wp.int32,
+                return_ctype=True,
+            ),
+            tile_col_group=wp.from_torch(
+                tile_col_group,
+                dtype=wp.int32,
+                return_ctype=True,
+            ),
+            tile_system=wp.from_torch(tile_system, dtype=wp.int32, return_ctype=True),
+            cutoff=2.0,
+            natom=natom,
+            max_pairs=max_pairs,
+            pair_counter=wp.from_torch(
+                pair_counter,
+                dtype=wp.int32,
+                return_ctype=True,
+            ),
+            coo_list=wp.from_torch(coo_list, dtype=wp.int32, return_ctype=True),
+            coo_shifts=wp.from_torch(
+                coo_shifts,
+                dtype=wp.int32,
+                return_ctype=True,
+            ),
+            wp_dtype=wp.float32,
+            device=device,
+            n_tiles=2,
+            rebuild_flags=wp.from_torch(
+                rebuild_flags,
+                dtype=wp.bool,
+                return_ctype=True,
+            ),
+            tile_offsets=wp.from_torch(
+                tile_offsets,
+                dtype=wp.int32,
+                return_ctype=True,
+            ),
+            tile_counts=wp.from_torch(
+                tile_counts,
+                dtype=wp.int32,
+                return_ctype=True,
+            ),
+            pair_offsets=wp.from_torch(
+                pair_offsets,
+                dtype=wp.int32,
+                return_ctype=True,
+            ),
+            pair_counts=wp.from_torch(
+                pair_counts,
+                dtype=wp.int32,
+                return_ctype=True,
+            ),
+        )
+
+        assert int(pair_counts[0].item()) > 0
+        assert int(pair_counts[1].item()) == 0
+        assert torch.all(coo_list[:8] != -77)
+        assert torch.all(coo_list[max_pairs:] == -77)
+        assert torch.all(coo_shifts[max_pairs:] == -77)
diff --git a/test/neighbors/test_cluster_tile_kernels.py b/test/neighbors/test_cluster_tile_kernels.py
index 626eb7b3..f3ca57f7 100644
--- a/test/neighbors/test_cluster_tile_kernels.py
+++ b/test/neighbors/test_cluster_tile_kernels.py
@@ -31,8 +31,12 @@
     TILE_GROUP_SIZE,
     build_cluster_tile_list,
     query_cluster_tile,
+    query_cluster_tile_coo,
+)
+from nvalchemiops.neighbors.cluster_tile.kernels import (
+    _bbox_distance_sq,
+    _coo_segment_is_valid,
 )
-from nvalchemiops.neighbors.cluster_tile.kernels import _bbox_distance_sq
 from nvalchemiops.neighbors.cluster_tile.launchers import (
     _compute_morton,
     _permute_gather_soa,
@@ -54,6 +58,18 @@ def _bbox_distance_sq_test_kernel(
     out[i] = _bbox_distance_sq(d[i], rg_ext[i], cg_ext[i])
 
 
+@wp.kernel(enable_backward=False)
+def _coo_segment_is_valid_test_kernel(
+    starts: wp.array(dtype=wp.int32),
+    stops: wp.array(dtype=wp.int32),
+    max_pairs: wp.int32,
+    out: wp.array(dtype=wp.int32),
+) -> None:
+    """Evaluate segmented COO interval validation for test values."""
+    i = wp.tid()
+    out[i] = _coo_segment_is_valid(starts[i], stops[i], max_pairs)
+
+
 def _mat33f_from_torch(mat: torch.Tensor):
     """Zero-copy view a ``(1, 3, 3)`` or ``(3, 3)`` torch tensor as a
     ``wp.array(dtype=wp.mat33f, shape=(1,))``.  Cluster-tile kernels read
@@ -124,6 +140,141 @@ def test_overlap_and_gap_cases(self, device):
         torch.testing.assert_close(out, expected)
 
 
+class TestSegmentedCooPhysicalBounds:
+    """Ensure caller offsets cannot escape physical COO output storage."""
+
+    def test_interval_validation_rejects_negative_reversed_and_oversized(
+        self,
+        device,
+    ):
+        """Only intervals inside the physical COO range are valid."""
+        starts = torch.tensor([0, -1, 8, 0], dtype=torch.int32, device=device)
+        stops = torch.tensor([8, 8, 0, 9], dtype=torch.int32, device=device)
+        out = torch.empty(4, dtype=torch.int32, device=device)
+
+        wp.launch(
+            kernel=_coo_segment_is_valid_test_kernel,
+            dim=4,
+            inputs=[
+                wp.from_torch(starts, dtype=wp.int32, return_ctype=True),
+                wp.from_torch(stops, dtype=wp.int32, return_ctype=True),
+                8,
+                wp.from_torch(out, dtype=wp.int32, return_ctype=True),
+            ],
+            device=device,
+        )
+
+        assert torch.equal(
+            out.cpu(),
+            torch.tensor([1, 0, 0, 0], dtype=torch.int32),
+        )
+
+    def test_single_segment_requires_full_physical_interval(self, device):
+        """A single-system segment must span the full physical COO buffer."""
+        natom = TILE_GROUP_SIZE
+        max_pairs = 8
+        backing_capacity = 1024
+        positions = torch.zeros((natom, 3), dtype=torch.float32, device=device)
+        sorted_atom_index = torch.arange(natom, dtype=torch.int32, device=device)
+        num_tiles = torch.ones(1, dtype=torch.int32, device=device)
+        tile_row_group = torch.zeros(1, dtype=torch.int32, device=device)
+        tile_col_group = torch.zeros(1, dtype=torch.int32, device=device)
+        cell = torch.eye(3, dtype=torch.float32, device=device) * 8.0
+        inv_cell = torch.linalg.inv(cell).contiguous()
+        pair_counter = torch.zeros(1, dtype=torch.int32, device=device)
+        coo_list = torch.full(
+            (backing_capacity, 2),
+            -77,
+            dtype=torch.int32,
+            device=device,
+        )
+        coo_shifts = torch.full(
+            (backing_capacity, 3),
+            -77,
+            dtype=torch.int32,
+            device=device,
+        )
+        rebuild_flags = torch.ones(1, dtype=torch.bool, device=device)
+        pair_offsets = torch.tensor(
+            [0, max_pairs // 2],
+            dtype=torch.int32,
+            device=device,
+        )
+        pair_counts = torch.zeros(1, dtype=torch.int32, device=device)
+
+        query_cluster_tile_coo(
+            sorted_atom_index=wp.from_torch(
+                sorted_atom_index,
+                dtype=wp.int32,
+                return_ctype=True,
+            ),
+            sorted_pos_x=wp.from_torch(
+                positions[:, 0],
+                dtype=wp.float32,
+                return_ctype=True,
+            ),
+            sorted_pos_y=wp.from_torch(
+                positions[:, 1],
+                dtype=wp.float32,
+                return_ctype=True,
+            ),
+            sorted_pos_z=wp.from_torch(
+                positions[:, 2],
+                dtype=wp.float32,
+                return_ctype=True,
+            ),
+            num_tiles=wp.from_torch(num_tiles, dtype=wp.int32, return_ctype=True),
+            tile_row_group=wp.from_torch(
+                tile_row_group,
+                dtype=wp.int32,
+                return_ctype=True,
+            ),
+            tile_col_group=wp.from_torch(
+                tile_col_group,
+                dtype=wp.int32,
+                return_ctype=True,
+            ),
+            cell=_mat33f_from_torch(cell),
+            inv_cell=_mat33f_from_torch(inv_cell),
+            cutoff=2.0,
+            natom=natom,
+            max_pairs=max_pairs,
+            pair_counter=wp.from_torch(
+                pair_counter,
+                dtype=wp.int32,
+                return_ctype=True,
+            ),
+            coo_list=wp.from_torch(coo_list, dtype=wp.int32, return_ctype=True),
+            coo_shifts=wp.from_torch(
+                coo_shifts,
+                dtype=wp.int32,
+                return_ctype=True,
+            ),
+            wp_dtype=wp.float32,
+            device=device,
+            n_tiles=1,
+            rebuild_flags=wp.from_torch(
+                rebuild_flags,
+                dtype=wp.bool,
+                return_ctype=True,
+            ),
+            pair_offsets=wp.from_torch(
+                pair_offsets,
+                dtype=wp.int32,
+                return_ctype=True,
+            ),
+            pair_counts=wp.from_torch(
+                pair_counts,
+                dtype=wp.int32,
+                return_ctype=True,
+            ),
+        )
+
+        assert int(pair_counts.item()) == 0
+        assert torch.all(coo_list == -77)
+        assert torch.all(coo_shifts == -77)
+
+
 class TestComputeMortonLauncher:
     """Smoke + invariant tests for the fused Morton + identity-init kernel."""