diff --git a/image_analysis/geom_refinement/BeamCenterFromSpots.cpp b/image_analysis/geom_refinement/BeamCenterFromSpots.cpp
index baf96e118..20dc3d8c7 100644
--- a/image_analysis/geom_refinement/BeamCenterFromSpots.cpp
+++ b/image_analysis/geom_refinement/BeamCenterFromSpots.cpp
@@ -6,7 +6,9 @@
 #include "../../common/JFJochMath.h"   // PI
 
 #include <algorithm>
+#include <array>
 #include <cmath>
+#include <future>
 #include <numeric>
 #include <tuple>
 
@@ -888,19 +890,32 @@ FindBeamCenterFromSpotSymmetry(const DiffractionExperiment &experiment,
     // nothing here can say which is the crystal's, so the scatter of the frame pairs - which stays
     // small for either of them - is not the uncertainty and this is. The spindle is re-fitted from
     // each start, so a fit that depends on where the search began is part of what is reported.
+    //
+    // The four run at once. Each takes its own copy of the geometry and only reads the spots, so they
+    // share nothing, and what is wanted from them is a MAX - which does not care in what order they
+    // finish. This is the whole of the estimator's cost: Estimate() is a brute-force grid over the
+    // spindle, and asking for the uncertainty runs it five times.
+    const std::array<std::pair<float, float>, 4> starts{{{CONSISTENCY_START_PXL, 0.0f},
+                                                         {-CONSISTENCY_START_PXL, 0.0f},
+                                                         {0.0f, CONSISTENCY_START_PXL},
+                                                         {0.0f, -CONSISTENCY_START_PXL}}};
+    std::vector<std::future<float>> restarts;
+    restarts.reserve(starts.size());
+    for (const auto &[dx, dy]: starts)
+        restarts.push_back(std::async(std::launch::async, [&, dx = dx, dy = dy] {
+            DiffractionGeometry from = geom;
+            from.BeamX_pxl(estimate->beam_x_pxl + dx).BeamY_pxl(estimate->beam_y_pxl + dy);
+            SpindleEstimate again;
+            if (const auto other = Estimate(from, *goniometer, frame_angle_deg, spots,
+                                            spindle_estimate ? &again : nullptr))
+                return std::hypot(other->beam_x_pxl - estimate->beam_x_pxl,
+                                  other->beam_y_pxl - estimate->beam_y_pxl);
+            return 0.0f;
+        }));
+
     float spread = 0.0f;
-    for (const auto &[dx, dy]: {std::pair<float, float>{CONSISTENCY_START_PXL, 0.0f},
-                                {-CONSISTENCY_START_PXL, 0.0f},
-                                {0.0f, CONSISTENCY_START_PXL},
-                                {0.0f, -CONSISTENCY_START_PXL}}) {
-        DiffractionGeometry from = geom;
-        from.BeamX_pxl(estimate->beam_x_pxl + dx).BeamY_pxl(estimate->beam_y_pxl + dy);
-        SpindleEstimate again;
-        if (const auto other = Estimate(from, *goniometer, frame_angle_deg, spots,
-                                        spindle_estimate ? &again : nullptr))
-            spread = std::max(spread, std::hypot(other->beam_x_pxl - estimate->beam_x_pxl,
-                                                 other->beam_y_pxl - estimate->beam_y_pxl));
-    }
+    for (auto &restart: restarts)
+        spread = std::max(spread, restart.get());
     estimate->sigma_pxl = std::max(estimate->sigma_pxl, spread);
     return estimate;
 }