Comparisons against EDBO+ were implemented according to instructions in the code accompanying its publication, using the same scripts but replacing the Pd-catalysed C-H arylation data set in the EDBO+ directory with the expanded C-H arylation virtual benchmark dataset generated in this work. To ensure fair comparison, EDBO+ was tested with both CPU and GPU (see Supplementary Information section 2). For noisy benchmarks in section 2.4, we perturbed objective values only for yield and catalyst turnover number for all datasets, excluding noiseless reaction costs in the C-H arylation virtual dataset. The noisy values were clamped at 0% and 100% for yield, and for turnover in the Suzuki Coupling virtual datasets, at 0 and the max observed values. The nested acquisition function described in section 2.5 first selected promising temperatures by ranking the temperatures with the highest average q-NEHVI acquisition function value, followed by standard acquisition function evaluation on the search space restricted to the obtained top performing temperatures. For the benchmarks constrained to 1 and 2 unique temperatures per batch, model initialisation data was constructed by randomly selecting 1, or 2 unique temperatures according to constraints. Then, Sobol samples were drawn from experimental conditions restricted to those temperatures to initialise optimisation process. All Wilcoxon statistical significance tests included in the Supplementary Information were implemented using SciPy .