diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp
index f5fa161aef..4a88a259ca 100644
--- a/cpp/src/branch_and_bound/branch_and_bound.cpp
+++ b/cpp/src/branch_and_bound/branch_and_bound.cpp
@@ -314,12 +314,13 @@ branch_and_bound_t<i_t, f_t>::branch_and_bound_t(
 template <typename i_t, typename f_t>
 f_t branch_and_bound_t<i_t, f_t>::get_lower_bound()
 {
-  f_t lower_bound = lower_bound_numerical_.load();
+  f_t lower_bound = lower_bound_numerical_.load(std::memory_order::relaxed);
 
   if (bfs_worker_pool_.is_initialized()) {
     for (i_t i = 0; i < bfs_worker_pool_.size(); ++i) {
       if (bfs_worker_pool_[i]->is_active) {
-        lower_bound = std::min(lower_bound, bfs_worker_pool_[i]->lower_bound.load());
+        lower_bound =
+          std::min(lower_bound, bfs_worker_pool_[i]->lower_bound.load(std::memory_order::relaxed));
         lower_bound = std::min(lower_bound, bfs_worker_pool_[i]->node_queue.get_lower_bound());
       }
     }
@@ -404,10 +405,11 @@ void branch_and_bound_t<i_t, f_t>::report_heuristic(f_t obj, heuristics_origin_t
     log_line += std::format(" {:>8.2f}", toc(exploration_stats_.start_time));
     settings_.log.printf("%s\n", log_line.c_str());
   } else {
-    if (solving_root_relaxation_.load()) {
+    if (solving_root_relaxation_.load(std::memory_order::relaxed)) {
       f_t user_obj   = compute_user_objective(original_lp_, obj);
-      f_t user_lower = compute_user_objective(original_lp_, root_lp_current_lower_bound_.load());
-      f_t user_gap   = user_relative_gap(user_obj, user_lower);
+      f_t user_lower = compute_user_objective(
+        original_lp_, root_lp_current_lower_bound_.load(std::memory_order::relaxed));
+      f_t user_gap              = user_relative_gap(user_obj, user_lower);
       std::string user_gap_text = to_percentage(user_gap);
       settings_.log.print_format(
         "New solution from primal heuristics. Objective {:+.6e}. Gap {}. Time {:.2f}\n",
@@ -572,10 +574,10 @@ bool branch_and_bound_t<i_t, f_t>::set_solution_from_heuristics(const std::vecto
     if (is_feasible) {
       std::lock_guard lock(mutex_upper_);
       if (accept_new_incumbent(obj)) {
-        f_t current_upper_bound = upper_bound_.load();
+        f_t current_upper_bound = upper_bound_.load(std::memory_order::relaxed);
         upper_bound_            = std::min(current_upper_bound, obj);
         incumbent_.set_incumbent_solution(obj, crushed_solution);
-        if (current_upper_bound > upper_bound_.load()) {
+        if (current_upper_bound > upper_bound_.load(std::memory_order::relaxed)) {
           report_heuristic(obj, origin);
           success = true;
         }
@@ -827,7 +829,7 @@ void branch_and_bound_t<i_t, f_t>::repair_heuristic_solutions()
       if (is_feasible) {
         mutex_upper_.lock();
         if (accept_new_incumbent(repaired_obj)) {
-          upper_bound_ = std::min(upper_bound_.load(), repaired_obj);
+          upper_bound_ = std::min(upper_bound_.load(std::memory_order::relaxed), repaired_obj);
           incumbent_.set_incumbent_solution(repaired_obj, repaired_solution);
           report_heuristic(repaired_obj, heuristics_origin_t::HEURISTICS);
 
@@ -903,35 +905,38 @@ void branch_and_bound_t<i_t, f_t>::set_final_solution(mip_solution_t<i_t, f_t>&
     settings_.heuristic_preemption_callback();
   }
 
-  f_t user_obj         = compute_user_objective(original_lp_, upper_bound_.load());
+  f_t user_obj =
+    compute_user_objective(original_lp_, upper_bound_.load(std::memory_order::relaxed));
   f_t user_bound       = compute_user_objective(original_lp_, lower_bound);
   f_t gap              = std::abs(user_obj - user_bound);
   f_t gap_rel          = user_relative_gap(user_obj, user_bound);
   bool is_maximization = original_lp_.obj_scale < 0.0;
 
-  settings_.log.print_format("Explored {} nodes ({} simplex iterations) in {:.2f}s.",
-                             exploration_stats_.nodes_explored.load(),
-                             exploration_stats_.total_simplex_iters.load(),
-                             toc(exploration_stats_.start_time));
+  settings_.log.print_format(
+    "Explored {} nodes ({} simplex iterations) in {:.2f}s.",
+    exploration_stats_.nodes_explored.load(std::memory_order::relaxed),
+    exploration_stats_.total_simplex_iters.load(std::memory_order::relaxed),
+    toc(exploration_stats_.start_time));
 
-  if (exploration_stats_.orbital_fixing_nodes.load() > 0 ||
-      exploration_stats_.orbital_conflict_nodes.load() > 0) {
+  if (exploration_stats_.orbital_fixing_nodes.load(std::memory_order::relaxed) > 0 ||
+      exploration_stats_.orbital_conflict_nodes.load(std::memory_order::relaxed) > 0) {
     settings_.log.print_format(
       "Orbital fixing applied at {} nodes, {} total variable fixings, "
       "{} nodes with conflicting orbits\n",
-      exploration_stats_.orbital_fixing_nodes.load(),
-      exploration_stats_.orbital_fixings_applied.load(),
-      exploration_stats_.orbital_conflict_nodes.load());
+      exploration_stats_.orbital_fixing_nodes.load(std::memory_order::relaxed),
+      exploration_stats_.orbital_fixings_applied.load(std::memory_order::relaxed),
+      exploration_stats_.orbital_conflict_nodes.load(std::memory_order::relaxed));
   }
-  if (exploration_stats_.lexical_reduction_nodes.load() > 0) {
+  if (exploration_stats_.lexical_reduction_nodes.load(std::memory_order::relaxed) > 0) {
     settings_.log.print_format(
       "Lexical reduction applied at {} nodes, {} total variable fixings, {} nodes pruned\n",
-      exploration_stats_.lexical_reduction_nodes.load(),
-      exploration_stats_.lexical_reduction_fixings_applied.load(),
-      exploration_stats_.lexical_reduction_pruned_nodes.load());
+      exploration_stats_.lexical_reduction_nodes.load(std::memory_order::relaxed),
+      exploration_stats_.lexical_reduction_fixings_applied.load(std::memory_order::relaxed),
+      exploration_stats_.lexical_reduction_pruned_nodes.load(std::memory_order::relaxed));
   }
-  if (integer_pivots_.load() > 0) {
-    settings_.log.print_format("Number of integer pivots: {}\n", integer_pivots_.load());
+  if (integer_pivots_.load(std::memory_order::relaxed) > 0) {
+    settings_.log.print_format("Number of integer pivots: {}\n",
+                               integer_pivots_.load(std::memory_order::relaxed));
   }
 
   if (gap <= settings_.absolute_mip_gap_tol || gap_rel <= settings_.relative_mip_gap_tol) {
@@ -950,18 +955,12 @@ void branch_and_bound_t<i_t, f_t>::set_final_solution(mip_solution_t<i_t, f_t>&
     } else {
       settings_.log.printf("Optimal solution found.\n");
     }
-    if (settings_.heuristic_preemption_callback != nullptr) {
-      settings_.heuristic_preemption_callback();
-    }
   }
 
   if (solver_status_ == mip_status_t::UNSET) {
     if (exploration_stats_.nodes_explored > 0 && exploration_stats_.nodes_unexplored == 0 &&
         upper_bound_ == inf) {
       solver_status_ = mip_status_t::INFEASIBLE;
-      if (settings_.heuristic_preemption_callback != nullptr) {
-        settings_.heuristic_preemption_callback();
-      }
     }
   }
 
@@ -970,6 +969,7 @@ void branch_and_bound_t<i_t, f_t>::set_final_solution(mip_solution_t<i_t, f_t>&
     solution.objective     = incumbent_.objective;
     solution.has_incumbent = true;
   }
+
   solution.lower_bound        = lower_bound;
   solution.nodes_explored     = exploration_stats_.nodes_explored;
   solution.simplex_iterations = exploration_stats_.total_simplex_iters;
@@ -1002,7 +1002,7 @@ void branch_and_bound_t<i_t, f_t>::add_feasible_solution(const lp_problem_t<i_t,
   mutex_upper_.lock();
   if (accept_new_incumbent(leaf_objective)) {
     incumbent_.set_incumbent_solution(leaf_objective, sol.empty() ? leaf_solution : sol);
-    fetch_min(upper_bound_, leaf_objective);
+    upper_bound_.fetch_min(leaf_objective, std::memory_order::relaxed);
 
     char symbol = feasible_solution_symbol(thread_type, settings_.diving_settings.show_type);
     report(lp, symbol, leaf_objective, get_lower_bound(), leaf_depth, 0);
@@ -1202,10 +1202,11 @@ struct nondeterministic_policy_t : tree_update_policy_t<i_t, f_t> {
   void on_numerical_issue(mip_node_t<i_t, f_t>* node) override
   {
     if (worker->search_strategy == search_strategy_t::BEST_FIRST) {
-      fetch_min(bnb.lower_bound_numerical_, node->lower_bound);
+      bnb.lower_bound_numerical_.fetch_min(node->lower_bound, std::memory_order::relaxed);
       log.printf("LP returned numerical issue on node %d. Best bound set to %+10.6e.\n",
                  node->node_id,
-                 compute_user_objective(bnb.original_lp_, bnb.lower_bound_numerical_.load()));
+                 compute_user_objective(
+                   bnb.original_lp_, bnb.lower_bound_numerical_.load(std::memory_order::relaxed)));
     }
   }
 
@@ -1444,14 +1445,12 @@ std::pair<node_status_t, branch_direction_t> branch_and_bound_t<i_t, f_t>::updat
   if (lp_status == dual_status_t::DUAL_UNBOUNDED) {
     node_ptr->lower_bound = inf;
     policy.graphviz(search_tree, node_ptr, "infeasible", 0.0);
-    search_tree.update(node_ptr, node_status_t::INFEASIBLE);
     status = node_status_t::INFEASIBLE;
 
   } else if (lp_status == dual_status_t::CUTOFF) {
     f_t leaf_obj          = compute_objective(leaf_problem, leaf_solution.x);
     node_ptr->lower_bound = upper_bound;
     policy.graphviz(search_tree, node_ptr, "cut off", leaf_obj);
-    search_tree.update(node_ptr, node_status_t::FATHOMED);
     status = node_status_t::FATHOMED;
 
   } else if (lp_status == dual_status_t::OPTIMAL) {
@@ -1485,7 +1484,6 @@ std::pair<node_status_t, branch_direction_t> branch_and_bound_t<i_t, f_t>::updat
     if (num_frac == 0) {
       policy.handle_integer_solution(node_ptr, leaf_obj, leaf_solution.x);
       policy.graphviz(search_tree, node_ptr, "integer feasible", leaf_obj);
-      search_tree.update(node_ptr, node_status_t::INTEGER_FEASIBLE);
       status = node_status_t::INTEGER_FEASIBLE;
 
     } else if (node_ptr->lower_bound <= upper_bound + abs_fathom_tol) {
@@ -1510,12 +1508,11 @@ std::pair<node_status_t, branch_direction_t> branch_and_bound_t<i_t, f_t>::updat
                          worker->leaf_vstatus,
                          leaf_problem,
                          log);
-      search_tree.update(node_ptr, node_status_t::HAS_CHILDREN);
+      node_ptr->set_status(node_status_t::HAS_CHILDREN);
       status = node_status_t::HAS_CHILDREN;
 
     } else {
       policy.graphviz(search_tree, node_ptr, "fathomed", node_ptr->lower_bound);
-      search_tree.update(node_ptr, node_status_t::FATHOMED);
       status = node_status_t::FATHOMED;
     }
   } else if (lp_status == dual_status_t::TIME_LIMIT) {
@@ -1527,11 +1524,13 @@ std::pair<node_status_t, branch_direction_t> branch_and_bound_t<i_t, f_t>::updat
   } else {
     policy.on_numerical_issue(node_ptr);
     policy.graphviz(search_tree, node_ptr, "numerical", 0.0);
-    search_tree.update(node_ptr, node_status_t::NUMERICAL);
     status = node_status_t::NUMERICAL;
   }
 
   policy.on_node_completed(node_ptr, status, round_dir);
+
+  // `set_status` may free `node_ptr` for an inactive status, so it must be the last access.
+  if (inactive_status(status)) { node_ptr->set_status(status); }
   return {status, round_dir};
 }
 
@@ -1639,7 +1638,8 @@ dual_status_t branch_and_bound_t<i_t, f_t>::solve_node_lp(
     lp_settings.concurrent_halt = &node_concurrent_halt_;
   }
   lp_settings.set_log(false);
-  lp_settings.cut_off    = compute_node_cutoff(worker->leaf_problem, upper_bound_.load());
+  lp_settings.cut_off =
+    compute_node_cutoff(worker->leaf_problem, upper_bound_.load(std::memory_order::relaxed));
   lp_settings.inside_mip = 2;
   lp_settings.time_limit = settings_.time_limit - toc(exploration_stats_.start_time);
   if (lp_settings.time_limit <= 0.0) { return dual_status_t::TIME_LIMIT; }
@@ -1747,7 +1747,7 @@ dual_status_t branch_and_bound_t<i_t, f_t>::solve_node_lp(
                                         fractional,
                                         integer_pivot_work);
           if (num_integer_increased > 0) {
-            integer_pivots_.fetch_add(num_integer_increased, std::memory_order_release);
+            integer_pivots_.fetch_add(num_integer_increased, std::memory_order::relaxed);
           }
         }
       }
@@ -1765,7 +1765,7 @@ template <typename i_t, typename f_t>
 void branch_and_bound_t<i_t, f_t>::plunge_with(bfs_worker_t<i_t, f_t>* worker,
                                                mip_node_t<i_t, f_t>* start_node)
 {
-  assert(worker != nullptr && worker->is_active.load());
+  assert(worker != nullptr && worker->is_active.load(std::memory_order::relaxed));
   assert(start_node != nullptr);
 
   // Stack holds at most 2 entries: the preferred child + its sibling.
@@ -1812,9 +1812,9 @@ void branch_and_bound_t<i_t, f_t>::plunge_with(bfs_worker_t<i_t, f_t>* worker,
     // - The lower bound of the parent is lower or equal to its children
     worker->lower_bound = node_ptr->lower_bound;
 
-    if (node_ptr->lower_bound > upper_bound_.load()) {
+    if (node_ptr->lower_bound > upper_bound_.load(std::memory_order::relaxed)) {
       search_tree_.graphviz_node(settings_.log, node_ptr, "cutoff", node_ptr->lower_bound);
-      search_tree_.update(node_ptr, node_status_t::FATHOMED);
+      node_ptr->set_status(node_status_t::FATHOMED);
       worker->recompute_basis  = true;
       worker->recompute_bounds = true;
       --exploration_stats_.nodes_unexplored;
@@ -1831,7 +1831,7 @@ void branch_and_bound_t<i_t, f_t>::plunge_with(bfs_worker_t<i_t, f_t>* worker,
     if (((nodes_since_last_log >= 1000 || abs_gap < 10 * settings_.absolute_mip_gap_tol) &&
          time_since_last_log >= 1) ||
         (time_since_last_log > 30) || now > settings_.time_limit) {
-      if (exploration_stats_.report_solver_progress.exchange(false)) {
+      if (exploration_stats_.report_solver_progress.exchange(false, std::memory_order::relaxed)) {
         report(worker->leaf_problem,
                ' ',
                upper_bound_,
@@ -1981,7 +1981,7 @@ void branch_and_bound_t<i_t, f_t>::launch_bfs_worker(bfs_worker_t<i_t, f_t>* wor
   bfs_worker_t<i_t, f_t>* idle_worker = bfs_worker_pool_.pop_idle_worker();
   if (!idle_worker) return;
 
-  assert(idle_worker->is_active.load() == false);
+  assert(idle_worker->is_active.load(std::memory_order::relaxed) == false);
   assert(idle_worker->node_queue.best_first_queue_size() == 0);
 
   // Pre-emptively set the lower bound of the idle worker for the top of the heap
@@ -2026,10 +2026,12 @@ template <typename i_t, typename f_t>
 void branch_and_bound_t<i_t, f_t>::best_first_search_with(bfs_worker_t<i_t, f_t>* worker)
 {
   f_t lower_bound = get_lower_bound();
-  f_t user_obj    = compute_user_objective(worker->leaf_problem, upper_bound_.load());
-  f_t user_lower  = compute_user_objective(worker->leaf_problem, lower_bound);
-  f_t abs_gap     = compute_user_abs_gap(worker->leaf_problem, upper_bound_.load(), lower_bound);
-  f_t rel_gap     = user_relative_gap(user_obj, user_lower);
+  f_t user_obj =
+    compute_user_objective(worker->leaf_problem, upper_bound_.load(std::memory_order::relaxed));
+  f_t user_lower = compute_user_objective(worker->leaf_problem, lower_bound);
+  f_t abs_gap    = compute_user_abs_gap(
+    worker->leaf_problem, upper_bound_.load(std::memory_order::relaxed), lower_bound);
+  f_t rel_gap = user_relative_gap(user_obj, user_lower);
   f_t steal_chance =
     settings_.bnb_steal_chance >= 0 ? settings_.bnb_steal_chance : MIP_DEFAULT_STEAL_CHANCE;
   node_queue_t<i_t, f_t>& node_queue = worker->node_queue;
@@ -2103,11 +2105,11 @@ void branch_and_bound_t<i_t, f_t>::best_first_search_with(bfs_worker_t<i_t, f_t>
     if (!start_node) continue;
     worker->lower_bound = start_node->lower_bound;
 
-    if (upper_bound_.load() < start_node->lower_bound) {
+    if (upper_bound_.load(std::memory_order::relaxed) < start_node->lower_bound) {
       // This node was put on the heap earlier but its lower bound is now greater than the
       // current upper bound
       search_tree_.graphviz_node(settings_.log, start_node, "cutoff", start_node->lower_bound);
-      search_tree_.update(start_node, node_status_t::FATHOMED);
+      start_node->set_status(node_status_t::FATHOMED);
       --exploration_stats_.nodes_unexplored;
       continue;
     }
@@ -2115,10 +2117,12 @@ void branch_and_bound_t<i_t, f_t>::best_first_search_with(bfs_worker_t<i_t, f_t>
     plunge_with(worker, start_node);
 
     lower_bound = get_lower_bound();
-    user_obj    = compute_user_objective(worker->leaf_problem, upper_bound_.load());
-    user_lower  = compute_user_objective(worker->leaf_problem, lower_bound);
-    abs_gap     = compute_user_abs_gap(worker->leaf_problem, upper_bound_.load(), lower_bound);
-    rel_gap     = user_relative_gap(user_obj, user_lower);
+    user_obj =
+      compute_user_objective(worker->leaf_problem, upper_bound_.load(std::memory_order::relaxed));
+    user_lower = compute_user_objective(worker->leaf_problem, lower_bound);
+    abs_gap    = compute_user_abs_gap(
+      worker->leaf_problem, upper_bound_.load(std::memory_order::relaxed), lower_bound);
+    rel_gap = user_relative_gap(user_obj, user_lower);
 
     if (abs_gap <= settings_.absolute_mip_gap_tol || rel_gap <= settings_.relative_mip_gap_tol) {
       solver_status_ = mip_status_t::OPTIMAL;
@@ -2186,7 +2190,7 @@ void branch_and_bound_t<i_t, f_t>::dive_with(diving_worker_t<i_t, f_t>* worker,
 
     worker->lower_bound = node_ptr->lower_bound;
 
-    if (node_ptr->lower_bound > upper_bound_.load()) {
+    if (node_ptr->lower_bound > upper_bound_.load(std::memory_order::relaxed)) {
       worker->recompute_basis  = true;
       worker->recompute_bounds = true;
       continue;
@@ -2277,7 +2281,7 @@ bool branch_and_bound_t<i_t, f_t>::launch_diving_worker(bfs_worker_t<i_t, f_t>*
     return false;
   }
 
-  if (upper_bound_.load() < diving_worker->start_node.lower_bound) {
+  if (upper_bound_.load(std::memory_order::relaxed) < diving_worker->start_node.lower_bound) {
     diving_worker_pool_.return_worker_to_pool(diving_worker);
     return false;
   }
@@ -2307,7 +2311,8 @@ bool branch_and_bound_t<i_t, f_t>::launch_diving_worker(bfs_worker_t<i_t, f_t>*
   diving_worker->set_active();
   ++bfs_worker->active_diving_workers;
 
-  assert(bfs_worker->active_diving_workers.load() <= bfs_worker->max_diving_workers);
+  assert(bfs_worker->active_diving_workers.load(std::memory_order::relaxed) <=
+         bfs_worker->max_diving_workers);
 
 #pragma omp task affinity(*diving_worker) priority(CUOPT_DEFAULT_TASK_PRIORITY) default(none) \
   firstprivate(diving_worker)
@@ -2388,9 +2393,10 @@ void branch_and_bound_t<i_t, f_t>::solve_submip(diving_worker_t<i_t, f_t>* worke
   bool is_root_heuristic = submip_settings.inside_root_node;
   i_t submip_level       = settings_.submip_settings.level + 1;
   f_t user_lower         = compute_user_objective(worker->leaf_problem, get_lower_bound());
-  f_t user_obj           = compute_user_objective(worker->leaf_problem, upper_bound_.load());
-  f_t rel_gap            = user_relative_gap(user_obj, user_lower);
-  i_t explored           = exploration_stats_.nodes_explored;
+  f_t user_obj =
+    compute_user_objective(worker->leaf_problem, upper_bound_.load(std::memory_order::relaxed));
+  f_t rel_gap  = user_relative_gap(user_obj, user_lower);
+  i_t explored = exploration_stats_.nodes_explored;
 
   submip_settings.print_presolve_stats                     = false;
   submip_settings.num_threads                              = 1;
@@ -2436,7 +2442,7 @@ void branch_and_bound_t<i_t, f_t>::solve_submip(diving_worker_t<i_t, f_t>* worke
     submip_settings.time_limit,
     submip_settings.node_limit,
     submip_settings.branch_and_bound_simplex_iteration_limit,
-    exploration_stats_.total_simplex_iters.load(),
+    exploration_stats_.total_simplex_iters.load(std::memory_order::relaxed),
     submip_settings.relative_mip_gap_tol);
 
   // The `worker->leaf_problem` is directly converted to an `user_problem_t`, meaning that
@@ -2520,8 +2526,9 @@ void branch_and_bound_t<i_t, f_t>::solve_submip(diving_worker_t<i_t, f_t>* worke
 
   // Even if we do not have a valid incumbent now, the upper bound can still be set by the early
   // heuristics.
-  if (std::isfinite(upper_bound_.load())) {
-    const f_t user_upper    = compute_user_objective(worker->leaf_problem, upper_bound_.load());
+  if (std::isfinite(upper_bound_.load(std::memory_order::relaxed))) {
+    const f_t user_upper =
+      compute_user_objective(worker->leaf_problem, upper_bound_.load(std::memory_order::relaxed));
     const f_t submip_cutoff = compute_presolved_objective(submip_bnb.original_lp_, user_upper);
     submip_bnb.set_initial_upper_bound(submip_cutoff);
   }
@@ -2535,7 +2542,8 @@ void branch_and_bound_t<i_t, f_t>::solve_submip(diving_worker_t<i_t, f_t>* worke
   } else {
     // This should only be called by the main solver.
     submip_bnb.set_submip_halt_callback([this](f_t, f_t submip_lower_bound) {
-      f_t user_upper = compute_user_objective(this->original_lp_, this->upper_bound_.load());
+      f_t user_upper = compute_user_objective(this->original_lp_,
+                                              this->upper_bound_.load(std::memory_order::relaxed));
       bool is_cutoff = original_lp_.obj_scale > 0 ? submip_lower_bound > user_upper
                                                   : user_upper > submip_lower_bound;
       bool is_solver_running =
@@ -2603,7 +2611,7 @@ void branch_and_bound_t<i_t, f_t>::solve_submip(diving_worker_t<i_t, f_t>* worke
     log_prefix,
     mip_status_to_string(submip_status),
     submip_solution.simplex_iterations,
-    exploration_stats_.total_simplex_iters.load(),
+    exploration_stats_.total_simplex_iters.load(std::memory_order::relaxed),
     presolve_time,
     submip_time);
 
@@ -2934,9 +2942,9 @@ void branch_and_bound_t<i_t, f_t>::mutation(diving_worker_t<i_t, f_t>* worker,
   DEBUG_SUBMIP(
     "{}success={}, infeasible={}, calls={}, fixrate={:.4g} ({:.0f}/{}), max_fixrate={:.4g}\n",
     submip_settings.log.log_prefix,
-    mutation_stats_.total_success.load(),
-    mutation_stats_.total_infeasible.load(),
-    mutation_stats_.total_calls.load(),
+    mutation_stats_.total_success.load(std::memory_order::relaxed),
+    mutation_stats_.total_infeasible.load(std::memory_order::relaxed),
+    mutation_stats_.total_calls.load(std::memory_order::relaxed),
     fixrate,
     fixrate * integer_list.size(),
     integer_list.size(),
@@ -3125,7 +3133,7 @@ void branch_and_bound_t<i_t, f_t>::recursive_submip(
       DEBUG_SUBMIP("{}Round {}: max iteration reached! {}/{}",
                    submip_settings.log.log_prefix,
                    round,
-                   stats.total_simplex_iters.load(),
+                   stats.total_simplex_iters.load(std::memory_order::relaxed),
                    simplex_iter_limit)
       break;
     }
@@ -3146,12 +3154,12 @@ void branch_and_bound_t<i_t, f_t>::recursive_submip(
     node.lower_bound = leaf_obj;
 
     snap_to_lattice(worker->leaf_problem, &node, leaf_obj);
-    if (leaf_obj > upper_bound_.load()) {
+    if (leaf_obj > upper_bound_.load(std::memory_order::relaxed)) {
       DEBUG_SUBMIP("{}Round {}: reached cutoff point. obj={:.4g}. upper_bound={:.4g}",
                    submip_settings.log.log_prefix,
                    round,
                    leaf_obj,
-                   upper_bound_.load())
+                   upper_bound_.load(std::memory_order::relaxed))
       break;
     }
 
@@ -3163,7 +3171,7 @@ void branch_and_bound_t<i_t, f_t>::recursive_submip(
                    submip_settings.log.log_prefix,
                    round,
                    leaf_obj,
-                   upper_bound_.load())
+                   upper_bound_.load(std::memory_order::relaxed))
       break;
     }
 
@@ -3230,9 +3238,9 @@ void branch_and_bound_t<i_t, f_t>::recursive_submip(
     "{}success={}, infeasible={}, calls={}, fixrate={:.4g} ({:.0f}/{}), max_fixrate={:.4g}, "
     "min_fixrate={:.4g}\n",
     submip_settings.log.log_prefix,
-    submip_stats.total_success.load(),
-    submip_stats.total_infeasible.load(),
-    submip_stats.total_calls.load(),
+    submip_stats.total_success.load(std::memory_order::relaxed),
+    submip_stats.total_infeasible.load(std::memory_order::relaxed),
+    submip_stats.total_calls.load(std::memory_order::relaxed),
     fixrate,
     fixrate * num_integers,
     num_integers,
@@ -3695,7 +3703,7 @@ typename branch_and_bound_t<i_t, f_t>::cut_pass_action_t branch_and_bound_t<i_t,
                                   integer_pivot_work);
     root_integer_pivot_work_ += integer_pivot_work;
     if (num_integer_increased > 0) {
-      integer_pivots_.fetch_add(num_integer_increased, std::memory_order_release);
+      integer_pivots_.fetch_add(num_integer_increased, std::memory_order::relaxed);
     }
   }
   if (num_integer_increased > 0) {
@@ -3799,20 +3807,21 @@ typename branch_and_bound_t<i_t, f_t>::cut_pass_action_t branch_and_bound_t<i_t,
   }
 
   if (settings_.reduced_cost_strengthening >= 1 &&
-      (upper_bound_.load() < last_upper_bound ||
-       upper_bound_.load() <= reduced_cost_bounds.get_max_objective())) {
+      (upper_bound_.load(std::memory_order::relaxed) < last_upper_bound ||
+       upper_bound_.load(std::memory_order::relaxed) <= reduced_cost_bounds.get_max_objective())) {
     mutex_upper_.lock();
-    last_upper_bound              = upper_bound_.load();
+    last_upper_bound              = upper_bound_.load(std::memory_order::relaxed);
     std::vector<f_t> lower_bounds = original_lp_.lower;
     std::vector<f_t> upper_bounds = original_lp_.upper;
-    i_t new_bounds = reduced_cost_bounds.update_bounds_from_new_incumbent(original_lp_,
-                                                                          settings_,
-                                                                          root_relax_soln_,
-                                                                          root_objective_,
-                                                                          upper_bound_.load(),
-                                                                          var_types_,
-                                                                          lower_bounds,
-                                                                          upper_bounds);
+    i_t new_bounds                = reduced_cost_bounds.update_bounds_from_new_incumbent(
+      original_lp_,
+      settings_,
+      root_relax_soln_,
+      root_objective_,
+      upper_bound_.load(std::memory_order::relaxed),
+      var_types_,
+      lower_bounds,
+      upper_bounds);
     mutex_upper_.unlock();
     mutex_original_lp_.lock();
     original_lp_.lower = lower_bounds;
@@ -3937,7 +3946,7 @@ typename branch_and_bound_t<i_t, f_t>::cut_pass_action_t branch_and_bound_t<i_t,
       original_lp_, settings_, var_types_, root_objective_, root_relax_soln_.z, root_vstatus_);
     settings_.log.debug("New reduced cost objective %e (current %e)\n",
                         reduced_cost_bounds.get_max_objective(),
-                        upper_bound_.load());
+                        upper_bound_.load(std::memory_order::relaxed));
     if (settings_.primal_degenerate_pivots != 0) {
       pivot_to_improve_reduced_cost_strengthening(original_lp_,
                                                   lp_settings,
@@ -3955,7 +3964,7 @@ typename branch_and_bound_t<i_t, f_t>::cut_pass_action_t branch_and_bound_t<i_t,
     }
     settings_.log.debug("After pivoting: new reduced cost objective %e (current %e)\n",
                         reduced_cost_bounds.get_max_objective(),
-                        upper_bound_.load());
+                        upper_bound_.load(std::memory_order::relaxed));
   }
 
   // Refresh fractional info after re-solving with cuts; the pre-cut count is stale.
@@ -4010,18 +4019,20 @@ typename branch_and_bound_t<i_t, f_t>::cut_pass_action_t branch_and_bound_t<i_t,
   num_fractional = fractional_variables(settings_, root_relax_soln_.x, var_types_, fractional);
 
   if (num_fractional == 0) {
-    upper_bound_ = root_objective_;
     mutex_upper_.lock();
+    upper_bound_ = root_objective_;
     incumbent_.set_incumbent_solution(root_objective_, root_relax_soln_.x);
     mutex_upper_.unlock();
   }
-  f_t obj = upper_bound_.load();
+  f_t obj = upper_bound_.load(std::memory_order::relaxed);
   report(original_lp_, ' ', obj, root_objective_, 0, num_fractional);
 
-  f_t user_obj   = compute_user_objective(original_lp_, upper_bound_.load());
+  f_t user_obj =
+    compute_user_objective(original_lp_, upper_bound_.load(std::memory_order::relaxed));
   f_t user_lower = compute_user_objective(original_lp_, root_objective_);
   f_t rel_gap    = user_relative_gap(user_obj, user_lower);
-  f_t abs_gap    = compute_user_abs_gap(original_lp_, upper_bound_.load(), root_objective_);
+  f_t abs_gap    = compute_user_abs_gap(
+    original_lp_, upper_bound_.load(std::memory_order::relaxed), root_objective_);
   if (rel_gap < settings_.relative_mip_gap_tol || abs_gap < settings_.absolute_mip_gap_tol) {
     if (num_fractional == 0) { set_solution_at_root(solution, cut_info); }
     set_final_solution(solution, root_objective_);
@@ -4137,7 +4148,7 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
   lp_settings.inside_root_node                = true;
   lp_settings.concurrent_halt                 = get_root_concurrent_halt();
   lp_settings.dual_simplex_objective_callback = [this](f_t user_obj) {
-    root_lp_current_lower_bound_.store(user_obj);
+    root_lp_current_lower_bound_.store(user_obj, std::memory_order::relaxed);
   };
   std::vector<i_t> basic_list(original_lp_.num_rows);
   std::vector<i_t> nonbasic_list;
@@ -4308,7 +4319,7 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
     original_lp_, settings_, var_types_, root_objective_, root_relax_soln_.z, root_vstatus_);
   settings_.log.debug("New reduced cost objective %e (current %e)\n",
                       reduced_cost_bounds.get_max_objective(),
-                      upper_bound_.load());
+                      upper_bound_.load(std::memory_order::relaxed));
   if (settings_.primal_degenerate_pivots != 0) {
     pivot_to_improve_reduced_cost_strengthening(original_lp_,
                                                 lp_settings,
@@ -4326,7 +4337,7 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
   }
   settings_.log.debug("After pivoting: new reduced cost objective %e (current %e)\n",
                       reduced_cost_bounds.get_max_objective(),
-                      upper_bound_.load());
+                      upper_bound_.load(std::memory_order::relaxed));
 
   f_t pivot_out_integer_variables_start_time = tic();
   i_t num_integer_increased                  = 0;
@@ -4355,7 +4366,7 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
                                   integer_pivot_work);
     root_integer_pivot_work_ += integer_pivot_work;
     if (num_integer_increased > 0) {
-      integer_pivots_.fetch_add(num_integer_increased, std::memory_order_release);
+      integer_pivots_.fetch_add(num_integer_increased, std::memory_order::relaxed);
     }
   }
   if (num_integer_increased > 0) {
@@ -4454,10 +4465,12 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
       return mip_status_t::OPTIMAL;
     }
 
-    f_t user_obj   = compute_user_objective(original_lp_, upper_bound_.load());
+    f_t user_obj =
+      compute_user_objective(original_lp_, upper_bound_.load(std::memory_order::relaxed));
     f_t user_lower = compute_user_objective(original_lp_, root_objective_);
-    f_t abs_gap    = compute_user_abs_gap(original_lp_, upper_bound_.load(), root_objective_);
-    f_t rel_gap    = user_relative_gap(user_obj, user_lower);
+    f_t abs_gap    = compute_user_abs_gap(
+      original_lp_, upper_bound_.load(std::memory_order::relaxed), root_objective_);
+    f_t rel_gap = user_relative_gap(user_obj, user_lower);
 
     if (abs_gap <= settings_.absolute_mip_gap_tol || rel_gap <= settings_.relative_mip_gap_tol) {
       if (settings_.benchmark_info_ptr != nullptr) {
@@ -4592,18 +4605,19 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
   }
 
   if (settings_.reduced_cost_strengthening >= 2 &&
-      (upper_bound_.load() < last_upper_bound ||
-       upper_bound_.load() <= reduced_cost_bounds.get_max_objective())) {
+      (upper_bound_.load(std::memory_order::relaxed) < last_upper_bound ||
+       upper_bound_.load(std::memory_order::relaxed) <= reduced_cost_bounds.get_max_objective())) {
     std::vector<f_t> lower_bounds = original_lp_.lower;
     std::vector<f_t> upper_bounds = original_lp_.upper;
-    i_t num_changed = reduced_cost_bounds.update_bounds_from_new_incumbent(original_lp_,
-                                                                           settings_,
-                                                                           root_relax_soln_,
-                                                                           root_objective_,
-                                                                           upper_bound_.load(),
-                                                                           var_types_,
-                                                                           lower_bounds,
-                                                                           upper_bounds);
+    i_t num_changed               = reduced_cost_bounds.update_bounds_from_new_incumbent(
+      original_lp_,
+      settings_,
+      root_relax_soln_,
+      root_objective_,
+      upper_bound_.load(std::memory_order::relaxed),
+      var_types_,
+      lower_bounds,
+      upper_bounds);
     if (num_changed > 0) {
       settings_.log.printf("Updated %d integer bounds using reduced cost strengthening\n",
                            num_changed);
@@ -4762,12 +4776,12 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
 
       // We need to clear the queue and use the info in the search tree for the lower bound
       while (worker->node_queue.best_first_queue_size() > 0 &&
-             worker->node_queue.get_lower_bound() > upper_bound_.load()) {
+             worker->node_queue.get_lower_bound() > upper_bound_.load(std::memory_order::relaxed)) {
         mip_node_t<i_t, f_t>* start_node = worker->node_queue.pop();
         // This node was put on the heap earlier but its lower bound is now greater than the
         // current upper bound
         search_tree_.graphviz_node(settings_.log, start_node, "cutoff", start_node->lower_bound);
-        search_tree_.update(start_node, node_status_t::FATHOMED);
+        start_node->set_status(node_status_t::FATHOMED);
         --exploration_stats_.nodes_unexplored;
       }
 
@@ -4778,16 +4792,16 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
   }
 
   DEBUG_SUBMIP("RINS: success={}, infeasible={}, empty={}, calls={}",
-               rins_stats_.total_success.load(),
-               rins_stats_.total_infeasible.load(),
-               rins_stats_.total_empty.load(),
-               rins_stats_.total_calls.load());
+               rins_stats_.total_success.load(std::memory_order::relaxed),
+               rins_stats_.total_infeasible.load(std::memory_order::relaxed),
+               rins_stats_.total_empty.load(std::memory_order::relaxed),
+               rins_stats_.total_calls.load(std::memory_order::relaxed));
 
   DEBUG_SUBMIP("RENS: success={}, infeasible={}, empty={}, calls={}",
-               rens_stats_.total_success.load(),
-               rens_stats_.total_infeasible.load(),
-               rens_stats_.total_empty.load(),
-               rens_stats_.total_calls.load());
+               rens_stats_.total_success.load(std::memory_order::relaxed),
+               rens_stats_.total_infeasible.load(std::memory_order::relaxed),
+               rens_stats_.total_empty.load(std::memory_order::relaxed),
+               rens_stats_.total_calls.load(std::memory_order::relaxed));
 
   set_final_solution(solution, lower_bound);
   return solver_status_;
@@ -5097,7 +5111,7 @@ void branch_and_bound_t<i_t, f_t>::run_deterministic_bfs_loop(
       if (node->lower_bound > upper_bound) {
         worker.current_node = nullptr;
         worker.record_fathomed(node, node->lower_bound);
-        search_tree.update(node, node_status_t::FATHOMED);
+        node->set_status(node_status_t::FATHOMED);
         --exploration_stats_.nodes_unexplored;
         continue;
       }
@@ -5170,7 +5184,7 @@ void branch_and_bound_t<i_t, f_t>::deterministic_sync_callback()
     std::vector<uint64_t> state_data;
     state_data.push_back(static_cast<uint64_t>(exploration_stats_.nodes_explored));
     state_data.push_back(static_cast<uint64_t>(exploration_stats_.nodes_unexplored));
-    f_t ub = upper_bound_.load();
+    f_t ub = upper_bound_.load(std::memory_order::relaxed);
     f_t lb = deterministic_compute_lower_bound();
     state_data.push_back(std::bit_cast<uint64_t>(ub));
     state_data.push_back(std::bit_cast<uint64_t>(lb));
@@ -5210,7 +5224,7 @@ void branch_and_bound_t<i_t, f_t>::deterministic_sync_callback()
   }
 
   f_t lower_bound = deterministic_compute_lower_bound();
-  f_t upper_bound = upper_bound_.load();
+  f_t upper_bound = upper_bound_.load(std::memory_order::relaxed);
   f_t user_obj    = compute_user_objective(original_lp_, upper_bound);
   f_t user_lower  = compute_user_objective(original_lp_, lower_bound);
   f_t abs_gap     = compute_user_abs_gap(original_lp_, upper_bound, lower_bound);
@@ -5326,8 +5340,8 @@ node_status_t branch_and_bound_t<i_t, f_t>::solve_node_deterministic(
 
   if (!feasible) {
     node_ptr->lower_bound = std::numeric_limits<f_t>::infinity();
-    search_tree.update(node_ptr, node_status_t::INFEASIBLE);
     worker.record_infeasible(node_ptr);
+    node_ptr->set_status(node_status_t::INFEASIBLE);
     --exploration_stats_.nodes_unexplored;
     ++exploration_stats_.nodes_explored;
     worker.recompute_bounds_and_basis = true;
@@ -5410,14 +5424,14 @@ void branch_and_bound_t<i_t, f_t>::deterministic_process_worker_solutions(
                const queued_integer_solution_t<i_t, f_t>* b) { return *a < *b; });
 
   f_t deterministic_lower = deterministic_compute_lower_bound();
-  f_t current_upper       = upper_bound_.load();
+  f_t current_upper       = upper_bound_.load(std::memory_order::relaxed);
 
   for (const auto* sol : all_solutions) {
     if (sol->objective < current_upper) {
       f_t user_obj         = compute_user_objective(original_lp_, sol->objective);
       f_t user_lower       = compute_user_objective(original_lp_, deterministic_lower);
-      i_t nodes_explored   = exploration_stats_.nodes_explored.load();
-      i_t nodes_unexplored = exploration_stats_.nodes_unexplored.load();
+      i_t nodes_explored   = exploration_stats_.nodes_explored.load(std::memory_order::relaxed);
+      i_t nodes_unexplored = exploration_stats_.nodes_unexplored.load(std::memory_order::relaxed);
 
       search_strategy_t worker_type = get_worker_type(pool, sol->worker_id);
       report(original_lp_,
@@ -5430,7 +5444,7 @@ void branch_and_bound_t<i_t, f_t>::deterministic_process_worker_solutions(
 
       bool improved = false;
       if (accept_new_incumbent(sol->objective)) {
-        upper_bound_ = std::min(upper_bound_.load(), sol->objective);
+        upper_bound_ = std::min(upper_bound_.load(std::memory_order::relaxed), sol->objective);
         incumbent_.set_incumbent_solution(sol->objective, sol->solution);
         current_upper = sol->objective;
         improved      = true;
@@ -5593,7 +5607,7 @@ void branch_and_bound_t<i_t, f_t>::deterministic_sort_replay_events(
       // Process heuristic solution at its correct work unit timestamp position
       f_t new_upper = std::numeric_limits<f_t>::infinity();
       if (accept_new_incumbent(hsol.objective)) {
-        upper_bound_ = std::min(upper_bound_.load(), hsol.objective);
+        upper_bound_ = std::min(upper_bound_.load(std::memory_order::relaxed), hsol.objective);
         incumbent_.set_incumbent_solution(hsol.objective, hsol.solution);
         new_upper = hsol.objective;
       }
@@ -5620,14 +5634,14 @@ void branch_and_bound_t<i_t, f_t>::deterministic_sort_replay_events(
   deterministic_merge_pseudo_cost_updates(*deterministic_workers_);
 
   for (const auto& worker : *deterministic_workers_) {
-    fetch_min(lower_bound_numerical_, worker.local_lower_bound_ceiling);
+    lower_bound_numerical_.fetch_min(worker.local_lower_bound_ceiling, std::memory_order::relaxed);
   }
 }
 
 template <typename i_t, typename f_t>
 void branch_and_bound_t<i_t, f_t>::deterministic_prune_worker_nodes_vs_incumbent()
 {
-  f_t upper_bound = upper_bound_.load();
+  f_t upper_bound = upper_bound_.load(std::memory_order::relaxed);
 
   for (auto& worker : *deterministic_workers_) {
     // Check nodes in plunge stack - filter in place
@@ -5635,7 +5649,7 @@ void branch_and_bound_t<i_t, f_t>::deterministic_prune_worker_nodes_vs_incumbent
       std::deque<mip_node_t<i_t, f_t>*> surviving;
       for (auto* node : worker.plunge_stack) {
         if (node->lower_bound >= upper_bound) {
-          search_tree_.update(node, node_status_t::FATHOMED);
+          node->set_status(node_status_t::FATHOMED);
           --exploration_stats_.nodes_unexplored;
         } else {
           surviving.push_back(node);
@@ -5649,7 +5663,7 @@ void branch_and_bound_t<i_t, f_t>::deterministic_prune_worker_nodes_vs_incumbent
       std::vector<mip_node_t<i_t, f_t>*> surviving;
       for (auto* node : worker.backlog.data()) {
         if (node->lower_bound >= upper_bound) {
-          search_tree_.update(node, node_status_t::FATHOMED);
+          node->set_status(node_status_t::FATHOMED);
           --exploration_stats_.nodes_unexplored;
         } else {
           surviving.push_back(node);
@@ -5724,7 +5738,7 @@ template <typename i_t, typename f_t>
 f_t branch_and_bound_t<i_t, f_t>::deterministic_compute_lower_bound()
 {
   // Compute lower bound from BFS worker local structures only
-  f_t lower_bound = lower_bound_numerical_.load();
+  f_t lower_bound = lower_bound_numerical_.load(std::memory_order::relaxed);
 
   // Check all BFS worker queues
   for (const auto& worker : *deterministic_workers_) {
@@ -5746,7 +5760,7 @@ f_t branch_and_bound_t<i_t, f_t>::deterministic_compute_lower_bound()
 
   // Tree is exhausted
   if (lower_bound == std::numeric_limits<f_t>::infinity() && incumbent_.has_incumbent) {
-    lower_bound = upper_bound_.load();
+    lower_bound = upper_bound_.load(std::memory_order::relaxed);
   }
 
   return lower_bound;
@@ -5763,7 +5777,7 @@ void branch_and_bound_t<i_t, f_t>::deterministic_populate_diving_heap()
   const int num_diving                  = deterministic_diving_workers_->size();
   constexpr int target_nodes_per_worker = 10;
   const int target_total                = num_diving * target_nodes_per_worker;
-  f_t cutoff                            = upper_bound_.load();
+  f_t cutoff                            = upper_bound_.load(std::memory_order::relaxed);
 
   // Collect candidate nodes from BFS worker backlog heaps
   std::vector<std::pair<mip_node_t<i_t, f_t>*, f_t>> candidates;
diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp
index 58fe285b9c..4144eba7f4 100644
--- a/cpp/src/branch_and_bound/branch_and_bound.hpp
+++ b/cpp/src/branch_and_bound/branch_and_bound.hpp
@@ -157,7 +157,7 @@ class branch_and_bound_t {
   void set_initial_pseudocost(const pseudo_costs_t<i_t, f_t>& parent_pc,
                               const std::vector<i_t>& reduced_to_original);
 
-  f_t get_upper_bound() const { return upper_bound_.load(); }
+  f_t get_upper_bound() const { return upper_bound_.load(std::memory_order::relaxed); }
   bool has_solver_space_incumbent() const { return incumbent_.has_incumbent; }
 
   // Repair a low-quality solution from the heuristics.
diff --git a/cpp/src/branch_and_bound/mip_node.hpp b/cpp/src/branch_and_bound/mip_node.hpp
index 21d74c0e1d..c2bb020d52 100644
--- a/cpp/src/branch_and_bound/mip_node.hpp
+++ b/cpp/src/branch_and_bound/mip_node.hpp
@@ -80,8 +80,9 @@ class mip_node_t {
       objective_estimate(std::numeric_limits<f_t>::infinity()),
       packed_vstatus(0)
   {
-    children[0] = nullptr;
-    children[1] = nullptr;
+    children[0]  = nullptr;
+    children[1]  = nullptr;
+    num_children = 0;
   }
 
   mip_node_t(f_t root_lower_bound, const std::vector<simplex::variable_status_t>& basis)
@@ -96,8 +97,9 @@ class mip_node_t {
       objective_estimate(std::numeric_limits<f_t>::infinity()),
       packed_vstatus(compress_vstatus(basis))
   {
-    children[0] = nullptr;
-    children[1] = nullptr;
+    children[0]  = nullptr;
+    children[1]  = nullptr;
+    num_children = 0;
   }
 
   mip_node_t(const simplex::lp_problem_t<i_t, f_t>& problem,
@@ -126,6 +128,7 @@ class mip_node_t {
                                                                     : problem.upper[branch_var];
     children[0]      = nullptr;
     children[1]      = nullptr;
+    num_children     = 0;
   }
 
   void get_variable_bounds(std::vector<f_t>& lower,
@@ -177,54 +180,24 @@ class mip_node_t {
     children[1] = std::move(up_child);
     // When we add children we no longer need to store our basis
     packed_vstatus = {};
+    num_children   = 2;
   }
 
-  bool is_inactive() const
-  {
-    if (inactive_status(status)) { return true; }
-    if ((children[0] != nullptr && inactive_status(children[0]->status)) &&
-        (children[1] != nullptr && inactive_status(children[1]->status))) {
-      return true;
-    }
-    if (children[0] == nullptr && inactive_status(children[1]->status)) { return true; }
-    if (children[1] == nullptr && inactive_status(children[0]->status)) { return true; }
-    return false;
-  }
-
-  void update_bound()
-  {
-    if (children[0] != nullptr && children[1] != nullptr) {
-      if (inactive_status(children[0]->status) && inactive_status(children[1]->status)) {
-        lower_bound = std::min(children[0]->lower_bound, children[1]->lower_bound);
-      }
-    }
-    if (children[0] != nullptr && children[1] == nullptr) {
-      if (inactive_status(children[0]->status)) { lower_bound = children[0]->lower_bound; }
-    }
-    if (children[1] != nullptr && children[0] == nullptr) {
-      if (inactive_status(children[1]->status)) { lower_bound = children[1]->lower_bound; }
-    }
-  }
-
-  // outputs a stack containing inactive nodes in the tree that can be freed
-  void set_status(node_status_t node_status, std::vector<mip_node_t*>& stack)
+  void set_status(node_status_t node_status)
   {
     status = node_status;
-    if (inactive_status(status)) {
-      update_bound();
-      stack.push_back(this);
-      // Propagate to parent
-      mip_node_t* parent_ptr = parent;
-      while (parent_ptr != nullptr) {
-        if (parent_ptr->is_inactive()) {
-          parent_ptr->status = node_status_t::FATHOMED;
-          parent_ptr->update_bound();
-          stack.push_back(parent_ptr);
-        } else {
-          break;
-        }
-        parent_ptr = parent_ptr->parent;
-      }
+    if (!inactive_status(node_status)) { return; }
+
+    mip_node_t* parent_ptr = parent;
+    while (parent_ptr != nullptr) {
+      i_t old = parent_ptr->num_children.fetch_sub(1, std::memory_order::acq_rel);
+      if (old != 1) { return; }  // The sibling is still active.
+      parent_ptr->status = node_status_t::FATHOMED;
+      parent_ptr->lower_bound =
+        std::min(parent_ptr->children[0]->lower_bound, parent_ptr->children[1]->lower_bound);
+      parent_ptr->children[0].reset();
+      parent_ptr->children[1].reset();
+      parent_ptr = parent_ptr->parent;
     }
   }
 
@@ -302,8 +275,9 @@ class mip_node_t {
   f_t fractional_val;
   i_t integer_infeasible;
 
-  mip_node_t<i_t, f_t>* parent;
+  mip_node_t* parent;
   std::unique_ptr<mip_node_t> children[2];
+  omp_atomic_t<i_t> num_children;
 
   std::vector<uint8_t> packed_vstatus;
 
@@ -344,16 +318,6 @@ class mip_node_t {
   }
 };
 
-template <typename i_t, typename f_t>
-void remove_fathomed_nodes(std::vector<mip_node_t<i_t, f_t>*>& stack)
-{
-  for (int i = 0; i < stack.size(); ++i) {
-    for (int child = 0; child < 2; ++child) {
-      if (stack[i]->children[child] != nullptr) { stack[i]->children[child].reset(); }
-    }
-  }
-}
-
 template <typename i_t, typename f_t>
 class search_tree_t {
  public:
@@ -361,14 +325,6 @@ class search_tree_t {
 
   search_tree_t(mip_node_t<i_t, f_t>&& node) : root(std::move(node)), num_nodes(0) {}
 
-  void update(mip_node_t<i_t, f_t>* node_ptr, node_status_t status)
-  {
-    std::lock_guard<omp_mutex_t> lock(mutex);
-    std::vector<mip_node_t<i_t, f_t>*> stack;
-    node_ptr->set_status(status, stack);
-    remove_fathomed_nodes(stack);
-  }
-
   void branch(mip_node_t<i_t, f_t>* parent_node,
               const i_t branch_var,
               const f_t fractional_val,
@@ -377,7 +333,7 @@ class search_tree_t {
               const simplex::lp_problem_t<i_t, f_t>& original_lp,
               simplex::logger_t& log)
   {
-    i_t id = num_nodes.fetch_add(2);
+    i_t id = num_nodes.fetch_add(2, std::memory_order::relaxed);
 
     auto down_child = std::make_unique<mip_node_t<i_t, f_t>>(original_lp,
                                                              parent_node,
@@ -443,7 +399,6 @@ class search_tree_t {
   }
 
   mip_node_t<i_t, f_t> root;
-  omp_mutex_t mutex;
   omp_atomic_t<i_t> num_nodes;
 
   static constexpr bool write_graphviz = false;
diff --git a/cpp/src/branch_and_bound/node_queue.hpp b/cpp/src/branch_and_bound/node_queue.hpp
index 47c9cb22c7..c5114e863d 100644
--- a/cpp/src/branch_and_bound/node_queue.hpp
+++ b/cpp/src/branch_and_bound/node_queue.hpp
@@ -44,7 +44,7 @@ class heap_t {
     buffer.emplace_back(std::forward<Args>(args)...);
     std::push_heap(buffer.begin(), buffer.end(), comp);
     ++num_entries_;
-    assert(num_entries_.load() == buffer.size());
+    assert(num_entries_.load(std::memory_order::relaxed) == buffer.size());
   }
 
   T pop()
@@ -53,7 +53,7 @@ class heap_t {
     std::pop_heap(buffer.begin(), buffer.end(), comp);
     T node = std::move(buffer.back());
     buffer.pop_back();
-    assert(num_entries_.load() == buffer.size());
+    assert(num_entries_.load(std::memory_order::relaxed) == buffer.size());
     return node;
   }
 
@@ -154,7 +154,8 @@ class node_queue_t {
 
   f_t get_lower_bound()
   {
-    return best_first_heap_.empty() ? std::numeric_limits<f_t>::infinity() : lower_bound_.load();
+    return best_first_heap_.empty() ? std::numeric_limits<f_t>::infinity()
+                                    : lower_bound_.load(std::memory_order::relaxed);
   }
 
  private:
diff --git a/cpp/src/branch_and_bound/worker.hpp b/cpp/src/branch_and_bound/worker.hpp
index 2be97a8179..dcf2252a68 100644
--- a/cpp/src/branch_and_bound/worker.hpp
+++ b/cpp/src/branch_and_bound/worker.hpp
@@ -301,11 +301,11 @@ class diving_worker_t : public branch_and_bound_worker_t<i_t, f_t> {
   // Set this node inactive
   void set_inactive()
   {
-    if (!this->is_active.load()) { return; }
+    if (!this->is_active.load(std::memory_order::relaxed)) { return; }
     this->is_active = false;
 
     if (bfs_worker) {
-      assert(bfs_worker->active_diving_workers.load() > 0);
+      assert(bfs_worker->active_diving_workers.load(std::memory_order::relaxed) > 0);
       --bfs_worker->active_diving_workers;
     }
   }
diff --git a/cpp/src/branch_and_bound/worker_pool.hpp b/cpp/src/branch_and_bound/worker_pool.hpp
index 2360e8920d..2e67339eb4 100644
--- a/cpp/src/branch_and_bound/worker_pool.hpp
+++ b/cpp/src/branch_and_bound/worker_pool.hpp
@@ -53,7 +53,7 @@ class worker_pool_t {
       workers_[i]->symmetry_ptr = symmetry;
     }
 
-    is_initialized_ = true;
+    is_initialized_.store(true, std::memory_order::release);
   }
 
   WorkerType* pop_idle_worker()
@@ -65,7 +65,8 @@ class worker_pool_t {
       i_t idx = idle_workers_.front();
       idle_workers_.pop_front();
       num_idle_workers_--;
-      assert(idle_workers_.size() == static_cast<size_t>(num_idle_workers_.load()));
+      assert(idle_workers_.size() ==
+             static_cast<size_t>(num_idle_workers_.load(std::memory_order::relaxed)));
       assert(idx >= 0 && static_cast<size_t>(idx) < workers_.size());
       return workers_[idx].get();
     }
@@ -75,13 +76,14 @@ class worker_pool_t {
   {
     assert(worker != nullptr);
     worker->set_inactive();
-    assert(!worker->is_active.load());
+    assert(!worker->is_active.load(std::memory_order::relaxed));
 
-    if (!is_initialized_) return;
+    if (!is_initialized()) return;
 
     std::lock_guard lock(mutex_);
     assert(workers_[worker->worker_id].get() == worker);
-    assert(static_cast<size_t>(num_idle_workers_.load()) == idle_workers_.size());
+    assert(static_cast<size_t>(num_idle_workers_.load(std::memory_order::relaxed)) ==
+           idle_workers_.size());
     assert(idle_workers_.size() <= workers_.size());
 
     idle_workers_.push_back(worker->worker_id);
@@ -101,14 +103,14 @@ class worker_pool_t {
     return workers_[id].get();
   }
 
-  bool is_initialized() const { return is_initialized_; }
+  bool is_initialized() const { return is_initialized_.load(std::memory_order::acquire); }
 
   i_t num_idle() const { return num_idle_workers_; }
   i_t size() const { return workers_.size(); }
 
  private:
   std::vector<std::unique_ptr<WorkerType>> workers_;
-  bool is_initialized_ = false;
+  omp_atomic_t<bool> is_initialized_{false};
 
   omp_mutex_t mutex_;
   circular_deque_t<i_t> idle_workers_;
diff --git a/cpp/src/mip_heuristics/local_search/local_search.cu b/cpp/src/mip_heuristics/local_search/local_search.cu
index f4d2d90908..8d7b1d381c 100644
--- a/cpp/src/mip_heuristics/local_search/local_search.cu
+++ b/cpp/src/mip_heuristics/local_search/local_search.cu
@@ -111,13 +111,12 @@ void local_search_t<i_t, f_t>::start_cpufj_scratch_threads(population_t<i_t, f_t
         f_t obj, const std::vector<f_t>& h_vec, double /*work_units*/) {
         population.add_external_solution(h_vec, obj, solution_origin_t::CPUFJ);
         (void)problem_ptr;
-        if (obj < this->local_search_best_obj) {
+        if (obj < this->local_search_best_obj.fetch_min(obj, std::memory_order::relaxed)) {
           CUOPT_LOG_TRACE("******* New local search best obj %g, best overall %g",
                           problem_ptr->get_user_obj_from_solver_obj(obj),
                           problem_ptr->get_user_obj_from_solver_obj(
                             population.is_feasible() ? population.best_feasible().get_objective()
                                                      : std::numeric_limits<f_t>::max()));
-          this->local_search_best_obj = obj;
         }
       };
     counter++;
@@ -229,13 +228,12 @@ void local_search_t<i_t, f_t>::start_cpufj_lptopt_scratch_threads(
   scratch_cpu_fj_on_lp_opt->improvement_callback =
     [this, &population](f_t obj, const std::vector<f_t>& h_vec, double /*work_units*/) {
       population.add_external_solution(h_vec, obj, solution_origin_t::CPUFJ);
-      if (obj < this->local_search_best_obj) {
+      if (obj < this->local_search_best_obj.fetch_min(obj, std::memory_order::relaxed)) {
         CUOPT_LOG_DEBUG("******* New local search best obj %g, best overall %g",
                         context.problem_ptr->get_user_obj_from_solver_obj(obj),
                         context.problem_ptr->get_user_obj_from_solver_obj(
                           population.is_feasible() ? population.best_feasible().get_objective()
                                                    : std::numeric_limits<f_t>::max()));
-        this->local_search_best_obj = obj;
       }
     };
 
diff --git a/cpp/src/mip_heuristics/root_heuristics.hpp b/cpp/src/mip_heuristics/root_heuristics.hpp
index 784ba9a282..422b4d735a 100644
--- a/cpp/src/mip_heuristics/root_heuristics.hpp
+++ b/cpp/src/mip_heuristics/root_heuristics.hpp
@@ -305,8 +305,8 @@ struct root_heuristics_t {
     // On the first pass we use a thread to generate the clique table
     const i_t cut_generation     = cut_pass == 0 ? 2 : 1;
     const i_t persistent_workers = persistent_lanes_.size();
-    return std::max(i_t{0},
-                    max_workers_ - worker_count_->load() - persistent_workers - cut_generation);
+    const i_t worker_count       = worker_count_->load(std::memory_order::relaxed);
+    return std::max(i_t{0}, max_workers_ - worker_count - persistent_workers - cut_generation);
   }
 
   std::shared_ptr<cut_pass_heuristics_t<i_t, f_t>> create_new_cut_pass_heuristic(
diff --git a/cpp/src/utilities/omp_helpers.hpp b/cpp/src/utilities/omp_helpers.hpp
index 3fd31d169e..b40e2f1acc 100644
--- a/cpp/src/utilities/omp_helpers.hpp
+++ b/cpp/src/utilities/omp_helpers.hpp
@@ -87,58 +87,58 @@ class omp_atomic_t {
 
   T operator=(T new_val)
   {
-    store(new_val);
+    store(new_val, std::memory_order::relaxed);
     return new_val;
   }
 
-  operator T() const { return load(); }
-  T operator+=(T inc) { return fetch_add(inc) + inc; }
-  T operator-=(T inc) { return fetch_sub(inc) - inc; }
+  operator T() const { return load(std::memory_order::relaxed); }
+  T operator+=(T inc) { return fetch_add(inc, std::memory_order::relaxed) + inc; }
+  T operator-=(T inc) { return fetch_sub(inc, std::memory_order::relaxed) - inc; }
 
   // In theory, this should be enabled only for integers,
   // but it works for any numerical types.
-  T operator++() { return fetch_add(T(1)) + 1; }
-  T operator++(int) { return fetch_add(T(1)); }
-  T operator--() { return fetch_sub(T(1)) - 1; }
-  T operator--(int) { return fetch_sub(T(1)); }
+  T operator++() { return fetch_add(T(1), std::memory_order::relaxed) + 1; }
+  T operator++(int) { return fetch_add(T(1), std::memory_order::relaxed); }
+  T operator--() { return fetch_sub(T(1), std::memory_order::relaxed) - 1; }
+  T operator--(int) { return fetch_sub(T(1), std::memory_order::relaxed); }
 
   // Possible values for memory order: relaxed, acquire, seq_cst
-  T load(std::memory_order memory_order = std::memory_order::seq_cst) const
+  T load(std::memory_order memory_order) const
   {
     T res;
-    if (memory_order == std::memory_order::relaxed) {
-#pragma omp atomic read relaxed
+    if (memory_order == std::memory_order::seq_cst) {
+#pragma omp atomic read seq_cst
       res = val;
     } else if (memory_order == std::memory_order::acquire) {
 #pragma omp atomic read acquire
       res = val;
     } else {
-#pragma omp atomic read
+#pragma omp atomic read relaxed
       res = val;
     }
     return res;
   }
 
   // Possible values for memory order: relaxed, release, seq_cst
-  void store(T new_val, std::memory_order memory_order = std::memory_order::seq_cst)
+  void store(T new_val, std::memory_order memory_order)
   {
-    if (memory_order == std::memory_order::relaxed) {
-#pragma omp atomic write relaxed
+    if (memory_order == std::memory_order::seq_cst) {
+#pragma omp atomic write seq_cst
       val = new_val;
     } else if (memory_order == std::memory_order::release) {
 #pragma omp atomic write release
       val = new_val;
     } else {
-#pragma omp atomic write
+#pragma omp atomic write relaxed
       val = new_val;
     }
   }
 
-  T exchange(T other, std::memory_order memory_order = std::memory_order::seq_cst)
+  T exchange(T other, std::memory_order memory_order)
   {
     T old;
-    if (memory_order == std::memory_order::relaxed) {
-#pragma omp atomic capture relaxed
+    if (memory_order == std::memory_order::seq_cst) {
+#pragma omp atomic capture seq_cst
       {
         old = val;
         val = other;
@@ -162,7 +162,7 @@ class omp_atomic_t {
         val = other;
       }
     } else {
-#pragma omp atomic capture
+#pragma omp atomic capture relaxed
       {
         old = val;
         val = other;
@@ -171,11 +171,11 @@ class omp_atomic_t {
     return old;
   }
 
-  T fetch_add(T inc, std::memory_order memory_order = std::memory_order::seq_cst)
+  T fetch_add(T inc, std::memory_order memory_order)
   {
     T old;
-    if (memory_order == std::memory_order::relaxed) {
-#pragma omp atomic capture relaxed
+    if (memory_order == std::memory_order::seq_cst) {
+#pragma omp atomic capture seq_cst
       {
         old = val;
         val += inc;
@@ -199,7 +199,7 @@ class omp_atomic_t {
         val += inc;
       }
     } else {
-#pragma omp atomic capture
+#pragma omp atomic capture relaxed
       {
         old = val;
         val += inc;
@@ -208,7 +208,83 @@ class omp_atomic_t {
     return old;
   }
 
-  T fetch_sub(T inc) { return fetch_add(-inc); }
+  T fetch_sub(T inc, std::memory_order memory_order) { return fetch_add(-inc, memory_order); }
+
+  T fetch_min(T other, std::memory_order memory_order)
+  {
+    T old;
+
+    if (memory_order == std::memory_order::seq_cst) {
+#pragma omp atomic compare capture seq_cst
+      {
+        old = val;
+        if (other < val) { val = other; }
+      }
+    } else if (memory_order == std::memory_order::acquire) {
+#pragma omp atomic compare capture acquire
+      {
+        old = val;
+        if (other < val) { val = other; }
+      }
+    } else if (memory_order == std::memory_order::release) {
+#pragma omp atomic compare capture release
+      {
+        old = val;
+        if (other < val) { val = other; }
+      }
+    } else if (memory_order == std::memory_order::acq_rel) {
+#pragma omp atomic compare capture acq_rel
+      {
+        old = val;
+        if (other < val) { val = other; }
+      }
+    } else {
+#pragma omp atomic compare capture relaxed
+      {
+        old = val;
+        if (other < val) { val = other; }
+      }
+    }
+
+    return old;
+  }
+
+  T fetch_max(T other, std::memory_order memory_order)
+  {
+    T old;
+    if (memory_order == std::memory_order::seq_cst) {
+#pragma omp atomic compare capture seq_cst
+      {
+        old = val;
+        if (other > val) { val = other; }
+      }
+    } else if (memory_order == std::memory_order::acquire) {
+#pragma omp atomic compare capture acquire
+      {
+        old = val;
+        if (other > val) { val = other; }
+      }
+    } else if (memory_order == std::memory_order::release) {
+#pragma omp atomic compare capture release
+      {
+        old = val;
+        if (other > val) { val = other; }
+      }
+    } else if (memory_order == std::memory_order::acq_rel) {
+#pragma omp atomic compare capture acq_rel
+      {
+        old = val;
+        if (other > val) { val = other; }
+      }
+    } else {
+#pragma omp atomic compare capture relaxed
+      {
+        old = val;
+        if (other > val) { val = other; }
+      }
+    }
+    return old;
+  }
 
   // Get the underlying value without atomics
   T& underlying() { return val; }
@@ -216,36 +292,8 @@ class omp_atomic_t {
 
  private:
   T val;
-
-  friend double fetch_min(omp_atomic_t<double>& atomic_var, double other);
-  friend double fetch_max(omp_atomic_t<double>& atomic_var, double other);
 };
 
-// Free non-template functions are necessary because of a clang 20 bug
-// when omp atomic compare is used within a templated context.
-// see https://github.com/llvm/llvm-project/issues/127466
-inline double fetch_min(omp_atomic_t<double>& atomic_var, double other)
-{
-  double old;
-#pragma omp atomic compare capture
-  {
-    old = atomic_var.val;
-    if (other < atomic_var.val) { atomic_var.val = other; }
-  }
-  return old;
-}
-
-inline double fetch_max(omp_atomic_t<double>& atomic_var, double other)
-{
-  double old;
-#pragma omp atomic compare capture
-  {
-    old = atomic_var.val;
-    if (other > atomic_var.val) { atomic_var.val = other; }
-  }
-  return old;
-}
-
 }  // namespace cuopt
 
 #endif