diff --git a/RELEASING.md b/RELEASING.md
index 278692a7..99c5113a 100644
--- a/RELEASING.md
+++ b/RELEASING.md
@@ -1,3 +1,5 @@
+Note that this file is going to be AI Generated. If you're working on this, have your agent do the job. 
+
 # Releasing pele
 
 Release tested commits from `main` using matching version tags and GitHub Releases.
@@ -60,6 +62,54 @@ rattler-build build -r conda-recipe/recipe.yaml -c conda-forge \
   --variant c_stdlib=sysroot --variant c_stdlib_version=2.17
 ```
 
+## Deployment gotchas
+
+### Floating point arithmetic differs across platforms
+
+A fixed random seed does not guarantee identical numerical results across
+platforms. CPU instructions, compilers and math libraries can round intermediate
+operations differently. Non-convex optimization can amplify differences in the
+last bits and reach different stationary points, even when the initial coordinates
+and random-generator state agree exactly.
+
+For example, the tested macOS ARM NumPy 2.5.4 wheel for CPython 3.14 uses a
+fused multiply-add to compute
+`x1*x1 + x2*x2` in its legacy Gaussian conversion. Fusing the multiply and add
+rounds that combined operation once; separate operations round the product first.
+Replaying this arithmetic reproduced all 54 recorded macOS Gaussian values,
+including six that differed from Linux despite identical generator states.
+[NumPy's conversion source](https://github.com/numpy/numpy/blob/main/numpy/random/src/legacy/legacy-distributions.c#L18-L38)
+shows the expression involved.
+
+The [seeded LJ18 investigation](https://github.com/pele-python/pele/actions/runs/38177436370)
+also recorded a difference of one unit in the last place (1 ULP) in the first
+native minimization step. The resulting macOS saddle had exactly one unstable physical
+mode and a force RMS of `8.08e-5`, satisfying the `1e-4` stopping tolerance.
+Its transverse curvature search nevertheless exhausted 100 iterations with an
+RMS of `1.96e-5`, above its separate `1e-5` tolerance. At the same fixed coordinates,
+that search converged after 103 iterations when given a larger budget. A small
+iteration cap can therefore reject a valid candidate on one platform while
+passing on another. The [trace on PR #183](https://github.com/pele-python/pele/actions/runs/38178661821)
+reproduced this ARM result; its [Intel macOS wheel trace](https://github.com/pele-python/pele/actions/runs/38178661825)
+found the same cause, with the fixed-coordinate searches needing up to 117 iterations.
+
+A later [Intel macOS Python 3.12 trajectory trace](https://github.com/pele-python/pele/actions/runs/38184087581)
+reused one wheel across two runners and repeated the original seeded accounting
+fixture. All 44 starting configurations agreed bit for bit, but four searches
+failed: three after 21 genuine positive-curvature decisions, and one after
+exhausting 100 outer iterations with two unstable modes and a large force RMS.
+None of the 1,064 recorded primary snapshots misclassified the curvature sign,
+and evaluation accounting remained exact. In one positive-curvature failure,
+independent replay confirmed a fresh force RMS of `4.85e-4`, above `1e-4`, and
+48 positive projected physical curvatures, with the lowest `+0.10566`.
+These were legitimate failed searches, unlike the valid saddle rejected by the
+100-iteration cap. The trace reproduces the abort pattern; the earlier CI
+failure's exact coordinates were not recorded. Accounting tests therefore start
+near the existing known saddle, with a physical perturbation that still requires
+refinement, and retain their success and exact-count assertions. A loosely
+quenched random configuration does not guarantee a successful global
+transition-state search within a finite iteration budget.
+
 ## Configure publishing accounts
 
 Create accounts on [PyPI](https://pypi.org/) and [TestPyPI](https://test.pypi.org/).
diff --git a/pele/landscape/tests/test_double_ended_connect.py b/pele/landscape/tests/test_double_ended_connect.py
index 3c1a83c4..5698139d 100644
--- a/pele/landscape/tests/test_double_ended_connect.py
+++ b/pele/landscape/tests/test_double_ended_connect.py
@@ -45,10 +45,24 @@ def test2(self):
         # Allow enough search iterations across compiler and runtime variants.
         connect = DoubleEndedConnect(
             m1, m2, pot, mindist, db, verbosity=10,
-            local_connect_params={"NEBparams": {"iter_density": 100}},
+            local_connect_params={"NEBparams": {"iter_density": 100},
+                                  "tsSearchParams": {"index_seed": 0}},
         )
         connect.connect()
         self.assertTrue(connect.success())
+        for transition_state in db.transition_states():
+            coords = transition_state.coords
+            xyz = coords.reshape(-1, 3)
+            xyz = xyz - xyz.mean(axis=0)
+            rigid = np.column_stack(
+                [np.tile(axis, len(xyz)) for axis in np.eye(3)]
+                + [np.cross(axis, xyz).ravel() for axis in np.eye(3)]
+            )
+            physical = np.linalg.qr(rigid, mode="complete")[0][:, 6:]
+            hessian = pot.getHessian(coords)
+            spectrum = np.linalg.eigvalsh(physical.T @ hessian @ physical)
+            self.assertLess(spectrum[0], 0.0)
+            self.assertGreater(spectrum[1], 0.0)
 
         path = connect.returnPath()
 
diff --git a/pele/systems/tests/test_frozen_lj.py b/pele/systems/tests/test_frozen_lj.py
index 5a9dcd78..48715260 100644
--- a/pele/systems/tests/test_frozen_lj.py
+++ b/pele/systems/tests/test_frozen_lj.py
@@ -1,5 +1,7 @@
 import unittest
 
+import numpy as np
+
 from pele.systems.ljcluster_frozen import LJClusterFrozen
 from pele.systems import LJCluster
 from pele.thermodynamics import get_thermodynamic_information
@@ -7,6 +9,8 @@
 
 class TestLJCluster(unittest.TestCase):
     def setUp(self):
+        self.addCleanup(np.random.set_state, np.random.get_state())
+        np.random.seed(0)
         self.prepare_system(frozen_atoms=[0, 2, 4, 6])
 
     def prepare_system(self, frozen_atoms=None):
@@ -56,6 +60,7 @@ def test_basinhopping(self):
             bh.run(1)
 
     def test_connect(self):
+        self.system.params.double_ended_connect.local_connect_params.tsSearchParams.index_seed = 0
         db = self.make_database()
         min1, min2 = db.minima()[:2]
         connect = self.system.get_double_ended_connect(min1, min2, db)
diff --git a/pele/transition_states/_transition_state_refinement.py b/pele/transition_states/_transition_state_refinement.py
index 3d19f178..ad25ea70 100644
--- a/pele/transition_states/_transition_state_refinement.py
+++ b/pele/transition_states/_transition_state_refinement.py
@@ -4,6 +4,7 @@
 from pele.optimize import Result
 from pele.optimize import mylbfgs
 from pele.transition_states import FindLowestEigenVector
+from pele.transition_states._orthogopt import orthogopt
 from pele.transition_states._dimer_translator import _DimerTranslator
 from pele.transition_states._transverse_walker import _TransverseWalker
 from pele.utils.hessian import get_smallest_eig
@@ -48,6 +49,9 @@ class FindTransitionState:
                                 rotational symmetry
             orthogZeroEigs=None : the vector is unchanged
 
+    index_seed : int or None
+        Seed for the private random generator used by the additional-mode screen.
+        None uses fresh entropy without consuming the global NumPy random stream.
     lowestEigenvectorQuenchParams : dict
         these parameters are passed to the quench routine for he lowest
         eigenvector search
@@ -96,6 +100,12 @@ class FindTransitionState:
     if the tolerance for tangent space search is lower than the total tolerance,
     then it will never finish
 
+    A converged candidate is screened for additional unstable physical modes using
+    two iterative searches from independent random starts. This uses linear memory;
+    a passing screen is not an exhaustive certificate of Hessian index. The default
+    secondary-search tolerance bounds the gradient norm rather than a dimension-
+    dependent RMS. An explicitly supplied quench tolerance retains its RMS meaning.
+
     See Also
     --------
     findTransitionState : function wrapper for this class
@@ -126,8 +136,11 @@ def __init__(
         check_negative=True,
         invert_gradient=False,
         hessian_diagonalization=False,
+        index_seed=None,
     ):
         self.pot = pot
+        self.index_seed = index_seed
+        self._index_random = np.random.default_rng(index_seed)
         self.coords = np.copy(coords)
         self.nfev = 0
         self.tol = tol
@@ -224,6 +237,7 @@ def params(cls, obj=None):
         params[
             "lowestEigenvectorQuenchParams"
         ] = obj.lowestEigenvectorQuenchParams.copy()
+        params["index_seed"] = obj.index_seed
         params["tol"] = obj.tol
         params["nsteps"] = obj.nsteps
         params["nfail_max"] = obj.nfail_max
@@ -285,6 +299,7 @@ def run(self):
 
         self._compute_gradients(coords)
         iend = 0
+        index_screened = False
         for i in range(self.nsteps):
             iend = i
             # get the lowest eigenvalue and eigenvector
@@ -385,7 +400,29 @@ def run(self):
                     stepnum=i,
                 )
             if rms < self.tol:
-                break
+                self._getLowestEigenVector(coords, i, gradient=grad)
+                if (not self.leig_result.success
+                        or not np.isfinite(self.eigenval) or self.eigenval >= 0):
+                    break
+                transverse = self._get_lowest_transverse_eigenvector(coords, grad)
+                if not transverse.success:
+                    res.message.append("could not complete transverse stability screen")
+                    break
+                if transverse.eigenval > 0:
+                    index_screened = True
+                    break
+                res.message.append("additional negative curvature found")
+                if i + 1 == self.nsteps:
+                    break
+                direction = transverse.eigenvec
+                if np.dot(grad, direction) > 0:
+                    direction = -direction
+                coords = coords + min(0.01, self.maxstep_tangent) * direction
+                self._compute_gradients(coords)
+                if self._transverse_walker is not None:
+                    self.nfev += self._transverse_walker.get_result().nfev
+                    self._transverse_walker = None
+                continue
             if self.nfail >= self.nfail_max:
                 logger.warning(
                     "stopping findTransitionState.  too many failures in eigenvector search %s",
@@ -410,8 +447,12 @@ def run(self):
                     )
                     break
 
-        # done.  do one last eigenvector search because coords may have changed
-        self._getLowestEigenVector(coords, iend)
+        # A candidate passing the screen already has its final eigenpair checked.
+        if not index_screened:
+            self._getLowestEigenVector(coords, iend)
+        E = self.get_energy()
+        grad = self.get_gradient()
+        rms = np.linalg.norm(grad) * self.rmsnorm
 
         # print some data
         if self.verbosity > 0 or self.iprint > 0:
@@ -424,7 +465,7 @@ def run(self):
                 self.eigenval,
             )
 
-        success = True
+        success = index_screened
         # check if results make sense
         if self.eigenval >= 0.0:
             if self.verbosity > 2:
@@ -459,6 +500,89 @@ def run(self):
         res.nfev = self.nfev
         return res
 
+    def _get_lowest_transverse_eigenvector(self, coords, gradient):
+        """Screen for another unstable physical mode without forming a Hessian."""
+        orthogonalize = self.orthogZeroEigs
+        if orthogonalize == 0:
+            orthogonalize = orthogopt
+
+        def physical(vector):
+            vector = vector.copy()
+            if orthogonalize is not None:
+                vector = orthogonalize(vector, coords)
+            return vector
+
+        lowest = physical(self.eigenvec)
+        norm = np.linalg.norm(lowest)
+        if not np.isfinite(norm) or norm < 1e-12:
+            return Result(success=False)
+        lowest /= norm
+
+        def project(vector, position):
+            vector = physical(vector)
+            return vector - np.dot(vector, lowest) * lowest
+
+        params = self.lowestEigenvectorQuenchParams.copy()
+        # Unit eigenvectors need a dimension-independent gradient-norm threshold.
+        params.setdefault("tol", 1e-5 / np.sqrt(coords.size))
+        step = params.get("dx", 1e-6)
+
+        def hessian_product(vector):
+            _, plus = self.pot.getEnergyGradient(coords + step * vector)
+            _, minus = self.pot.getEnergyGradient(coords - step * vector)
+            self.nfev += 2
+            return (plus - minus) / (2 * step)
+
+        # Central differences also distinguish a negative mode from rigid-mode leakage.
+        primary_product = hessian_product(lowest)
+        primary_curvature = float(np.dot(lowest, primary_product))
+        primary_residual = np.linalg.norm(physical(primary_product) - primary_curvature * lowest)
+        if (not np.isfinite(primary_residual)
+                or not np.isfinite(primary_curvature)
+                or primary_curvature >= -primary_residual):
+            return Result(success=False)
+
+        # shortcut: two iterative starts can miss modes; use exhaustive spectral validation if certification is required.
+        # A local generator preserves the connection search's random sequence.
+        random = self._index_random
+        for _ in range(2):
+            guess = project(random.normal(size=coords.size), coords)
+            if np.linalg.norm(guess) < 1e-12:
+                for axis in range(coords.size):
+                    guess = np.zeros(coords.size)
+                    guess[axis] = 1.0
+                    guess = project(guess, coords)
+                    if np.linalg.norm(guess) >= 1e-12:
+                        break
+                else:
+                    # Only the checked unstable physical direction remains.
+                    return Result(success=True, eigenval=np.inf)
+            optimizer = FindLowestEigenVector(
+                coords, self.pot, eigenvec0=guess, orthogZeroEigs=project,
+                gradient=gradient, **params,
+            )
+            result = optimizer.run()
+            self.nfev += result.nfev
+            if not result.success or not np.isfinite(result.eigenval):
+                result.success = False
+                return result
+
+            direction = result.eigenvec
+            product = hessian_product(direction)
+            curvature = float(np.dot(direction, product))
+            residual = np.linalg.norm(project(product, coords) - curvature * direction)
+            coupling = 0.5 * (np.dot(lowest, product) + np.dot(direction, primary_product))
+            restricted = np.array([[primary_curvature, coupling], [coupling, curvature]])
+            result.eigenval = curvature
+            if (not np.all(np.isfinite(restricted)) or not np.isfinite(residual)
+                    or abs(curvature) <= residual):
+                result.success = False
+                return result
+            if curvature < 0:
+                result.success = np.linalg.eigvalsh(restricted)[-1] < 0
+                return result
+        return result
+
     def _get_lowest_eigenvector_RR(self, coords, gradient=None):
         """get the lowest eigenvector using Reyleigh Ritz minimization"""
         if "nsteps" in self.lowestEigenvectorQuenchParams:
@@ -614,6 +738,8 @@ def _minimizeTangentSpace(self, coords, energy=None, gradient=None):
 
     def _update_max_uphill_step(self, Fold, stepsize):
         """use a trust radius to update the maximum uphill step size"""
+        if stepsize == 0:
+            return
         Fnew = np.dot(self.eigenvec, self.get_gradient())
         # EPER=MIN(DABS(1.0D0-(FOBNEW-FOB)/(PSTEP*EVALMIN)),DABS(1.0D0-(-FOBNEW-FOB)/(PSTEP*EVALMIN)))
         a1 = 1.0 - (Fnew - Fold) / (stepsize * self.eigenval)
diff --git a/pele/transition_states/_transverse_walker.py b/pele/transition_states/_transverse_walker.py
index c5503ea3..c406b88b 100644
--- a/pele/transition_states/_transverse_walker.py
+++ b/pele/transition_states/_transverse_walker.py
@@ -111,10 +111,17 @@ def __init__(
             gradient=transverse_gradient,
             **minimizer_kwargs
         )
+        self._initial_H0 = self.walker.H0
 
     def update_eigenvec(self, eigenvec, eigenval):
-        """update the vecotr"""
+        """update the vector"""
+        old_vector = self.tspot.vector
         self.tspot.update_vector(eigenvec)
+        if not (np.array_equal(old_vector, self.tspot.vector) or
+                np.array_equal(old_vector, -self.tspot.vector)):
+            # Secant history from another projection can move along the excluded mode.
+            self.walker.reset()
+            self.walker.H0 = self._initial_H0
 
     def update_maxstep(self, maxstep):
         """update the maximum step size of the minimizer"""
diff --git a/pele/transition_states/tests/lj13_index_two.xyz b/pele/transition_states/tests/lj13_index_two.xyz
new file mode 100644
index 00000000..bf4aba4a
--- /dev/null
+++ b/pele/transition_states/tests/lj13_index_two.xyz
@@ -0,0 +1,15 @@
+13
+Recorded index-two candidate from the connection test
+LJ 1.2122346883612989 -0.6929443307790788 0.45087668210298265
+LJ 1.3493575684234813 0.30563112715902663 -0.035345239198378668
+LJ -0.99532473343649819 -0.57767505504430305 0.1432786388422807
+LJ 1.5784817299229652 -0.96002378359696505 -0.57925284369424734
+LJ 0.6186264596851152 0.62569254249252071 -0.82016815175178082
+LJ 1.2218868343284028 -0.16768392951581501 -1.3429024291359211
+LJ -0.26761773809639516 -0.96910258148931672 -0.6026477818204794
+LJ 0.65753316464245781 -0.41046513866176954 -0.45514205783651035
+LJ 0.64174869020119463 -1.4557093976098245 -0.14256864354711837
+LJ -0.27316555555765437 0.13206557898505247 -0.32581341182977808
+LJ 0.076108603335255084 -0.19908656303161248 -1.3538815467639136
+LJ 0.08024908140762875 -0.65300286334907787 0.4309591072748527
+LJ 2.1375687167827615 -0.01370930555883952 -0.73458586264199499
diff --git a/pele/transition_states/tests/test_find_transition_state_hef.py b/pele/transition_states/tests/test_find_transition_state_hef.py
index 7966ec1d..14b9cf02 100644
--- a/pele/transition_states/tests/test_find_transition_state_hef.py
+++ b/pele/transition_states/tests/test_find_transition_state_hef.py
@@ -24,6 +24,8 @@ def getEnergyGradient(self, coords):
 
 class TestFindTransitionState(unittest.TestCase):
     def setUp(self):
+        self.addCleanup(np.random.set_state, np.random.get_state())
+        np.random.seed(0)
         self.system = LJCluster(18)
         self.pot = PotWrapper(self.system.get_potential())
 
@@ -31,7 +33,7 @@ def make_dimer(self, x, **kwargs):
         return FindTransitionState(
             x,
             self.pot,
-            orthogZeroEigs=self.system.get_orthogonalize_to_zero_eigenvectors(),
+            index_seed=0, orthogZeroEigs=self.system.get_orthogonalize_to_zero_eigenvectors(),
             **kwargs
             #                                verbosity=10
         )
@@ -40,13 +42,32 @@ def test_params(self):
         params = FindTransitionState.params()
 
     def test1(self):
-        x = self.system.get_random_configuration()
+        path = os.path.join(os.path.dirname(__file__), "lj18_ts.xyz")
+        with open(path) as source:
+            x = read_xyz(source).coords.ravel()
+        perturbation = np.sin(np.arange(x.size) + 1.0)
+        perturbation = self.system.get_orthogonalize_to_zero_eigenvectors()(perturbation, x)
+        x = x + 0.01 * perturbation / np.linalg.norm(perturbation)
         dimer = self.make_dimer(x)
         res = dimer.run()
         self.assertTrue(res.success)
         self.assertEqual(res.nfev, self.pot.nfev)
         self.assertGreater(res.nfev, 1)
 
+        # Check physical curvature independently of the iterative index screen.
+        _, _, hessian = self.pot.pot.getEnergyGradientHessian(res.coords)
+        centered = res.coords.reshape(-1, 3)
+        centered = centered - centered.mean(axis=0)
+        rigid = np.column_stack([
+            np.tile(axis, len(centered)) for axis in np.eye(3)
+        ] + [np.cross(np.tile(axis, (len(centered), 1)), centered).ravel()
+             for axis in np.eye(3)])
+        self.assertEqual(np.linalg.matrix_rank(rigid), 6)
+        basis = np.linalg.qr(rigid, mode="complete")[0][:, 6:]
+        eigenvalues = np.linalg.eigvalsh(basis.T @ hessian @ basis)
+        self.assertLess(eigenvalues[0], -1e-5)
+        self.assertGreater(eigenvalues[1], 1e-5)
+
     def test2(self):
         # get the path of the file directory
         path = os.path.dirname(os.path.abspath(__file__))
@@ -75,8 +96,9 @@ def make_dimer(self, x, **kwargs):
         return FindTransitionState(
             x,
             self.pot,
-            orthogZeroEigs=self.system.get_orthogonalize_to_zero_eigenvectors(),
+            index_seed=0, orthogZeroEigs=self.system.get_orthogonalize_to_zero_eigenvectors(),
             invert_gradient=True,
+            **kwargs,
         )
 
 
@@ -87,11 +109,17 @@ def setUp(self):
         np.random.seed(0)
         self.system = LJCluster(18)
         self.pot = _PotWrapper(self.system.get_potential())
-        self.x = self.system.get_random_minimized_configuration(tol=10.0).coords
+        # Test evaluation accounting with actual refinement in a controlled saddle basin.
+        path = os.path.join(os.path.dirname(__file__), "lj18_ts.xyz")
+        with open(path) as source:
+            coords = read_xyz(source).coords.ravel()
+        perturbation = np.sin(np.arange(coords.size) + 1.0)
+        perturbation = self.system.get_orthogonalize_to_zero_eigenvectors()(perturbation, coords)
+        self.x = coords + 0.01 * perturbation / np.linalg.norm(perturbation)
 
     def do_check(self, **kwargs):
         self.pot.nfev = 0
-        opt = FindTransitionState(self.x, self.pot, **kwargs)
+        opt = FindTransitionState(self.x, self.pot, index_seed=0, **kwargs)
         ret = opt.run()
         self.assertEqual(ret.nfev, self.pot.nfev)
         self.assertTrue(ret.success)
diff --git a/pele/transition_states/tests/test_generalized_dimer.py b/pele/transition_states/tests/test_generalized_dimer.py
index 0f67e820..3195a01d 100644
--- a/pele/transition_states/tests/test_generalized_dimer.py
+++ b/pele/transition_states/tests/test_generalized_dimer.py
@@ -1,5 +1,6 @@
 import unittest
 import os
+import numpy as np
 
 from pele.systems import LJCluster
 from pele.transition_states._generalized_dimer import GeneralizedDimer
@@ -8,6 +9,7 @@
 
 class TestGeneralizedDimer(unittest.TestCase):
     def setUp(self):
+        np.random.seed(0)
         self.system = LJCluster(18)
         self.pot = self.system.get_potential()
 
diff --git a/pele/transition_states/tests/test_transition_state_index.py b/pele/transition_states/tests/test_transition_state_index.py
new file mode 100644
index 00000000..b78d2734
--- /dev/null
+++ b/pele/transition_states/tests/test_transition_state_index.py
@@ -0,0 +1,165 @@
+import unittest
+from unittest import mock
+
+import numpy as np
+
+from pele.potentials import BasePotential
+from pele.transition_states import FindTransitionState
+
+
+class _QuadraticPotential(BasePotential):
+    def __init__(self, hessian):
+        self.hessian = np.asarray(hessian, dtype=float)
+        self.nfev = 0
+
+    def getEnergyGradient(self, coords):
+        self.nfev += 1
+        gradient = self.hessian @ coords
+        return 0.5 * np.dot(coords, gradient), gradient
+
+    def getHessian(self, coords):
+        raise AssertionError("Default index validation must not form a dense Hessian")
+
+
+class TestTransitionStateIndexEdges(unittest.TestCase):
+    def test_rejects_rotated_index_two_with_stable_initial_rr_direction(self):
+        # Eigenvalues are -2, -1, +1; the first transverse probe points along +1.
+        pot = _QuadraticPotential([
+            [-2.0, 0.0, 0.0],
+            [0.0, -0.918372375740774, -0.3957172974185572],
+            [0.0, -0.3957172974185572, 0.918372375740774],
+        ])
+        ret = FindTransitionState(
+            np.zeros(3), pot, eigenvec0=np.array([1.0, 0.0, 0.0]),
+            index_seed=0, orthogZeroEigs=None, nsteps=1, verbosity=0,
+        ).run()
+        self.assertFalse(ret.success)
+
+    def test_random_start_finds_index_two_mode_missed_by_fixed_guesses(self):
+        # Spectrum -2, -1, +1, +1; both old fixed probes miss the -1 mode.
+        pot = _QuadraticPotential([
+            [-2.0, 0.0, 0.0, 0.0],
+            [0.0, -0.07552173356870506, -0.3733085508688484, 0.9246281380138601],
+            [0.0, -0.3733085508688484, 0.870426352344002, 0.3209340913541529],
+            [0.0, 0.9246281380138601, 0.3209340913541529, 0.20509538122470305],
+        ])
+        np.testing.assert_allclose(np.linalg.eigvalsh(pot.hessian),
+                                   [-2.0, -1.0, 1.0, 1.0], atol=1e-14)
+        ret = FindTransitionState(
+            np.zeros(4), pot, eigenvec0=np.array([1.0, 0.0, 0.0, 0.0]),
+            index_seed=1, orthogZeroEigs=None, nsteps=1, verbosity=0,
+        ).run()
+        self.assertEqual(ret.nfev, pot.nfev)
+        self.assertFalse(ret.success)
+
+    def test_weak_high_dimensional_mode_is_not_lost_to_rms_stopping(self):
+        dimension = 4096
+        hidden = np.random.default_rng(314159).normal(size=dimension)
+        hidden[0] = 0.0
+        hidden /= np.linalg.norm(hidden)
+
+        class RankTwoPotential(BasePotential):
+            def getEnergyGradient(self, coords):
+                # Spectrum -2, -0.01, +0.01... without allocating a Hessian.
+                gradient = 0.01 * coords - 0.02 * np.dot(hidden, coords) * hidden
+                gradient[0] -= 2.01 * coords[0]
+                return 0.5 * np.dot(coords, gradient), gradient
+
+        primary = np.zeros(dimension)
+        primary[0] = 1.0
+        search = FindTransitionState(
+            np.zeros(dimension), RankTwoPotential(), eigenvec0=primary,
+            index_seed=0, orthogZeroEigs=None, verbosity=0,
+        )
+        result = search._get_lowest_transverse_eigenvector(
+            np.zeros(dimension), np.zeros(dimension),
+        )
+        self.assertTrue(result.success)
+        self.assertLess(result.eigenval, 0.0)
+
+    def test_accepts_index_one_when_zero_modes_leave_no_transverse_direction(self):
+        def remove_zero_mode(vector, coords):
+            return np.array([vector[0], 0.0])
+
+        ret = FindTransitionState(
+            np.zeros(2), _QuadraticPotential([[-2.0, 0.0], [0.0, 0.0]]),
+            eigenvec0=np.array([1.0, 0.0]), index_seed=0, orthogZeroEigs=remove_zero_mode,
+            nsteps=1, verbosity=0,
+        ).run()
+        self.assertTrue(ret.success)
+        self.assertLess(ret.rms, 1e-4)
+
+    def test_accepts_one_dimensional_index_one_stationary_point(self):
+        ret = FindTransitionState(
+            np.zeros(1), _QuadraticPotential([[-2.0]]),
+            eigenvec0=np.ones(1), index_seed=0, orthogZeroEigs=None, nsteps=1, verbosity=0,
+        ).run()
+        self.assertTrue(ret.success)
+        self.assertAlmostEqual(ret.eigenval, -2.0)
+        self.assertLess(ret.rms, 1e-4)
+
+    def test_accepts_stiff_index_one_stationary_point(self):
+        # One negative mode and eleven positive modes with disparate stiffnesses.
+        curvatures = np.r_[-2.0, np.geomspace(1.0, 1000.0, 11)]
+        pot = _QuadraticPotential(np.diag(curvatures))
+        primary = np.zeros(curvatures.size)
+        primary[0] = 1.0
+        ret = FindTransitionState(
+            np.zeros(curvatures.size), pot, eigenvec0=primary,
+            index_seed=0, orthogZeroEigs=None, nsteps=1, verbosity=0,
+        ).run()
+        self.assertEqual(ret.nfev, pot.nfev)
+        self.assertTrue(ret.success, msg=repr(ret))
+
+    def test_success_requires_a_converged_returned_primary_eigenpair(self):
+        search = FindTransitionState(
+            np.zeros(1), _QuadraticPotential([[-2.0]]),
+            eigenvec0=np.ones(1), index_seed=0, orthogZeroEigs=None, nsteps=1, verbosity=0,
+        )
+        find_primary = search._get_lowest_eigenvector_RR
+        calls = 0
+
+        def fail_redundant_refresh(*args, **kwargs):
+            nonlocal calls
+            result = find_primary(*args, **kwargs)
+            calls += 1
+            if calls >= 3:
+                result.success = False
+            return result
+
+        with mock.patch.object(search, "_get_lowest_eigenvector_RR",
+                               side_effect=fail_redundant_refresh):
+            ret = search.run()
+        # Skipping the redundant refresh is valid; accepting its failed result is not.
+        self.assertTrue(not ret.success or search.leig_result.success)
+
+    def test_recovery_failure_returns_energy_and_gradient_at_returned_coords(self):
+        pot = _QuadraticPotential(np.diag([-2.0, -1.0, 1.0]))
+        search = FindTransitionState(
+            np.zeros(3), pot, eigenvec0=np.array([1.0, 0.0, 0.0]),
+            index_seed=0, orthogZeroEigs=None, demand_initial_negative_vec=True,
+            nsteps=2, verbosity=0,
+        )
+        search.npositive_max = 0
+        find_primary = search._getLowestEigenVector
+
+        def abort_after_recovery(coords, iteration, gradient=None):
+            overlap = find_primary(coords, iteration, gradient=gradient)
+            if iteration > 0:
+                search.eigenval = 1.0
+            return overlap
+
+        with mock.patch.object(search, "_getLowestEigenVector",
+                               side_effect=abort_after_recovery):
+            ret = search.run()
+        self.assertFalse(ret.success)
+        self.assertGreater(np.linalg.norm(ret.coords), 1e-5)
+        self.assertEqual(ret.nfev, pot.nfev)
+        energy, gradient = pot.getEnergyGradient(ret.coords)
+        self.assertAlmostEqual(ret.energy, energy, places=12)
+        np.testing.assert_allclose(ret.grad, gradient, rtol=1e-12, atol=1e-12)
+        self.assertAlmostEqual(ret.rms, np.linalg.norm(gradient) / np.sqrt(3.0), places=12)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/pele/transition_states/tests/test_transition_state_search.py b/pele/transition_states/tests/test_transition_state_search.py
index f6aae047..bc78f4d2 100644
--- a/pele/transition_states/tests/test_transition_state_search.py
+++ b/pele/transition_states/tests/test_transition_state_search.py
@@ -1,11 +1,13 @@
 import unittest
+from pathlib import Path
 
 
 import numpy as np
 import matplotlib.pyplot as plt
 
 from pele.transition_states import FindTransitionState, findTransitionState
-from pele.potentials import BasePotential
+from pele.potentials import BasePotential, LJ
+from pele.utils.xyz import read_xyz
 from pele.potentials.tests import _base_test
 from pele.potentials.tests._base_test import assert_arrays_almost_equal
 
@@ -75,6 +77,8 @@ def print_event(coords=None, **kwargs):
 
 class TestFindTransitionStateSimplePot(unittest.TestCase):
     def setUp(self):
+        self.addCleanup(np.random.set_state, np.random.get_state())
+        np.random.seed(0)
         self.pot = SimpleTSPot()
         self.x0 = np.array([0.1, 0.1])
         self.xts = np.zeros(self.x0.size)
@@ -85,7 +89,7 @@ def test1(self):
         opt = FindTransitionState(
             self.x0,
             self.pot,
-            orthogZeroEigs=None,
+            index_seed=0, orthogZeroEigs=None,
             # iprint=1,
             # verbosity=10, event=print_event,
             # tol=1e-3,
@@ -99,7 +103,7 @@ def test1(self):
         self.assertEqual(ret.nfev + 1, self.pot.nfev)
 
     def test_wrapper(self):
-        ret = findTransitionState(self.x0, self.pot, orthogZeroEigs=None)
+        ret = findTransitionState(self.x0, self.pot, index_seed=0, orthogZeroEigs=None)
         self.assertTrue(ret.success)
         assert_arrays_almost_equal(self, ret.coords, self.xts, places=3)
 
@@ -112,7 +116,7 @@ def event(**kwargs):
         opt = FindTransitionState(
             self.x0,
             self.pot,
-            orthogZeroEigs=None,
+            index_seed=0, orthogZeroEigs=None,
             tangentSpaceQuenchParams=dict(maxstep=1.0),
             event=event,
         )
@@ -126,7 +130,7 @@ def test_from_near_minimum(self):
         opt = FindTransitionState(
             x0,
             self.pot,
-            orthogZeroEigs=None,
+            index_seed=0, orthogZeroEigs=None,
             iprint=1,
             verbosity=10,  # event=print_event,
             # tol=1e-3,
@@ -145,7 +149,7 @@ def test_from_near_minimum_demand_negative_eigenvalue(self):
         opt = FindTransitionState(
             x0,
             self.pot,
-            orthogZeroEigs=None,
+            index_seed=0, orthogZeroEigs=None,
             demand_initial_negative_vec=True,
             iprint=1,
             verbosity=10,  # event=print_event,
@@ -166,7 +170,7 @@ def test1(self):
         opt = FindTransitionState(
             x0,
             pot,
-            orthogZeroEigs=None,
+            index_seed=0, orthogZeroEigs=None,
             iprint=1,
             verbosity=10,  # event=print_event,
             hessian_diagonalization=True,
@@ -176,5 +180,50 @@ def test1(self):
         print(ret)
 
 
+class IndexTwoPotential(BasePotential):
+    def getEnergyGradient(self, coords):
+        curvature = np.array([-2.0, -1.0, 1.0])
+        return 0.5 * np.dot(curvature * coords, coords), curvature * coords
+
+    def getHessian(self, coords):
+        raise AssertionError("Default index validation must not form a dense Hessian")
+
+
+class RecoverableIndexTwoPotential(BasePotential):
+    def getEnergyGradient(self, coords):
+        x, y, z = coords
+        return (-x*x + 0.25*(y*y - 1.0)**2 + 0.5*z*z,
+                np.array([-2.0*x, y*(y*y - 1.0), z]))
+
+
+class TestTransitionStateIndex(unittest.TestCase):
+    def test_rejects_index_two_stationary_point(self):
+        ret = FindTransitionState(
+            np.zeros(3), IndexTwoPotential(), eigenvec0=np.array([1., 0., 0.]),
+            index_seed=0, orthogZeroEigs=None, nsteps=1, verbosity=0,
+        ).run()
+        self.assertFalse(ret.success)
+
+    def test_rejects_recorded_lj13_index_two_saddle(self):
+        with Path(__file__).with_name("lj13_index_two.xyz").open() as source:
+            coords = read_xyz(source).coords.ravel()
+        pot = LJ()
+        values, vectors = np.linalg.eigh(pot.getHessian(coords))
+        self.assertLess(values[1], -1.0)
+        ret = FindTransitionState(
+            coords, pot, eigenvec0=vectors[:, 0], index_seed=0, nsteps=1, verbosity=0,
+        ).run()
+        self.assertFalse(ret.success)
+
+    def test_escapes_extra_unstable_mode(self):
+        ret = FindTransitionState(
+            np.zeros(3), RecoverableIndexTwoPotential(),
+            eigenvec0=np.array([1., 0., 0.]), index_seed=0, orthogZeroEigs=None, verbosity=0,
+        ).run()
+        self.assertTrue(ret.success)
+        self.assertAlmostEqual(abs(ret.coords[1]), 1.0, delta=1e-3)
+        self.assertLess(ret.rms, 1e-4)
+
+
 if __name__ == "__main__":
     unittest.main()
diff --git a/pele/transition_states/tests/test_transverse_walker.py b/pele/transition_states/tests/test_transverse_walker.py
index 1aa8ea4e..ba6c7236 100644
--- a/pele/transition_states/tests/test_transverse_walker.py
+++ b/pele/transition_states/tests/test_transverse_walker.py
@@ -23,7 +23,82 @@ def test(self):
         xnew = ret.coords
 
         dx = xnew - x0
-        self.assertLess(np.dot(dx, self.evec), 1e-6)
+        self.assertLess(abs(np.dot(dx, self.evec)), 1e-6)
+
+
+class TestTransverseWalkerHistory(unittest.TestCase):
+    def test_unchanged_plane_keeps_the_same_trajectory(self):
+        class QuadraticPotential:
+            def getEnergyGradient(self, coords):
+                gradient = np.array([1.0, 10.0, 100.0]) * coords
+                return 0.5 * np.dot(coords, gradient), gradient
+
+        pot = QuadraticPotential()
+        coords = np.array([1.0, -0.7, 0.4])
+        axis = np.array([0.0, 0.0, 1.0])
+        baseline = _TransverseWalker(coords, pot, axis, tol=1e-12)
+        expected = baseline.run(4).coords.copy()
+        for sign in (1.0, -1.0):
+            with self.subTest(sign=sign):
+                opt = _TransverseWalker(coords, pot, axis, tol=1e-12)
+                current = opt.run(3).coords.copy()
+                energy, gradient = pot.getEnergyGradient(current)
+                opt.update_eigenvec(sign * axis, -1.0)
+                opt.update_coords(current, energy, gradient)
+                np.testing.assert_allclose(opt.run(1).coords, expected,
+                                           rtol=0.0, atol=1e-14)
+
+    def test_rotation_restores_configured_initial_scale(self):
+        class QuadraticPotential:
+            def getEnergyGradient(self, coords):
+                gradient = np.array([1.0, 2.0, 3.0]) * coords
+                return 0.5 * np.dot(coords, gradient), gradient
+
+        pot = QuadraticPotential()
+        opt = _TransverseWalker(
+            np.array([1.0, -0.7, 0.4]), pot, np.array([0.0, 0.0, 1.0]),
+            tol=1e-12, H0=0.7, maxstep=10.0,
+        )
+        opt.run(3)
+        axis = np.array([1.0, 0.0, 0.0])
+        coords = np.array([-0.5, 1.0, 1.0])
+        energy, gradient = pot.getEnergyGradient(coords)
+        opt.update_eigenvec(axis, -1.0)
+        opt.update_coords(coords, energy, gradient)
+        fresh = _TransverseWalker(coords, pot, axis, tol=1e-12,
+                                  H0=0.7, maxstep=10.0)
+        np.testing.assert_allclose(opt.run(1).coords, fresh.run(1).coords,
+                                   rtol=0.0, atol=1e-14)
+
+    def test_motion_stays_in_updated_plane(self):
+        class QuadraticPotential:
+            hessian = np.array([[3.0, 1.0, 0.5],
+                                [1.0, 2.0, 0.25],
+                                [0.5, 0.25, 1.0]])
+
+            def getEnergyGradient(self, coords):
+                gradient = self.hessian @ coords
+                return 0.5 * np.dot(coords, gradient), gradient
+
+        pot = QuadraticPotential()
+        opt = _TransverseWalker(
+            np.array([1.0, -0.7, 0.4]), pot, np.array([0.0, 0.0, 1.0]),
+            tol=1e-12,
+        )
+        # Build valid secant history in the original transverse plane.
+        opt.run(3)
+        axis = np.array([1.0, 0.0, 0.0])
+        coords = np.array([-0.5, 1.0, 1.0])
+        energy, gradient = pot.getEnergyGradient(coords)
+        opt.update_eigenvec(axis, -1.0)
+        opt.update_coords(coords, energy, gradient)
+
+        start = coords.copy()
+        for _ in range(5):
+            previous = coords.copy()
+            coords = opt.run(1).coords.copy()
+            self.assertLess(abs(np.dot(coords - previous, axis)), 1e-14)
+        self.assertGreater(np.linalg.norm(coords - start), 1e-5)
 
 
 class TestTransverseWalker_NFEV(unittest.TestCase):