diff --git a/docs/source/how_to/how_to_multistart.ipynb b/docs/source/how_to/how_to_multistart.ipynb
index 25d44641c..b20d4183b 100644
--- a/docs/source/how_to/how_to_multistart.ipynb
+++ b/docs/source/how_to/how_to_multistart.ipynb
@@ -7,11 +7,23 @@
    "source": [
     "# How to do multistart optimizations\n",
     "\n",
-    "## How to turn on multistart\n",
+    "Sometimes you want to make sure that your optimization is robust to the initial\n",
+    "parameter values, i.e. that it does not get stuck at a local optimum. This is where\n",
+    "multistart comes in handy.\n",
     "\n",
-    "Turning on multistart literally just means adding `multistart=True` when you call `maximize` or `minimize` and adding sampling bounds to the params DataFrame. Of course, you can configure every aspect of the multistart optimization if you want. But if you don't, we pick good defaults for you. \n",
     "\n",
-    "Let's look at the well known \"sphere\" example again:"
+    "## What does multistart (not) do\n",
+    "\n",
+    "In short, multistart iteratively runs local optimizations from different initial\n",
+    "conditions. If enough local optimization convergence to the same point, it stops.\n",
+    "Importantly, it cannot guarantee that the result is the global optimum, but it can\n",
+    "increase your confidence in the result.\n",
+    "\n",
+    "## TL;DR\n",
+    "\n",
+    "To activate multistart at the default options, pass `multistart=True` to the `minimize`\n",
+    "or `maximize` function, as well as finite bounds on the parameters (which are used to\n",
+    "sample the initial points). The default options are discussed below."
    ]
   },
   {
@@ -21,209 +33,445 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "import numpy as np\n",
     "import optimagic as om\n",
-    "import pandas as pd"
+    "\n",
+    "\n",
+    "def fun(x):\n",
+    "    return x @ x\n",
+    "\n",
+    "\n",
+    "x0 = np.arange(7) - 4\n",
+    "\n",
+    "bounds = om.Bounds(\n",
+    "    lower=np.full_like(x0, -5),\n",
+    "    upper=np.full_like(x0, 10),\n",
+    ")\n",
+    "\n",
+    "algo_options = {\"stopping_maxfun\": 1_000}\n",
+    "\n",
+    "res = om.minimize(\n",
+    "    fun=fun,\n",
+    "    x0=x0,\n",
+    "    algorithm=\"scipy_neldermead\",\n",
+    "    algo_options=algo_options,\n",
+    "    bounds=bounds,\n",
+    "    multistart=True,\n",
+    ")"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "id": "2",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "def sphere(params):\n",
-    "    return params[\"value\"] @ params[\"value\"]"
+    "In this example, we limited each local optimization to 1_000 function evaluations. In\n",
+    "general, it is a good idea to limit the number of iterations and function evaluations\n",
+    "for the local optimization. Because of the iterative nature of multistart, this\n",
+    "limitation will usually not result in a precision issue."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3",
+   "metadata": {},
+   "source": [
+    "## What does multistart mean in optimagic?\n",
+    "\n",
+    "Our multistart optimizations are inspired by the [TikTak algorithm](https://github.com/serdarozkan/TikTak) and consist of the following steps:\n",
+    "\n",
+    "1. Draw a large exploration sample of parameter vectors randomly or using a\n",
+    "   low-discrepancy sequence.\n",
+    "1. Evaluate the objective function in parallel on the exploration sample.\n",
+    "1. Sort the parameter vectors from best to worst according to their objective function\n",
+    "   values. \n",
+    "1. Run local optimizations iteratively. That is, the first local optimization is started\n",
+    "   from the best parameter vector in the sample. All subsequent ones are started from a\n",
+    "   convex combination of the currently best known parameter vector and the next sample\n",
+    "   point. "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4",
+   "metadata": {},
+   "source": [
+    "## Visualizing multistart results\n",
+    "\n",
+    "To illustrate the multistart results, we will consider the optimization of a slightly\n",
+    "more complex objective function, compared to `fun` from above. We also limit the\n",
+    "number of exploration samples to 100."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "3",
+   "id": "5",
    "metadata": {},
    "outputs": [],
    "source": [
-    "params = pd.DataFrame()\n",
-    "params[\"value\"] = [1, 2, 3]\n",
-    "params[\"soft_lower_bound\"] = -5\n",
-    "params[\"soft_upper_bound\"] = 10\n",
-    "params"
+    "def alpine(x):\n",
+    "    return np.sum(np.abs(x * np.sin(x) + 0.1 * x))\n",
+    "\n",
+    "\n",
+    "res = om.minimize(\n",
+    "    alpine,\n",
+    "    x0=x0,\n",
+    "    algorithm=\"scipy_neldermead\",\n",
+    "    bounds=bounds,\n",
+    "    algo_options=algo_options,\n",
+    "    multistart=om.MultistartOptions(n_samples=100, seed=0),\n",
+    ")\n",
+    "\n",
+    "fig = om.criterion_plot(res, monotone=True)\n",
+    "fig.show(\"png\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6",
+   "metadata": {},
+   "source": [
+    "In the above image we see the optimization history for all of the local optimizations\n",
+    "that have been run by multistart. The turquoise line represents the history\n",
+    "corresponding to the local optimization that found the overall best parameter.\n",
+    "\n",
+    "We see that running a single optimization would not have sufficed, as some local\n",
+    "optimizations are stuck."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7",
+   "metadata": {},
+   "source": [
+    "## Multistart does not always run many optimization\n",
+    "\n",
+    "Since the local optimizations are run iteratively by multistart, it is possible that\n",
+    "only a handful of optimizations are actually run if all of them converge to the same\n",
+    "point. This convergence is determined by the `convergence_max_discoveries` option,\n",
+    "which defaults to 2. This means that if 2 local optimizations report the same point,\n",
+    "multistart will stop. Below we see that if we use the simpler objective function\n",
+    "(`fun`), and the `scipy_lbfgsb` algorithm, multistart runs only 2 local optimizations,\n",
+    "and then stops, as both of them converge to the same point. Note that, the\n",
+    "`scipy_lbfgsb` algorithm can solve this simple problem precisely, without reaching the\n",
+    "maximum number of function evaluations."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "4",
+   "id": "8",
    "metadata": {},
    "outputs": [],
    "source": [
     "res = om.minimize(\n",
-    "    fun=sphere,\n",
-    "    params=params,\n",
+    "    fun,\n",
+    "    x0=x0,\n",
     "    algorithm=\"scipy_lbfgsb\",\n",
-    "    multistart=True,\n",
+    "    bounds=bounds,\n",
+    "    algo_options=algo_options,\n",
+    "    multistart=om.MultistartOptions(n_samples=100, seed=0),\n",
     ")\n",
-    "res.params.round(5)"
+    "\n",
+    "fig = om.criterion_plot(res)\n",
+    "fig.show(\"png\")"
    ]
   },
   {
    "cell_type": "markdown",
-   "id": "5",
+   "id": "9",
    "metadata": {},
    "source": [
-    "## Understanding multistart results\n",
+    "## How to configure multistart\n",
+    "\n",
+    "Configuration of multistart can be done by passing an instance of\n",
+    "`optimagic.MultistartOptions` to `minimize` or `maximize`. Let's look at a few examples\n",
+    "configurations."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "10",
+   "metadata": {},
+   "source": [
+    "### How to run a specific number of optimizations\n",
+    "\n",
+    "To run a specific number of local optimizations, you need to set the `stopping_maxopt`\n",
+    "option. Note that this does not set the number of exploration samples, which is\n",
+    "controlled by the `n_samples` option. The number of exploration samples always needs\n",
+    "to be at least as large as the number of local optimizations.\n",
+    "\n",
+    "Note that, as long as `convergence_max_discoveries` is smaller than `stopping_maxopt`,\n",
+    "it is possible that a smaller number of local optimizations are run. To avoid this,\n",
+    "set `convergence_max_discoveries` to a value at least as large as `stopping_maxopt`.\n",
     "\n",
-    "The ``OptimizeResult`` result object of a multistart optimization has exactly the same structure as ``OptimizeResult`` from a standard optimization but with additional information.\n",
+    "To run, for example, 10 local optimizations from 15 exploration samples, do:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "11",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "res = om.minimize(\n",
+    "    alpine,\n",
+    "    x0=x0,\n",
+    "    algorithm=\"scipy_neldermead\",\n",
+    "    bounds=bounds,\n",
+    "    algo_options=algo_options,\n",
+    "    multistart=om.MultistartOptions(\n",
+    "        n_samples=15,\n",
+    "        stopping_maxopt=10,\n",
+    "        convergence_max_discoveries=10,\n",
+    "    ),\n",
+    ")\n",
     "\n",
-    "- `res.multistart_info[\"local_optima\"]` is a list with the results from all optimizations that were performed\n",
-    "- `res.multistart_info[\"start_parameters\"]` is a list with the start parameters from those optimizations \n",
-    "- `res.multistart_info[\"exploration_sample\"]` is a list with parameter vectors at which the criterion function was evaluated in an initial exploration phase. \n",
-    "- `res.multistart_info[\"exploration_results\"]` are the corresponding criterion values. "
+    "res.multistart_info.n_optimizations"
    ]
   },
   {
    "cell_type": "markdown",
-   "id": "6",
+   "id": "12",
    "metadata": {},
    "source": [
-    "## What does multistart mean in optimagic?\n",
+    "### How to set a custom exploration sample\n",
+    "\n",
+    "If you want to start the multistart algorithm with a custom exploration sample, you can\n",
+    "do so by passing a sequence of parameters to the `sample` option. Note that sequence\n",
+    "elements must be of the same type as your parameter.\n",
+    "\n",
+    "To generate a sample of 100 random parameters and run them through the multistart\n",
+    "algorithm, do:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "13",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rng = np.random.default_rng(12345)\n",
     "\n",
-    "The way we do multistart optimizations is inspired by the [TikTak algorithm](https://github.com/serdarozkan/TikTak). Our multistart optimizations consist of the following steps:\n",
+    "sample = [x0 + rng.uniform(-1, 1, size=len(x0)) for _ in range(100)]\n",
     "\n",
-    "1. Draw a large sample of parameter vectors, randomly or using a low-discrepancy sequence. The size, sampling method, distribution, and more things can be configured. \n",
-    "2. Evaluate the criterion function in parallel on all parameter vectors.\n",
-    "4. Sort the parameter vectors from best to worst. \n",
-    "5. Run local optimizations. The first local optimization is started from the best parameter vector in the sample. All subsequent ones are started from a convex combination of the currently best known parameter vector and the next sample point. "
+    "res = om.minimize(\n",
+    "    alpine,\n",
+    "    x0=x0,\n",
+    "    algorithm=\"scipy_neldermead\",\n",
+    "    bounds=bounds,\n",
+    "    algo_options=algo_options,\n",
+    "    multistart=om.MultistartOptions(sample=sample),\n",
+    ")"
    ]
   },
   {
    "cell_type": "markdown",
-   "id": "7",
+   "id": "14",
    "metadata": {},
    "source": [
-    "## How to configure mutlistart?\n",
+    "### How to run multistart in parallel\n",
+    "\n",
+    "\n",
+    "The multistart algorithm can be run in parallel by setting the `n_cores` option to a\n",
+    "value greater than 1. This will run the algorithm in batches. By default, the batch\n",
+    "size is set to `n_cores`, but can be controlled by setting the `batch_size` option. The\n",
+    "default batch evaluator is `joblib`, but can be controlled by setting the\n",
+    "`batch_evaluator` option to `\"pathos\"` or a custom callable.\n",
     "\n",
-    "As you can imagine from the above description, there are many details that can be configured. This can be done by adding a dictionary with `multistart_options` when calling `minimize` or `maximize`. Let's look at an extreme example where we manually set everything to it's default value:"
+    "To run the multistart algorithm in parallel, do:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "8",
+   "id": "15",
    "metadata": {},
    "outputs": [],
    "source": [
-    "options = {\n",
-    "    # Set the number of points at which criterion is evaluated\n",
-    "    # in the exploration phase\n",
-    "    \"n_samples\": 10 * len(params),\n",
-    "    # Pass in a DataFrame or array with a custom sample\n",
-    "    # for the exploration phase.\n",
-    "    \"sample\": None,\n",
-    "    # Determine number of optimizations, relative to n_samples\n",
-    "    \"share_optimizations\": 0.1,\n",
-    "    # Determine distribution from which sample is drawn\n",
-    "    \"sampling_distribution\": \"uniform\",\n",
-    "    # Determine sampling method. Allowed: [\"sobol\", \"random\",\n",
-    "    # \"halton\", \"hammersley\", \"korobov\", \"latin_hypercube\"]\n",
-    "    \"sampling_method\": \"sobol\",\n",
-    "    # Determine how start parameters for local optimizations are\n",
-    "    # calculated. Allowed: [\"tiktak\", \"linear\"] or a custom\n",
-    "    # function with arguments iteration, n_iterations, min_weight,\n",
-    "    # and max_weight\n",
-    "    \"mixing_weight_method\": \"tiktak\",\n",
-    "    # Determine bounds on mixing weights.\n",
-    "    \"mixing_weight_bounds\": (0.1, 0.995),\n",
-    "    # Determine after how many re-discoveries of the currently best\n",
-    "    # local optimum the multistart optimization converges.\n",
-    "    \"convergence.max_discoveries\": 2,\n",
-    "    # Determine the maximum relative distance two parameter vectors\n",
-    "    # can have to be considered equal for convergence purposes:\n",
-    "    \"convergence.relative_params_tolerance\": 0.01,\n",
-    "    # Determine how many cores are used\n",
-    "    \"n_cores\": 1,\n",
-    "    # Determine which batch_evaluator is used:\n",
-    "    \"batch_evaluator\": \"joblib\",\n",
-    "    # Determine the batch size. It must be larger than n_cores.\n",
-    "    # Setting the batch size larger than n_cores allows to reproduce\n",
-    "    # the exact results of a highly parallel optimization on a smaller\n",
-    "    # machine.\n",
-    "    \"batch_size\": 1,\n",
-    "    # Set the random seed:\n",
-    "    \"seed\": None,\n",
-    "    # Set how errors are handled during the exploration phase:\n",
-    "    \"exploration_error_handling\": \"continue\",\n",
-    "    # Set how errors are handled during the optimization phase:\n",
-    "    \"optimization_error_handling\": \"continue\",\n",
-    "}"
+    "res = om.minimize(\n",
+    "    alpine,\n",
+    "    x0=x0,\n",
+    "    algorithm=\"scipy_lbfgsb\",\n",
+    "    bounds=bounds,\n",
+    "    algo_options=algo_options,\n",
+    "    multistart=om.MultistartOptions(n_cores=2),\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "16",
+   "metadata": {},
+   "source": [
+    "## What to do if you do not have bounds\n",
+    "\n",
+    "Multistart requires finite bounds on the parameters. If your optimization problem is not\n",
+    "bounded, you can set soft lower and upper bounds. These bounds will only be used to\n",
+    "draw the exploration sample, and will not be used to constrain the local optimizations.\n",
+    "\n",
+    "To set soft bounds, do:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "9",
+   "id": "17",
    "metadata": {},
    "outputs": [],
    "source": [
     "res = om.minimize(\n",
-    "    fun=sphere,\n",
-    "    params=params,\n",
+    "    alpine,\n",
+    "    x0=x0,\n",
     "    algorithm=\"scipy_lbfgsb\",\n",
+    "    bounds=om.Bounds(soft_lower=np.full_like(x0, -3), soft_upper=np.full_like(x0, 8)),\n",
     "    multistart=True,\n",
-    "    multistart_options=options,\n",
-    ")\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "18",
+   "metadata": {},
+   "source": [
+    "## Understanding multistart results\n",
+    "\n",
+    "When activating multistart, the optimization result object corresponds to the local\n",
+    "optimization that found the best objective function value. The result object has the\n",
+    "additional attribute `multistart_info`, where all of the additional information is\n",
+    "stored. It has the following attributes:\n",
+    "\n",
+    "- `local_optima`: A list with the results from all local optimizations that were performed.\n",
+    "- `start_parameters`: A list with the start parameters from those optimizations \n",
+    "- `exploration_sample`: A list with parameter vectors at which the objective function was evaluated in an initial exploration phase. \n",
+    "- `exploration_results`: The corresponding objective values.\n",
+    "- `n_optimizations`: The number of local optimizations that were run.\n",
     "\n",
-    "res"
+    "To illustrate the multistart results, let us consider the optimization of the simple\n",
+    "`fun` objective function from above."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "10",
+   "id": "19",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "res = om.minimize(\n",
+    "    fun,\n",
+    "    x0=x0,\n",
+    "    algorithm=\"scipy_lbfgsb\",\n",
+    "    bounds=bounds,\n",
+    "    algo_options=algo_options,\n",
+    "    multistart=om.MultistartOptions(n_samples=100, convergence_max_discoveries=2),\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "20",
+   "metadata": {},
+   "source": [
+    "### Start parameters\n",
+    "\n",
+    "The start parameters are the parameter vectors from which the local optimizations were\n",
+    "started. Since the default number of `convergence_max_discoveries` is 2, and both\n",
+    "local optimizations were successfull, the start parameters have 2 rows."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "21",
    "metadata": {},
    "outputs": [],
    "source": [
-    "res.params.round(5)"
+    "res.multistart_info.start_parameters"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "22",
+   "metadata": {},
+   "source": [
+    "### Local Optima\n",
+    "\n",
+    "The local optima are the results from the local optimizations. Since in this example\n",
+    "only two local optimizations were run, the local optima list has two elements, each of\n",
+    "which is an optimization result object."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "11",
+   "id": "23",
    "metadata": {},
    "outputs": [],
    "source": [
-    "res.multistart_info[\"local_optima\"]"
+    "len(res.multistart_info.local_optima)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "24",
+   "metadata": {},
+   "source": [
+    "### Exploration sample\n",
+    "\n",
+    "The exploration sample is a list of parameter vectors at which the objective function\n",
+    "was evaluated. Above, we chose a random exploration sample of 100 parameter vectors."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "12",
+   "id": "25",
    "metadata": {},
    "outputs": [],
    "source": [
-    "res.multistart_info[\"start_parameters\"]"
+    "np.row_stack(res.multistart_info.exploration_sample).shape"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "26",
+   "metadata": {},
+   "source": [
+    "### Exploration results\n",
+    "\n",
+    "The exploration results are the objective function values at the exploration sample."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "13",
+   "id": "27",
    "metadata": {},
    "outputs": [],
    "source": [
-    "res.multistart_info[\"start_parameters\"]"
+    "len(res.multistart_info.exploration_results)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "28",
+   "metadata": {},
+   "source": [
+    "### Number of local optimizations"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "14",
+   "id": "29",
    "metadata": {},
    "outputs": [],
    "source": [
-    "res.multistart_info[\"exploration_results\"]"
+    "res.multistart_info.n_optimizations"
    ]
   }
  ],
diff --git a/docs/source/how_to/how_to_visualize_histories.ipynb b/docs/source/how_to/how_to_visualize_histories.ipynb
index 640be7a4e..37b3e43a1 100644
--- a/docs/source/how_to/how_to_visualize_histories.ipynb
+++ b/docs/source/how_to/how_to_visualize_histories.ipynb
@@ -177,12 +177,11 @@
     "\n",
     "\n",
     "res = om.minimize(\n",
-    "    sphere,\n",
-    "    params=np.arange(10),\n",
-    "    bounds=om.Bounds(soft_lower=np.full(10, -3), soft_upper=np.full(10, 10)),\n",
+    "    alpine,\n",
+    "    params=np.arange(7),\n",
+    "    bounds=om.Bounds(soft_lower=np.full(7, -3), soft_upper=np.full(7, 10)),\n",
     "    algorithm=\"scipy_neldermead\",\n",
-    "    multistart=True,\n",
-    "    multistart_options={\"n_samples\": 1000, \"convergence.max_discoveries\": 10},\n",
+    "    multistart=om.MultistartOptions(n_samples=100, convergence_max_discoveries=3),\n",
     ")"
    ]
   },
@@ -193,7 +192,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fig = om.criterion_plot(res, max_evaluations=3000)\n",
+    "fig = om.criterion_plot(res, max_evaluations=1000, monotone=True)\n",
     "fig.show(renderer=\"png\")"
    ]
   }
diff --git a/docs/source/tutorials/optimization_overview.ipynb b/docs/source/tutorials/optimization_overview.ipynb
index de250c213..bcdb06927 100644
--- a/docs/source/tutorials/optimization_overview.ipynb
+++ b/docs/source/tutorials/optimization_overview.ipynb
@@ -355,8 +355,7 @@
     "    params=np.arange(10),\n",
     "    algorithm=\"scipy_neldermead\",\n",
     "    bounds=bounds,\n",
-    "    multistart=True,\n",
-    "    multistart_options={\"convergence.max_discoveries\": 5},\n",
+    "    multistart=om.MultistartOptions(convergence_max_discoveries=5),\n",
     ")\n",
     "res.params.round(5)"
    ]
diff --git a/pyproject.toml b/pyproject.toml
index fc183ecea..9a60f537b 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -248,7 +248,6 @@ module = [
     "optimagic.optimization.optimization_logging",
     "optimagic.optimization.optimize_result",
     "optimagic.optimization.optimize",
-    "optimagic.optimization.process_multistart_sample",
     "optimagic.optimization.process_results",
     "optimagic.optimization.multistart",
     "optimagic.optimization.scipy_aliases",
diff --git a/src/optimagic/__init__.py b/src/optimagic/__init__.py
index 52a485cb7..11bb09b49 100644
--- a/src/optimagic/__init__.py
+++ b/src/optimagic/__init__.py
@@ -8,6 +8,7 @@
 from optimagic.benchmarking.run_benchmark import run_benchmark
 from optimagic.differentiation.derivatives import first_derivative, second_derivative
 from optimagic.logging.read_log import OptimizeLogReader
+from optimagic.optimization.multistart_options import MultistartOptions
 from optimagic.optimization.optimize import maximize, minimize
 from optimagic.optimization.optimize_result import OptimizeResult
 from optimagic.parameters.bounds import Bounds
@@ -48,5 +49,6 @@
     "OptimizeResult",
     "Bounds",
     "ScalingOptions",
+    "MultistartOptions",
     "__version__",
 ]
diff --git a/src/optimagic/deprecations.py b/src/optimagic/deprecations.py
index 3471a5a35..b746cb951 100644
--- a/src/optimagic/deprecations.py
+++ b/src/optimagic/deprecations.py
@@ -1,4 +1,5 @@
 import warnings
+from dataclasses import replace
 
 from optimagic.parameters.bounds import Bounds
 
@@ -71,6 +72,15 @@ def throw_scaling_options_future_warning():
     warnings.warn(msg, FutureWarning)
 
 
+def throw_multistart_options_future_warning():
+    msg = (
+        "Specifying multistart options via the argument `multistart_options` is "
+        "deprecated and will be removed in optimagic version 0.6.0 and later. You can "
+        "pass these options directly to the `multistart` argument instead."
+    )
+    warnings.warn(msg, FutureWarning)
+
+
 def replace_and_warn_about_deprecated_algo_options(algo_options):
     if not isinstance(algo_options, dict):
         return algo_options
@@ -138,3 +148,59 @@ def replace_and_warn_about_deprecated_bounds(
         bounds = Bounds(**old_bounds)
 
     return bounds
+
+
+def replace_and_warn_about_deprecated_multistart_options(options):
+    """Replace deprecated multistart options and warn about them.
+
+    Args:
+        options (MultistartOptions): The multistart options to replace.
+
+    Returns:
+        MultistartOptions: The replaced multistart options.
+
+    """
+    replacements = {}
+
+    if options.share_optimization is not None:
+        msg = (
+            "The share_optimization option is deprecated and will be removed in "
+            "version 0.6.0. Use stopping_maxopt instead to specify the number of "
+            "optimizations directly."
+        )
+        warnings.warn(msg, FutureWarning)
+
+    if options.convergence_relative_params_tolerance is not None:
+        msg = (
+            "The convergence_relative_params_tolerance option is deprecated and will "
+            "be removed in version 0.6.0. Use convergence_xtol_rel instead."
+        )
+        warnings.warn(msg, FutureWarning)
+        if options.convergence_xtol_rel is None:
+            replacements["convergence_xtol_rel"] = (
+                options.convergence_relative_params_tolerance
+            )
+
+    if options.optimization_error_handling is not None:
+        msg = (
+            "The optimization_error_handling option is deprecated and will be removed "
+            "in version 0.6.0. Setting this attribute also sets the error handling "
+            "for exploration. Use the new error_handling option to set the error "
+            "handling for both optimization and exploration."
+        )
+        warnings.warn(msg, FutureWarning)
+        if options.error_handling is None:
+            replacements["error_handling"] = options.optimization_error_handling
+
+    if options.exploration_error_handling is not None:
+        msg = (
+            "The exploration_error_handling option is deprecated and will be "
+            "removed in version 0.6.0. Setting this attribute also sets the error "
+            "handling for exploration. Use the new error_handling option to set the "
+            "error handling for both optimization and exploration."
+        )
+        warnings.warn(msg, FutureWarning)
+        if options.error_handling is None:
+            replacements["error_handling"] = options.exploration_error_handling
+
+    return replace(options, **replacements)
diff --git a/src/optimagic/exceptions.py b/src/optimagic/exceptions.py
index b41080d52..f4c5219b1 100644
--- a/src/optimagic/exceptions.py
+++ b/src/optimagic/exceptions.py
@@ -51,6 +51,10 @@ class InvalidScalingError(OptimagicError):
     """Exception for invalid user provided scaling."""
 
 
+class InvalidMultistartError(OptimagicError):
+    """Exception for invalid user provided multistart options."""
+
+
 class NotInstalledError(OptimagicError):
     """Exception when optional dependencies are needed but not installed."""
 
diff --git a/src/optimagic/optimization/create_optimization_problem.py b/src/optimagic/optimization/create_optimization_problem.py
index 8710c09aa..5cb27f5b1 100644
--- a/src/optimagic/optimization/create_optimization_problem.py
+++ b/src/optimagic/optimization/create_optimization_problem.py
@@ -15,6 +15,10 @@
 from optimagic.optimization.get_algorithm import (
     process_user_algorithm,
 )
+from optimagic.optimization.multistart_options import (
+    MultistartOptions,
+    pre_process_multistart,
+)
 from optimagic.optimization.scipy_aliases import (
     map_method_to_algorithm,
     split_fun_and_jac,
@@ -70,10 +74,7 @@ class OptimizationProblem:
     error_handling: Literal["raise", "continue"]
     error_penalty: dict[str, Any] | None
     scaling: ScalingOptions | None
-    # TODO: multistart will become None | MultistartOptions and multistart_options will
-    # be removed
-    multistart: bool
-    multistart_options: dict[str, Any] | None
+    multistart: MultistartOptions | None
     collect_history: bool
     skip_checks: bool
     direction: Literal["minimize", "maximize"]
@@ -100,7 +101,6 @@ def create_optimization_problem(
     error_penalty,
     scaling,
     multistart,
-    multistart_options,
     collect_history,
     skip_checks,
     # scipy aliases
@@ -126,6 +126,7 @@ def create_optimization_problem(
     soft_lower_bounds,
     soft_upper_bounds,
     scaling_options,
+    multistart_options,
 ):
     # ==================================================================================
     # error handling needed as long as fun is an optional argument (i.e. until
@@ -187,7 +188,13 @@ def create_optimization_problem(
 
     if scaling_options is not None:
         deprecations.throw_scaling_options_future_warning()
-        scaling = scaling_options if scaling is None else scaling
+        if scaling is True and scaling_options is not None:
+            scaling = scaling_options
+
+    if multistart_options is not None:
+        deprecations.throw_multistart_options_future_warning()
+        if multistart is True and multistart_options is not None:
+            multistart = multistart_options
 
     algo_options = replace_and_warn_about_deprecated_algo_options(algo_options)
 
@@ -303,6 +310,7 @@ def create_optimization_problem(
     # ==================================================================================
     bounds = pre_process_bounds(bounds)
     scaling = pre_process_scaling(scaling)
+    multistart = pre_process_multistart(multistart)
 
     fun_kwargs = {} if fun_kwargs is None else fun_kwargs
     constraints = [] if constraints is None else constraints
@@ -312,7 +320,6 @@ def create_optimization_problem(
     numdiff_options = {} if numdiff_options is None else numdiff_options
     log_options = {} if log_options is None else log_options
     error_penalty = {} if error_penalty is None else error_penalty
-    multistart_options = {} if multistart_options is None else multistart_options
     if logging:
         logging = Path(logging)
 
@@ -405,11 +412,8 @@ def create_optimization_problem(
         if not isinstance(scaling, ScalingOptions | None):
             raise ValueError("scaling must be a ScalingOptions object or None")
 
-        if not isinstance(multistart, bool):
-            raise ValueError("multistart must be a boolean")
-
-        if not isinstance(multistart_options, dict | None):
-            raise ValueError("multistart_options must be a dictionary or None")
+        if not isinstance(multistart, MultistartOptions | None):
+            raise ValueError("multistart must be a MultistartOptions object or None")
 
         if not isinstance(collect_history, bool):
             raise ValueError("collect_history must be a boolean")
@@ -446,7 +450,6 @@ def create_optimization_problem(
         error_penalty=error_penalty,
         scaling=scaling,
         multistart=multistart,
-        multistart_options=multistart_options,
         collect_history=collect_history,
         skip_checks=skip_checks,
         direction=direction,
diff --git a/src/optimagic/optimization/multistart.py b/src/optimagic/optimization/multistart.py
index 344d8d519..b141d5360 100644
--- a/src/optimagic/optimization/multistart.py
+++ b/src/optimagic/optimization/multistart.py
@@ -13,8 +13,10 @@
 
 import warnings
 from functools import partial
+from typing import Literal
 
 import numpy as np
+from numpy.typing import NDArray
 from scipy.stats import qmc, triang
 
 from optimagic.batch_evaluators import process_batch_evaluator
@@ -39,7 +41,7 @@ def run_multistart_optimization(
     database,
     error_handling,
 ):
-    steps = determine_steps(options["n_samples"], options["n_optimizations"])
+    steps = determine_steps(options.n_samples, stopping_maxopt=options.stopping_maxopt)
 
     scheduled_steps = log_scheduled_steps_and_get_ids(
         steps=steps,
@@ -47,18 +49,18 @@ def run_multistart_optimization(
         database=database,
     )
 
-    if options["sample"] is not None:
-        sample = options["sample"]
+    if options.sample is not None:
+        sample = options.sample
     else:
-        sample = draw_exploration_sample(
+        sample = _draw_exploration_sample(
             x=x,
             lower=lower_sampling_bounds,
             upper=upper_sampling_bounds,
             # -1 because we add start parameters
-            n_samples=options["n_samples"] - 1,
-            sampling_distribution=options["sampling_distribution"],
-            sampling_method=options["sampling_method"],
-            seed=options["seed"],
+            n_samples=options.n_samples - 1,
+            distribution=options.sampling_distribution,
+            method=options.sampling_method,
+            seed=options.seed,
         )
 
         sample = np.vstack([x.reshape(1, -1), sample])
@@ -79,10 +81,10 @@ def run_multistart_optimization(
         criterion,
         primary_key=primary_key,
         sample=sample,
-        batch_evaluator=options["batch_evaluator"],
-        n_cores=options["n_cores"],
+        batch_evaluator=options.batch_evaluator,
+        n_cores=options.n_cores,
         step_id=scheduled_steps[0],
-        error_handling=options["exploration_error_handling"],
+        error_handling=options.error_handling,
     )
 
     if logging:
@@ -97,14 +99,14 @@ def run_multistart_optimization(
     sorted_sample = exploration_res["sorted_sample"]
     sorted_values = exploration_res["sorted_values"]
 
-    n_optimizations = options["n_optimizations"]
-    if n_optimizations > len(sorted_sample):
-        n_skipped_steps = n_optimizations - len(sorted_sample)
-        n_optimizations = len(sorted_sample)
+    stopping_maxopt = options.stopping_maxopt
+    if stopping_maxopt > len(sorted_sample):
+        n_skipped_steps = stopping_maxopt - len(sorted_sample)
+        stopping_maxopt = len(sorted_sample)
         warnings.warn(
             "There are less valid starting points than requested optimizations. "
             "The number of optimizations has been reduced from "
-            f"{options['n_optimizations']} to {len(sorted_sample)}."
+            f"{options.stopping_maxopt} to {len(sorted_sample)}."
         )
         skipped_steps = scheduled_steps[-n_skipped_steps:]
         scheduled_steps = scheduled_steps[:-n_skipped_steps]
@@ -119,8 +121,8 @@ def run_multistart_optimization(
 
     batched_sample = get_batched_optimization_sample(
         sorted_sample=sorted_sample,
-        n_optimizations=n_optimizations,
-        batch_size=options["batch_size"],
+        stopping_maxopt=stopping_maxopt,
+        batch_size=options.batch_size,
     )
 
     state = {
@@ -134,8 +136,8 @@ def run_multistart_optimization(
     }
 
     convergence_criteria = {
-        "xtol": options["convergence_relative_params_tolerance"],
-        "max_discoveries": options["convergence_max_discoveries"],
+        "xtol": options.convergence_xtol_rel,
+        "max_discoveries": options.convergence_max_discoveries,
     }
 
     problem_functions = {
@@ -143,17 +145,11 @@ def run_multistart_optimization(
         for name, func in problem_functions.items()
     }
 
-    batch_evaluator = options["batch_evaluator"]
-
-    weight_func = partial(
-        options["mixing_weight_method"],
-        min_weight=options["mixing_weight_bounds"][0],
-        max_weight=options["mixing_weight_bounds"][1],
-    )
+    batch_evaluator = options.batch_evaluator
 
     opt_counter = 0
     for batch in batched_sample:
-        weight = weight_func(opt_counter, n_optimizations)
+        weight = options.weight_func(opt_counter, stopping_maxopt)
         starts = [weight * state["best_x"] + (1 - weight) * x for x in batch]
 
         arguments = [
@@ -165,8 +161,8 @@ def run_multistart_optimization(
             func=local_algorithm,
             arguments=arguments,
             unpack_symbol="**",
-            n_cores=options["n_cores"],
-            error_handling=options["optimization_error_handling"],
+            n_cores=options.n_cores,
+            error_handling=options.error_handling,
         )
 
         state, is_converged = update_convergence_state(
@@ -188,18 +184,20 @@ def run_multistart_optimization(
                     )
             break
 
-    raw_res = state["best_res"]
-    raw_res["multistart_info"] = {
+    multistart_info = {
         "start_parameters": state["start_history"],
         "local_optima": state["result_history"],
         "exploration_sample": sorted_sample,
         "exploration_results": exploration_res["sorted_values"],
     }
 
+    raw_res = state["best_res"]
+    raw_res["multistart_info"] = multistart_info
+
     return raw_res
 
 
-def determine_steps(n_samples, n_optimizations):
+def determine_steps(n_samples, stopping_maxopt):
     """Determine the number and type of steps for the multistart optimization.
 
     This is mainly used to write them to the log. The number of steps is also
@@ -207,7 +205,7 @@ def determine_steps(n_samples, n_optimizations):
 
     Args:
         n_samples (int): Number of exploration points for the multistart optimization.
-        n_optimizations (int): Number of local optimizations.
+        stopping_maxopt (int): Number of local optimizations.
 
 
     Returns:
@@ -222,7 +220,7 @@ def determine_steps(n_samples, n_optimizations):
     }
 
     steps = [exploration_step]
-    for i in range(n_optimizations):
+    for i in range(stopping_maxopt):
         optimization_step = {
             "type": "optimization",
             "status": "scheduled",
@@ -232,49 +230,36 @@ def determine_steps(n_samples, n_optimizations):
     return steps
 
 
-def draw_exploration_sample(
-    x,
-    lower,
-    upper,
-    n_samples,
-    sampling_distribution,
-    sampling_method,
-    seed,
-):
+def _draw_exploration_sample(
+    x: NDArray[np.float64],
+    lower: NDArray[np.float64],
+    upper: NDArray[np.float64],
+    n_samples: int,
+    distribution: Literal["uniform", "triangular"],
+    method: Literal["sobol", "random", "halton", "latin_hypercube"],
+    seed: int | np.random.Generator | None,
+) -> NDArray[np.float64]:
     """Get a sample of parameter values for the first stage of the tiktak algorithm.
 
     The sample is created randomly or using a low discrepancy sequence. Different
     distributions are available.
 
     Args:
-        x (np.ndarray): Internal parameter vector of shape (n_params,).
-        lower (np.ndarray): Vector of internal lower bounds of shape (n_params,).
-        upper (np.ndarray): Vector of internal upper bounds of shape (n_params,).
-        n_samples (int): Number of sample points on which one function evaluation
-            shall be performed. Default is 10 * n_params.
-        sampling_distribution (str): One of "uniform", "triangular". Default is
-            "uniform", as in the original tiktak algorithm.
-        sampling_method (str): One of "sobol", "halton", "latin_hypercube" or
-            "random". Default is sobol for problems with up to 200 parameters
-            and random for problems with more than 200 parameters.
-        seed (int): Random seed.
+        x: Internal parameter vector of shape (n_params,).
+        lower: Vector of internal lower bounds of shape (n_params,).
+        upper: Vector of internal upper bounds of shape (n_params,).
+        n_samples: Number of sample points.
+        distribution: The distribution from which the exploration sample is
+            drawn. Allowed are "uniform" and "triangular". Defaults to "uniform".
+        method: The method used to draw the exploration sample. Allowed are
+            "sobol", "random", "halton", and "latin_hypercube". Defaults to "sobol".
+        seed: Random number seed or generator.
 
     Returns:
-        np.ndarray: Numpy array of shape (n_samples, n_params).
-            Each row represents a vector of parameter values.
+        Array of shape (n_samples, n_params). Each row represents a vector of parameter
+            values.
 
     """
-    valid_rules = ["sobol", "halton", "latin_hypercube", "random"]
-    valid_distributions = ["uniform", "triangular"]
-
-    if sampling_method not in valid_rules:
-        raise ValueError(
-            f"Invalid rule: {sampling_method}. Must be one of\n\n{valid_rules}\n\n"
-        )
-
-    if sampling_distribution not in valid_distributions:
-        raise ValueError(f"Unsupported distribution: {sampling_distribution}")
-
     for name, bound in zip(["lower", "upper"], [lower, upper], strict=False):
         if not np.isfinite(bound).all():
             raise ValueError(
@@ -282,7 +267,7 @@ def draw_exploration_sample(
                 f"soft_{name}_bounds for all parameters."
             )
 
-    if sampling_method == "sobol":
+    if method == "sobol":
         # Draw `n` points from the open interval (lower, upper)^d.
         # Note that scipy uses the half-open interval [lower, upper)^d internally.
         # We apply a burn-in phase of 1, i.e. we skip the first point in the sequence
@@ -291,21 +276,21 @@ def draw_exploration_sample(
         _ = sampler.fast_forward(1)
         sample_unscaled = sampler.random(n=n_samples)
 
-    elif sampling_method == "halton":
+    elif method == "halton":
         sampler = qmc.Halton(d=len(lower), scramble=False, seed=seed)
         sample_unscaled = sampler.random(n=n_samples)
 
-    elif sampling_method == "latin_hypercube":
+    elif method == "latin_hypercube":
         sampler = qmc.LatinHypercube(d=len(lower), strength=1, seed=seed)
         sample_unscaled = sampler.random(n=n_samples)
 
-    elif sampling_method == "random":
+    elif method == "random":
         rng = get_rng(seed)
         sample_unscaled = rng.uniform(size=(n_samples, len(lower)))
 
-    if sampling_distribution == "uniform":
+    if distribution == "uniform":
         sample_scaled = qmc.scale(sample_unscaled, lower, upper)
-    elif sampling_distribution == "triangular":
+    elif distribution == "triangular":
         sample_scaled = triang.ppf(
             sample_unscaled,
             c=(x - lower) / (upper - lower),
@@ -402,7 +387,7 @@ def run_explorations(
     return out
 
 
-def get_batched_optimization_sample(sorted_sample, n_optimizations, batch_size):
+def get_batched_optimization_sample(sorted_sample, stopping_maxopt, batch_size):
     """Create a batched sample of internal parameters for the optimization phase.
 
     Note that in the end the optimizations will not be started from those parameter
@@ -412,7 +397,7 @@ def get_batched_optimization_sample(sorted_sample, n_optimizations, batch_size):
     Args:
         sorted_sample (np.ndarray): 2d numpy array with containing sorted internal
             parameter vectors.
-        n_optimizations (int): Number of optimizations to run. If sample is shorter
+        stopping_maxopt (int): Number of optimizations to run. If sample is shorter
             than that, optimizations are run on all entries of the sample.
         batch_size (int): Batch size.
 
@@ -421,12 +406,12 @@ def get_batched_optimization_sample(sorted_sample, n_optimizations, batch_size):
             The inner lists have length ``batch_size`` or shorter.
 
     """
-    n_batches = int(np.ceil(n_optimizations / batch_size))
+    n_batches = int(np.ceil(stopping_maxopt / batch_size))
 
     start = 0
     batched = []
     for _ in range(n_batches):
-        stop = min(start + batch_size, len(sorted_sample), n_optimizations)
+        stop = min(start + batch_size, len(sorted_sample), stopping_maxopt)
         batched.append(list(sorted_sample[start:stop]))
         start = stop
     return batched
@@ -536,19 +521,3 @@ def update_convergence_state(
     }
 
     return new_state, is_converged
-
-
-def _tiktak_weights(iteration, n_iterations, min_weight, max_weight):
-    return np.clip(np.sqrt(iteration / n_iterations), min_weight, max_weight)
-
-
-def _linear_weights(iteration, n_iterations, min_weight, max_weight):
-    unscaled = iteration / n_iterations
-    span = max_weight - min_weight
-    return min_weight + unscaled * span
-
-
-WEIGHT_FUNCTIONS = {
-    "tiktak": _tiktak_weights,
-    "linear": _linear_weights,
-}
diff --git a/src/optimagic/optimization/multistart_options.py b/src/optimagic/optimization/multistart_options.py
new file mode 100644
index 000000000..644f068f9
--- /dev/null
+++ b/src/optimagic/optimization/multistart_options.py
@@ -0,0 +1,436 @@
+from dataclasses import dataclass
+from functools import partial
+from typing import Callable, Literal, Sequence, TypedDict, cast
+
+import numpy as np
+from numpy.typing import NDArray
+from typing_extensions import NotRequired
+
+from optimagic.batch_evaluators import process_batch_evaluator
+from optimagic.deprecations import replace_and_warn_about_deprecated_multistart_options
+from optimagic.exceptions import InvalidMultistartError
+from optimagic.typing import PyTree
+
+# ======================================================================================
+# Public Options
+# ======================================================================================
+
+
+@dataclass(frozen=True)
+class MultistartOptions:
+    """Multistart options in optimization problems.
+
+    Attributes:
+        n_samples: The number of points at which the objective function is evaluated
+            during the exploration phase. If None, n_samples is set to 100 times the
+            number of parameters.
+        stopping_maxopt: The maximum number of local optimizations to run. Defaults to
+            10% of n_samples. This number may not be reached if multistart converges
+            earlier.
+        sampling_distribution: The distribution from which the exploration sample is
+            drawn. Allowed are "uniform" and "triangular". Defaults to "uniform".
+        sampling_method: The method used to draw the exploration sample. Allowed are
+            "sobol", "random", "halton", and "latin_hypercube". Defaults to "random".
+        sample: A sequence of PyTrees or None. If None, a sample is drawn from the
+            sampling distribution.
+        mixing_weight_method: The method used to determine the mixing weight, i,e, how
+            start parameters for local optimizations are calculated. Allowed are
+            "tiktak" and "linear", or a custom callable. Defaults to "tiktak".
+        mixing_weight_bounds: The lower and upper bounds for the mixing weight.
+            Defaults to (0.1, 0.995).
+        convergence_max_discoveries: The maximum number of discoveries for convergence.
+            Determines after how many re-descoveries of the currently best local
+            optima the multistart algorithm stops. Defaults to 2.
+        convergence_xtol_rel: The relative tolerance in parameters
+            for convergence. Determines the maximum relative distance two parameter
+            vecctors can have to be considered equal. Defaults to 0.01.
+        n_cores: The number of cores to use for parallelization. Defaults to 1.
+        batch_evaluator: The evaluator to use for batch evaluation. Allowed are "joblib"
+            and "pathos", or a custom callable.
+        batch_size: The batch size for batch evaluation. Must be larger than n_cores
+            or None.
+        seed: The seed for the random number generator.
+        error_handling: The error handling for exploration and optimization errors.
+            Allowed are "raise" and "continue".
+
+    Raises:
+        InvalidMultistartError: If the multistart options cannot be processed, e.g.
+            because they do not have the correct type.
+
+    """
+
+    n_samples: int | None = None
+    stopping_maxopt: int | None = None
+    sampling_distribution: Literal["uniform", "triangular"] = "uniform"
+    sampling_method: Literal["sobol", "random", "halton", "latin_hypercube"] = "random"
+    sample: Sequence[PyTree] | None = None
+    mixing_weight_method: (
+        Literal["tiktak", "linear"] | Callable[[int, int, float, float], float]
+    ) = "tiktak"
+    mixing_weight_bounds: tuple[float, float] = (0.1, 0.995)
+    convergence_xtol_rel: float | None = None
+    convergence_max_discoveries: int = 2
+    n_cores: int = 1
+    # TODO: Add more informative type hint for batch_evaluator
+    batch_evaluator: Literal["joblib", "pathos"] | Callable = "joblib"  # type: ignore
+    batch_size: int | None = None
+    seed: int | np.random.Generator | None = None
+    error_handling: Literal["raise", "continue"] | None = None
+    # Deprecated attributes
+    share_optimization: float | None = None
+    convergence_relative_params_tolerance: float | None = None
+    optimization_error_handling: Literal["raise", "continue"] | None = None
+    exploration_error_handling: Literal["raise", "continue"] | None = None
+
+    def __post_init__(self) -> None:
+        _validate_attribute_types_and_values(self)
+
+
+class MultistartOptionsDict(TypedDict):
+    n_samples: NotRequired[int | None]
+    stopping_maxopt: NotRequired[int | None]
+    sampling_distribution: NotRequired[Literal["uniform", "triangular"]]
+    sampling_method: NotRequired[
+        Literal["sobol", "random", "halton", "latin_hypercube"]
+    ]
+    sample: NotRequired[Sequence[PyTree] | None]
+    mixing_weight_method: NotRequired[
+        Literal["tiktak", "linear"] | Callable[[int, int, float, float], float]
+    ]
+    mixing_weight_bounds: NotRequired[tuple[float, float]]
+    convergence_xtol_rel: NotRequired[float | None]
+    convergence_max_discoveries: NotRequired[int]
+    n_cores: NotRequired[int]
+    batch_evaluator: NotRequired[Literal["joblib", "pathos"] | Callable]  # type: ignore
+    batch_size: NotRequired[int | None]
+    seed: NotRequired[int | np.random.Generator | None]
+    error_handling: NotRequired[Literal["raise", "continue"] | None]
+    # Deprecated attributes
+    share_optimization: NotRequired[float | None]
+    convergence_relative_params_tolerance: NotRequired[float | None]
+    optimization_error_handling: NotRequired[Literal["raise", "continue"] | None]
+    exploration_error_handling: NotRequired[Literal["raise", "continue"] | None]
+
+
+def pre_process_multistart(
+    multistart: bool | MultistartOptions | MultistartOptionsDict | None,
+) -> MultistartOptions | None:
+    """Convert all valid types of multistart to a optimagic.MultistartOptions.
+
+    This just harmonizes multiple ways of specifying multistart options into a single
+    format. It performs runime type checks, but it does not check whether multistart
+    options are consistent with other option choices.
+
+    Args:
+        multistart: The user provided multistart options.
+        n_params: The number of parameters in the optimization problem.
+
+    Returns:
+        The multistart options in the optimagic format.
+
+    Raises:
+        InvalidMultistartError: If the multistart options cannot be processed, e.g.
+            because they do not have the correct type.
+
+    """
+    if isinstance(multistart, bool):
+        multistart = MultistartOptions() if multistart else None
+    elif isinstance(multistart, MultistartOptions) or multistart is None:
+        pass
+    else:
+        try:
+            multistart = MultistartOptions(**multistart)
+        except (KeyboardInterrupt, SystemExit):
+            raise
+        except Exception as e:
+            if isinstance(e, InvalidMultistartError):
+                raise e
+            raise InvalidMultistartError(
+                f"Invalid multistart options of type: {type(multistart)}. Multistart "
+                "options must be of type optimagic.MultistartOptions, a dictionary "
+                "with valid keys, None, or a boolean."
+            ) from e
+
+    if multistart is not None:
+        multistart = replace_and_warn_about_deprecated_multistart_options(multistart)
+        # The replace and warn function cannot be typed due to circular imports, but
+        # we know that the return type is MultistartOptions
+        multistart = cast(MultistartOptions, multistart)
+
+    return multistart
+
+
+def _validate_attribute_types_and_values(options: MultistartOptions) -> None:
+    if options.n_samples is not None and (
+        not isinstance(options.n_samples, int) or options.n_samples < 1
+    ):
+        raise InvalidMultistartError(
+            f"Invalid number of samples: {options.n_samples}. Number of samples "
+            "must be a positive integer or None."
+        )
+
+    if options.stopping_maxopt is not None and (
+        not isinstance(options.stopping_maxopt, int) or options.stopping_maxopt < 0
+    ):
+        raise InvalidMultistartError(
+            f"Invalid number of optimizations: {options.stopping_maxopt}. Number of "
+            "optimizations must be a positive integer or None."
+        )
+
+    if (
+        options.n_samples is not None
+        and options.stopping_maxopt is not None
+        and options.n_samples < options.stopping_maxopt
+    ):
+        raise InvalidMultistartError(
+            f"Invalid number of samples: {options.n_samples}. Number of samples "
+            "must be at least as large as the number of optimizations."
+        )
+
+    if options.sampling_distribution not in ("uniform", "triangular"):
+        raise InvalidMultistartError(
+            f"Invalid sampling distribution: {options.sampling_distribution}. Sampling "
+            f"distribution must be one of ('uniform', 'triangular')."
+        )
+
+    if options.sampling_method not in ("sobol", "random", "halton", "latin_hypercube"):
+        raise InvalidMultistartError(
+            f"Invalid sampling method: {options.sampling_method}. Sampling method "
+            f"must be one of ('sobol', 'random', 'halton', 'latin_hypercube')."
+        )
+
+    if not isinstance(options.sample, Sequence | None):
+        raise InvalidMultistartError(
+            f"Invalid sample: {options.sample}. Sample must be a sequence of "
+            "parameters."
+        )
+
+    if not callable(
+        options.mixing_weight_method
+    ) and options.mixing_weight_method not in ("tiktak", "linear"):
+        raise InvalidMultistartError(
+            f"Invalid mixing weight method: {options.mixing_weight_method}. Mixing "
+            "weight method must be Callable or one of ('tiktak', 'linear')."
+        )
+
+    if (
+        not isinstance(options.mixing_weight_bounds, tuple)
+        or len(options.mixing_weight_bounds) != 2
+        or not set(type(x) for x in options.mixing_weight_bounds) <= {int, float}
+    ):
+        raise InvalidMultistartError(
+            f"Invalid mixing weight bounds: {options.mixing_weight_bounds}. Mixing "
+            "weight bounds must be a tuple of two numbers."
+        )
+
+    if options.convergence_xtol_rel is not None and (
+        not isinstance(options.convergence_xtol_rel, int | float)
+        or options.convergence_xtol_rel < 0
+    ):
+        raise InvalidMultistartError(
+            "Invalid relative params tolerance:"
+            f"{options.convergence_xtol_rel}. Relative params "
+            "tolerance must be a number."
+        )
+
+    if (
+        not isinstance(options.convergence_max_discoveries, int | float)
+        or options.convergence_max_discoveries < 1
+    ):
+        raise InvalidMultistartError(
+            f"Invalid max discoveries: {options.convergence_max_discoveries}. Max "
+            "discoveries must be a positive integer or infinity."
+        )
+
+    if not isinstance(options.n_cores, int) or options.n_cores < 1:
+        raise InvalidMultistartError(
+            f"Invalid number of cores: {options.n_cores}. Number of cores "
+            "must be a positive integer."
+        )
+
+    if not callable(options.batch_evaluator) and options.batch_evaluator not in (
+        "joblib",
+        "pathos",
+    ):
+        raise InvalidMultistartError(
+            f"Invalid batch evaluator: {options.batch_evaluator}. Batch evaluator "
+            "must be a Callable or one of 'joblib' or 'pathos'."
+        )
+
+    if options.batch_size is not None and (
+        not isinstance(options.batch_size, int) or options.batch_size < options.n_cores
+    ):
+        raise InvalidMultistartError(
+            f"Invalid batch size: {options.batch_size}. Batch size "
+            "must be a positive integer larger than n_cores, or None."
+        )
+
+    if not isinstance(options.seed, int | np.random.Generator | None):
+        raise InvalidMultistartError(
+            f"Invalid seed: {options.seed}. Seed "
+            "must be an integer, a numpy random generator, or None."
+        )
+
+    if options.error_handling is not None and options.error_handling not in (
+        "raise",
+        "continue",
+    ):
+        raise InvalidMultistartError(
+            f"Invalid error handling: {options.error_handling}. Error handling must be "
+            "'raise' or 'continue'."
+        )
+
+
+# ======================================================================================
+# Internal Options
+# ======================================================================================
+
+
+def _tiktak_weights(
+    iteration: int, n_iterations: int, min_weight: float, max_weight: float
+) -> float:
+    return np.clip(np.sqrt(iteration / n_iterations), min_weight, max_weight)
+
+
+def _linear_weights(
+    iteration: int, n_iterations: int, min_weight: float, max_weight: float
+) -> float:
+    unscaled = iteration / n_iterations
+    span = max_weight - min_weight
+    return min_weight + unscaled * span
+
+
+WEIGHT_FUNCTIONS = {
+    "tiktak": _tiktak_weights,
+    "linear": _linear_weights,
+}
+
+
+@dataclass(frozen=True)
+class InternalMultistartOptions:
+    """Multistart options used internally in optimagic.
+
+    Compared to `MultistartOptions`, this data class has stricter types and combines
+    some of the attributes. It is generated at runtime using a `MultistartOptions`
+    instance and the function `get_internal_multistart_options_from_public`.
+
+    """
+
+    n_samples: int
+    weight_func: Callable[[int, int], float]
+    convergence_xtol_rel: float
+    convergence_max_discoveries: int
+    sampling_distribution: Literal["uniform", "triangular"]
+    sampling_method: Literal["sobol", "random", "halton", "latin_hypercube"]
+    sample: NDArray[np.float64] | None
+    seed: int | np.random.Generator | None
+    n_cores: int
+    # TODO: Add more informative type hint for batch_evaluator
+    batch_evaluator: Callable  # type: ignore
+    batch_size: int
+    error_handling: Literal["raise", "continue"]
+    stopping_maxopt: int
+
+    def __post_init__(self) -> None:
+        must_be_at_least_1 = [
+            "n_samples",
+            "stopping_maxopt",
+            "n_cores",
+            "batch_size",
+            "convergence_max_discoveries",
+        ]
+
+        for attr in must_be_at_least_1:
+            if getattr(self, attr) < 1:
+                raise InvalidMultistartError(f"{attr} must be at least 1.")
+
+        if self.batch_size < self.n_cores:
+            raise InvalidMultistartError("batch_size must be at least n_cores.")
+
+        if self.convergence_xtol_rel < 0:
+            raise InvalidMultistartError("convergence_xtol_rel must be at least 0.")
+
+
+def get_internal_multistart_options_from_public(
+    options: MultistartOptions,
+    params: PyTree,
+    params_to_internal: Callable[[PyTree], NDArray[np.float64]],
+) -> InternalMultistartOptions:
+    """Get internal multistart options from public multistart options.
+
+    Args:
+        options: The pre-processed multistart options.
+        params: The parameters of the optimization problem.
+        params_to_internal: A function that converts parameters to internal parameters.
+
+    Returns:
+        InternalMultistartOptions: The updated options with runtime defaults.
+
+    """
+    x = params_to_internal(params)
+
+    if options.sample is not None:
+        sample = np.array([params_to_internal(x) for x in list(options.sample)])
+        n_samples = len(options.sample)
+    else:
+        sample = None
+        n_samples = options.n_samples  # type: ignore
+
+    batch_size = options.n_cores if options.batch_size is None else options.batch_size
+    batch_evaluator = process_batch_evaluator(options.batch_evaluator)
+
+    if callable(options.mixing_weight_method):
+        weight_func = options.mixing_weight_method
+    else:
+        _weight_method = WEIGHT_FUNCTIONS[options.mixing_weight_method]
+
+    weight_func = partial(
+        _weight_method,
+        min_weight=options.mixing_weight_bounds[0],
+        max_weight=options.mixing_weight_bounds[1],
+    )
+
+    if n_samples is None:
+        if options.stopping_maxopt is None:
+            n_samples = 100 * len(x)
+        else:
+            n_samples = 10 * options.stopping_maxopt
+
+    if options.share_optimization is None:
+        share_optimization = 0.1
+    else:
+        share_optimization = options.share_optimization
+
+    if options.stopping_maxopt is None:
+        stopping_maxopt = max(1, int(share_optimization * n_samples))
+    else:
+        stopping_maxopt = options.stopping_maxopt
+
+    # Set defaults resulting from deprecated attributes
+    if options.error_handling is not None:
+        error_handling = options.error_handling
+    else:
+        error_handling = "continue"
+
+    if options.convergence_xtol_rel is not None:
+        convergence_xtol_rel = options.convergence_xtol_rel
+    else:
+        convergence_xtol_rel = 0.01
+
+    return InternalMultistartOptions(
+        # Attributes taken directly from MultistartOptions
+        convergence_max_discoveries=options.convergence_max_discoveries,
+        n_cores=options.n_cores,
+        sampling_distribution=options.sampling_distribution,
+        sampling_method=options.sampling_method,
+        seed=options.seed,
+        # Updated attributes
+        sample=sample,
+        n_samples=n_samples,
+        weight_func=weight_func,
+        error_handling=error_handling,
+        convergence_xtol_rel=convergence_xtol_rel,
+        stopping_maxopt=stopping_maxopt,
+        batch_evaluator=batch_evaluator,
+        batch_size=batch_size,
+    )
diff --git a/src/optimagic/optimization/optimize.py b/src/optimagic/optimization/optimize.py
index 65fcb6753..4dd3e09e2 100644
--- a/src/optimagic/optimization/optimize.py
+++ b/src/optimagic/optimization/optimize.py
@@ -16,7 +16,6 @@
 import warnings
 from pathlib import Path
 
-from optimagic.batch_evaluators import process_batch_evaluator
 from optimagic.exceptions import (
     InvalidFunctionError,
     InvalidKwargsError,
@@ -40,12 +39,13 @@
     internal_criterion_and_derivative_template,
 )
 from optimagic.optimization.multistart import (
-    WEIGHT_FUNCTIONS,
     run_multistart_optimization,
 )
+from optimagic.optimization.multistart_options import (
+    get_internal_multistart_options_from_public,
+)
 from optimagic.optimization.optimization_logging import log_scheduled_steps_and_get_ids
 from optimagic.optimization.optimize_result import OptimizeResult
-from optimagic.optimization.process_multistart_sample import process_multistart_sample
 from optimagic.optimization.process_results import process_internal_optimizer_result
 from optimagic.parameters.bounds import Bounds
 from optimagic.parameters.conversion import (
@@ -75,7 +75,6 @@ def maximize(
     error_penalty=None,
     scaling=False,
     multistart=False,
-    multistart_options=None,
     collect_history=True,
     skip_checks=False,
     # scipy aliases
@@ -101,6 +100,7 @@ def maximize(
     soft_lower_bounds=None,
     soft_upper_bounds=None,
     scaling_options=None,
+    multistart_options=None,
 ):
     """Maximize fun using algorithm subject to constraints.
 
@@ -136,7 +136,6 @@ def maximize(
         error_penalty=error_penalty,
         scaling=scaling,
         multistart=multistart,
-        multistart_options=multistart_options,
         collect_history=collect_history,
         skip_checks=skip_checks,
         # scipy aliases
@@ -162,6 +161,7 @@ def maximize(
         soft_lower_bounds=soft_lower_bounds,
         soft_upper_bounds=soft_upper_bounds,
         scaling_options=scaling_options,
+        multistart_options=multistart_options,
     )
     return _optimize(problem)
 
@@ -186,7 +186,6 @@ def minimize(
     error_penalty=None,
     scaling=False,
     multistart=False,
-    multistart_options=None,
     collect_history=True,
     skip_checks=False,
     # scipy aliases
@@ -212,6 +211,7 @@ def minimize(
     soft_lower_bounds=None,
     soft_upper_bounds=None,
     scaling_options=None,
+    multistart_options=None,
 ):
     """Minimize criterion using algorithm subject to constraints.
 
@@ -248,7 +248,6 @@ def minimize(
         error_penalty=error_penalty,
         scaling=scaling,
         multistart=multistart,
-        multistart_options=multistart_options,
         collect_history=collect_history,
         skip_checks=skip_checks,
         # scipy aliases
@@ -274,6 +273,7 @@ def minimize(
         soft_lower_bounds=soft_lower_bounds,
         soft_upper_bounds=soft_upper_bounds,
         scaling_options=scaling_options,
+        multistart_options=multistart_options,
     )
     return _optimize(problem)
 
@@ -347,7 +347,7 @@ def _optimize(problem: OptimizationProblem) -> OptimizeResult:
         primary_key=problem.algo_info.primary_criterion_entry,
         scaling=problem.scaling,
         derivative_eval=used_deriv,
-        add_soft_bounds=problem.multistart,
+        add_soft_bounds=problem.multistart is not None,
     )
 
     # ==================================================================================
@@ -373,7 +373,7 @@ def _optimize(problem: OptimizationProblem) -> OptimizeResult:
     # Do some things that require internal parameters or bounds
     # ==================================================================================
 
-    if converter.has_transforming_constraints and problem.multistart:
+    if converter.has_transforming_constraints and problem.multistart is not None:
         raise NotImplementedError(
             "multistart optimizations are not yet compatible with transforming "
             "constraints."
@@ -453,7 +453,7 @@ def _optimize(problem: OptimizationProblem) -> OptimizeResult:
     # ==================================================================================
     # Do actual optimization
     # ==================================================================================
-    if not problem.multistart:
+    if problem.multistart is None:
         steps = [{"type": "optimization", "name": "optimization"}]
 
         step_ids = log_scheduled_steps_and_get_ids(
@@ -464,10 +464,9 @@ def _optimize(problem: OptimizationProblem) -> OptimizeResult:
 
         raw_res = internal_algorithm(**problem_functions, x=x, step_id=step_ids[0])
     else:
-        multistart_options = _fill_multistart_options_with_defaults(
-            options=problem.multistart_options,
+        multistart_options = get_internal_multistart_options_from_public(
+            options=problem.multistart,
             params=problem.params,
-            x=x,
             params_to_internal=converter.params_to_internal,
         )
 
@@ -603,53 +602,3 @@ def _fill_numdiff_options_with_defaults(numdiff_options, lower_bounds, upper_bou
 
     numdiff_options = {**default_numdiff_options, **numdiff_options}
     return numdiff_options
-
-
-def _setdefault(candidate, default):
-    out = default if candidate is None else candidate
-    return out
-
-
-def _fill_multistart_options_with_defaults(options, params, x, params_to_internal):
-    """Fill options for multistart optimization with defaults."""
-    defaults = {
-        "sample": None,
-        "n_samples": 10 * len(x),
-        "share_optimizations": 0.1,
-        "sampling_distribution": "uniform",
-        "sampling_method": "sobol" if len(x) <= 200 else "random",
-        "mixing_weight_method": "tiktak",
-        "mixing_weight_bounds": (0.1, 0.995),
-        "convergence_relative_params_tolerance": 0.01,
-        "convergence_max_discoveries": 2,
-        "n_cores": 1,
-        "batch_evaluator": "joblib",
-        "seed": None,
-        "exploration_error_handling": "continue",
-        "optimization_error_handling": "continue",
-    }
-
-    options = {k.replace(".", "_"): v for k, v in options.items()}
-    out = {**defaults, **options}
-
-    if "batch_size" not in out:
-        out["batch_size"] = out["n_cores"]
-    else:
-        if out["batch_size"] < out["n_cores"]:
-            raise ValueError("batch_size must be at least as large as n_cores.")
-
-    out["batch_evaluator"] = process_batch_evaluator(out["batch_evaluator"])
-
-    if isinstance(out["mixing_weight_method"], str):
-        out["mixing_weight_method"] = WEIGHT_FUNCTIONS[out["mixing_weight_method"]]
-
-    if out["sample"] is not None:
-        out["sample"] = process_multistart_sample(
-            out["sample"], params, params_to_internal
-        )
-        out["n_samples"] = len(out["sample"])
-
-    out["n_optimizations"] = max(1, int(out["n_samples"] * out["share_optimizations"]))
-    del out["share_optimizations"]
-
-    return out
diff --git a/src/optimagic/optimization/optimize_result.py b/src/optimagic/optimization/optimize_result.py
index afa4d3687..16b6eb47a 100644
--- a/src/optimagic/optimization/optimize_result.py
+++ b/src/optimagic/optimization/optimize_result.py
@@ -1,6 +1,6 @@
 import warnings
 from dataclasses import dataclass, field
-from typing import Any, Dict
+from typing import Any, Dict, Optional
 
 import numpy as np
 import pandas as pd
@@ -60,7 +60,7 @@ class OptimizeResult:
 
     convergence_report: Dict | None = None
 
-    multistart_info: Dict | None = None
+    multistart_info: Optional["MultistartInfo"] = None
     algorithm_output: Dict = field(default_factory=dict)
 
     # ==================================================================================
@@ -198,6 +198,32 @@ def to_pickle(self, path):
         to_pickle(self, path=path)
 
 
+@dataclass(frozen=True)
+class MultistartInfo:
+    """Information about the multistart optimization.
+
+    Attributes:
+        start_parameters: List of start parameters for each optimization.
+        local_optima: List of optimization results.
+        exploration_sample: List of parameters used for exploration.
+        exploration_results: List of function values corresponding to exploration.
+        n_optimizations: Number of local optimizations that were run.
+
+    """
+
+    start_parameters: list[PyTree]
+    local_optima: list[OptimizeResult]
+    exploration_sample: list[PyTree]
+    exploration_results: list[float]
+
+    def __getitem__(self, item):
+        return getattr(self, item)
+
+    @property
+    def n_optimizations(self) -> int:
+        return len(self.local_optima)
+
+
 def _format_convergence_report(report, algorithm):
     report = pd.DataFrame.from_dict(report)
     columns = ["one_step", "five_steps"]
diff --git a/src/optimagic/optimization/process_multistart_sample.py b/src/optimagic/optimization/process_multistart_sample.py
deleted file mode 100644
index 274215998..000000000
--- a/src/optimagic/optimization/process_multistart_sample.py
+++ /dev/null
@@ -1,72 +0,0 @@
-import numpy as np
-import pandas as pd
-
-
-def process_multistart_sample(raw_sample, params, params_to_internal):
-    """Process a user provided multistart sample.
-
-    Args:
-        raw_sample (list, pd.DataFrame or np.ndarray): A user provided sample of
-            external start parameters.
-        params (pytree): User provided start parameters.
-        params_to_internal (callable): A function that converts external parameters
-            to internal ones.
-
-
-    Returns:
-        np.ndarray: 2d numpy array where each row is an internal parameter vector.
-
-    """
-    is_df_params = isinstance(params, pd.DataFrame) and "value" in params
-    is_np_params = isinstance(params, np.ndarray) and params.ndim == 1
-
-    if isinstance(raw_sample, pd.DataFrame):
-        if not is_df_params:
-            msg = (
-                "User provided multistart samples can only be a DataFrame if "
-                "params is a DataFrame with 'value' column."
-            )
-            raise ValueError(msg)
-        elif not raw_sample.columns.equals(params.index):
-            msg = (
-                "If you provide a custom sample as DataFrame the columns of that "
-                "DataFrame and the index of params must be equal."
-            )
-            raise ValueError(msg)
-
-        list_sample = [params.assign(value=row) for _, row in raw_sample.iterrows()]
-
-    elif isinstance(raw_sample, np.ndarray):
-        if not is_np_params:
-            msg = (
-                "User provided multistart samples can only be a numpy array if params "
-                "is a 1d numpy array."
-            )
-            raise ValueError(msg)
-        elif raw_sample.ndim != 2:
-            msg = (
-                "If user provided multistart samples are a numpy array, the array "
-                "must be two dimensional."
-            )
-            raise ValueError(msg)
-        elif raw_sample.shape[1] != len(params):
-            msg = (
-                "If user provided multistart samples are a numpy array, the number of "
-                "columns must be equal to the number of parameters."
-            )
-            raise ValueError(msg)
-
-        list_sample = list(raw_sample)
-
-    elif not isinstance(raw_sample, (list, tuple)):
-        msg = (
-            "User provided multistart samples must be a list, tuple, numpy array or "
-            "DataFrame."
-        )
-        raise TypeError(msg)
-    else:
-        list_sample = list(raw_sample)
-
-    sample = np.array([params_to_internal(x) for x in list_sample])
-
-    return sample
diff --git a/src/optimagic/optimization/process_results.py b/src/optimagic/optimization/process_results.py
index c1f09516c..e69c7d32f 100644
--- a/src/optimagic/optimization/process_results.py
+++ b/src/optimagic/optimization/process_results.py
@@ -1,7 +1,7 @@
 import numpy as np
 
 from optimagic.optimization.convergence_report import get_convergence_report
-from optimagic.optimization.optimize_result import OptimizeResult
+from optimagic.optimization.optimize_result import MultistartInfo, OptimizeResult
 from optimagic.parameters.conversion import aggregate_func_output_to_value
 
 
@@ -143,13 +143,12 @@ def _process_multistart_info(info, converter, primary_key, fixed_kwargs, skip_ch
     else:
         exploration_res = [switch_sign(res) for res in info["exploration_results"]]
 
-    out = {
-        "start_parameters": starts,
-        "local_optima": optima,
-        "exploration_sample": sample,
-        "exploration_results": exploration_res,
-    }
-    return out
+    return MultistartInfo(
+        start_parameters=starts,
+        local_optima=optima,
+        exploration_sample=sample,
+        exploration_results=exploration_res,
+    )
 
 
 def _dummy_result_from_traceback(candidate, fixed_kwargs):  # noqa: ARG001
diff --git a/src/optimagic/parameters/scaling.py b/src/optimagic/parameters/scaling.py
index fab36aee3..03b6c2b62 100644
--- a/src/optimagic/parameters/scaling.py
+++ b/src/optimagic/parameters/scaling.py
@@ -18,12 +18,19 @@ class ScalingOptions:
         magnitude: A factor by which the scaled parameters are multiplied to adjust
             their magnitude. Must be a positive number. Default is 1.0.
 
+    Raises:
+        InvalidScalingError: If scaling options cannot be processed, e.g. because they
+            do not have the correct type.
+
     """
 
     method: Literal["start_values", "bounds"] = "start_values"
     clipping_value: float = 0.1
     magnitude: float = 1.0
 
+    def __post_init__(self) -> None:
+        _validate_attribute_types_and_values(self)
+
 
 class ScalingOptionsDict(TypedDict):
     method: NotRequired[Literal["start_values", "bounds"]]
@@ -47,8 +54,8 @@ def pre_process_scaling(
         The scaling options in the optimagic format.
 
     Raises:
-        InvalidScalingOptionsError: If scaling options cannot be processed, e.g. because
-            they do not have the correct type.
+        InvalidScalingError: If scaling options cannot be processed, e.g. because they
+            do not have the correct type.
 
     """
     if isinstance(scaling, bool):
@@ -61,6 +68,8 @@ def pre_process_scaling(
         except (KeyboardInterrupt, SystemExit):
             raise
         except Exception as e:
+            if isinstance(e, InvalidScalingError):
+                raise e
             raise InvalidScalingError(
                 f"Invalid scaling options of type: {type(scaling)}. Scaling options "
                 "must be of type optimagic.ScalingOptions, a dictionary with a subset "
@@ -68,23 +77,27 @@ def pre_process_scaling(
                 "boolean."
             ) from e
 
-    if isinstance(scaling, ScalingOptions):
-        if scaling.method not in ("start_values", "bounds"):
-            raise InvalidScalingError(
-                f"Invalid scaling method: {scaling.method}. Valid methods are "
-                "'start_values' and 'bounds'."
-            )
-
-        if not isinstance(scaling.clipping_value, (int, float)):
-            raise InvalidScalingError(
-                f"Invalid clipping value: {scaling.clipping_value}. Clipping value "
-                "must be a number."
-            )
+    return scaling
 
-        if not isinstance(scaling.magnitude, (int, float)) or scaling.magnitude <= 0:
-            raise InvalidScalingError(
-                f"Invalid scaling magnitude: {scaling.magnitude}. Scaling magnitude "
-                "must be a positive number."
-            )
 
-    return scaling
+def _validate_attribute_types_and_values(options: ScalingOptions) -> None:
+    if options.method not in ("start_values", "bounds"):
+        raise InvalidScalingError(
+            f"Invalid scaling method: {options.method}. Valid methods are "
+            "'start_values' and 'bounds'."
+        )
+
+    if (
+        not isinstance(options.clipping_value, int | float)
+        or options.clipping_value <= 0
+    ):
+        raise InvalidScalingError(
+            f"Invalid clipping value: {options.clipping_value}. Clipping value "
+            "must be a positive number."
+        )
+
+    if not isinstance(options.magnitude, int | float) or options.magnitude <= 0:
+        raise InvalidScalingError(
+            f"Invalid scaling magnitude: {options.magnitude}. Scaling magnitude "
+            "must be a positive number."
+        )
diff --git a/tests/optimagic/optimization/test_history_collection.py b/tests/optimagic/optimization/test_history_collection.py
index 29d6851a7..b98b69578 100644
--- a/tests/optimagic/optimization/test_history_collection.py
+++ b/tests/optimagic/optimization/test_history_collection.py
@@ -34,7 +34,7 @@ def test_history_collection_with_parallelization(algorithm, tmp_path):
         params=np.arange(5),
         algorithm=algorithm,
         bounds=Bounds(lower=lb, upper=ub),
-        algo_options={"n_cores": 2, "stopping.max_iterations": 3},
+        algo_options={"n_cores": 2, "stopping_maxiter": 3},
         logging=logging,
         log_options={"if_database_exists": "replace", "fast_logging": True},
     ).history
diff --git a/tests/optimagic/optimization/test_multistart.py b/tests/optimagic/optimization/test_multistart.py
index b170e3c7a..8f22e72ac 100644
--- a/tests/optimagic/optimization/test_multistart.py
+++ b/tests/optimagic/optimization/test_multistart.py
@@ -5,9 +5,7 @@
 import pytest
 from numpy.testing import assert_array_almost_equal as aaae
 from optimagic.optimization.multistart import (
-    _linear_weights,
-    _tiktak_weights,
-    draw_exploration_sample,
+    _draw_exploration_sample,
     get_batched_optimization_sample,
     run_explorations,
     update_convergence_state,
@@ -42,13 +40,13 @@ def test_draw_exploration_sample(dist, rule, lower, upper):
 
     for _ in range(2):
         results.append(
-            draw_exploration_sample(
+            _draw_exploration_sample(
                 x=np.ones_like(lower) * 0.5,
                 lower=lower,
                 upper=upper,
                 n_samples=3,
-                sampling_distribution=dist,
-                sampling_method=rule,
+                distribution=dist,
+                method=rule,
                 seed=1234,
             )
         )
@@ -92,7 +90,7 @@ def _dummy(x, **kwargs):
 def test_get_batched_optimization_sample():
     calculated = get_batched_optimization_sample(
         sorted_sample=np.arange(12).reshape(6, 2),
-        n_optimizations=5,
+        stopping_maxopt=5,
         batch_size=4,
     )
     expected = [[[0, 1], [2, 3], [4, 5], [6, 7]], [[8, 9]]]
@@ -108,17 +106,6 @@ def test_get_batched_optimization_sample():
             assert calc_entry.tolist() == exp_entry
 
 
-def test_linear_weights():
-    calculated = _linear_weights(5, 10, 0.4, 0.8)
-    expected = 0.6
-    assert np.allclose(calculated, expected)
-
-
-def test_tiktak_weights():
-    assert np.allclose(0.3, _tiktak_weights(0, 10, 0.3, 0.8))
-    assert np.allclose(0.8, _tiktak_weights(10, 10, 0.3, 0.8))
-
-
 @pytest.fixture()
 def current_state():
     state = {
diff --git a/tests/optimagic/optimization/test_multistart_options.py b/tests/optimagic/optimization/test_multistart_options.py
new file mode 100644
index 000000000..20785fdb9
--- /dev/null
+++ b/tests/optimagic/optimization/test_multistart_options.py
@@ -0,0 +1,162 @@
+import numpy as np
+import pytest
+from optimagic.exceptions import InvalidMultistartError
+from optimagic.optimization.multistart_options import (
+    MultistartOptions,
+    _linear_weights,
+    _tiktak_weights,
+    get_internal_multistart_options_from_public,
+    pre_process_multistart,
+)
+
+
+def test_pre_process_multistart_trivial_case():
+    multistart = MultistartOptions(n_samples=10, convergence_max_discoveries=55)
+    got = pre_process_multistart(multistart)
+    assert got == multistart
+
+
+def test_pre_process_multistart_none_case():
+    assert pre_process_multistart(None) is None
+
+
+def test_pre_process_multistart_false_case():
+    assert pre_process_multistart(False) is None
+
+
+def test_pre_process_multistart_dict_case():
+    got = pre_process_multistart(
+        multistart={
+            "n_samples": 10,
+            "convergence_max_discoveries": 55,
+        }
+    )
+    assert got == MultistartOptions(
+        n_samples=10,
+        convergence_max_discoveries=55,
+    )
+
+
+def test_pre_process_multistart_invalid_type():
+    with pytest.raises(InvalidMultistartError, match="Invalid multistart options"):
+        pre_process_multistart(multistart="invalid")
+
+
+def test_pre_process_multistart_invalid_dict_key():
+    with pytest.raises(InvalidMultistartError, match="Invalid multistart options"):
+        pre_process_multistart(multistart={"invalid": "invalid"})
+
+
+def test_pre_process_multistart_invalid_dict_value():
+    with pytest.raises(InvalidMultistartError, match="Invalid number of samples"):
+        pre_process_multistart(multistart={"n_samples": "invalid"})
+
+
+@pytest.mark.parametrize("value", ["invalid", -1])
+def test_multistart_options_invalid_n_samples_value(value):
+    with pytest.raises(InvalidMultistartError, match="Invalid number of samples"):
+        MultistartOptions(n_samples=value)
+
+
+@pytest.mark.parametrize("value", ["invalid", -1])
+def test_multistart_options_invalid_stopping_maxopt(value):
+    with pytest.raises(InvalidMultistartError, match="Invalid number of optimizations"):
+        MultistartOptions(stopping_maxopt=value)
+
+
+def test_multistart_options_stopping_maxopt_less_than_n_samples():
+    with pytest.raises(InvalidMultistartError, match="Invalid number of samples"):
+        MultistartOptions(n_samples=1, stopping_maxopt=2)
+
+
+def test_multistart_options_invalid_sampling_distribution():
+    with pytest.raises(InvalidMultistartError, match="Invalid sampling distribution"):
+        MultistartOptions(sampling_distribution="invalid")
+
+
+def test_multistart_options_invalid_sampling_method():
+    with pytest.raises(InvalidMultistartError, match="Invalid sampling method"):
+        MultistartOptions(sampling_method="invalid")
+
+
+def test_multistart_options_invalid_mixing_weight_method():
+    with pytest.raises(InvalidMultistartError, match="Invalid mixing weight method"):
+        MultistartOptions(mixing_weight_method="invalid")
+
+
+@pytest.mark.parametrize("value", [("a", "b"), (1, 2, 3), {"a": 1.0, "b": 3.0}])
+def test_multistart_options_invalid_mixing_weight_bounds(value):
+    with pytest.raises(InvalidMultistartError, match="Invalid mixing weight bounds"):
+        MultistartOptions(mixing_weight_bounds=value)
+
+
+def test_multistart_options_invalid_convergence_xtol_rel():
+    with pytest.raises(InvalidMultistartError, match="Invalid relative params"):
+        MultistartOptions(convergence_xtol_rel="invalid")
+
+
+@pytest.mark.parametrize("value", ["invalid", -1])
+def test_multistart_options_invalid_convergence_max_discoveries(value):
+    with pytest.raises(InvalidMultistartError, match="Invalid max discoveries"):
+        MultistartOptions(convergence_max_discoveries=value)
+
+
+@pytest.mark.parametrize("value", ["invalid", -1])
+def test_multistart_options_invalid_n_cores(value):
+    with pytest.raises(InvalidMultistartError, match="Invalid number of cores"):
+        MultistartOptions(n_cores=value)
+
+
+@pytest.mark.parametrize("value", ["invalid", -1])
+def test_multistart_options_invalid_batch_size(value):
+    with pytest.raises(InvalidMultistartError, match="Invalid batch size"):
+        MultistartOptions(batch_size=value)
+
+
+def test_multistart_options_batch_size_smaller_than_n_cores():
+    with pytest.raises(InvalidMultistartError, match="Invalid batch size"):
+        MultistartOptions(batch_size=1, n_cores=2)
+
+
+def test_multistart_options_invalid_batch_evaluator():
+    with pytest.raises(InvalidMultistartError, match="Invalid batch evaluator"):
+        MultistartOptions(batch_evaluator="invalid")
+
+
+def test_multistart_options_invalid_seed():
+    with pytest.raises(InvalidMultistartError, match="Invalid seed"):
+        MultistartOptions(seed="invalid")
+
+
+def test_multistart_options_invalid_error_handling():
+    with pytest.raises(InvalidMultistartError, match="Invalid error handling"):
+        MultistartOptions(error_handling="invalid")
+
+
+def test_linear_weights():
+    calculated = _linear_weights(5, 10, 0.4, 0.8)
+    expected = 0.6
+    assert np.allclose(calculated, expected)
+
+
+def test_tiktak_weights():
+    assert np.allclose(0.3, _tiktak_weights(0, 10, 0.3, 0.8))
+    assert np.allclose(0.8, _tiktak_weights(10, 10, 0.3, 0.8))
+
+
+def test_get_internal_multistart_options_from_public_defaults():
+    options = MultistartOptions()
+
+    got = get_internal_multistart_options_from_public(
+        options,
+        params=np.arange(5),
+        params_to_internal=lambda x: x,
+    )
+
+    assert got.convergence_xtol_rel == 0.01
+    assert got.convergence_max_discoveries == options.convergence_max_discoveries
+    assert got.n_cores == options.n_cores
+    assert got.error_handling == "continue"
+    assert got.n_samples == 500
+    assert got.stopping_maxopt == 50
+    assert got.batch_size == 1
diff --git a/tests/optimagic/optimization/test_process_multistart_sample.py b/tests/optimagic/optimization/test_process_multistart_sample.py
deleted file mode 100644
index f610172c4..000000000
--- a/tests/optimagic/optimization/test_process_multistart_sample.py
+++ /dev/null
@@ -1,25 +0,0 @@
-import numpy as np
-import pandas as pd
-import pytest
-from numpy.testing import assert_array_almost_equal as aaae
-from optimagic.optimization.process_multistart_sample import process_multistart_sample
-
-samples = [
-    (
-        pd.DataFrame(np.ones((2, 3)), columns=["a", "b", "c"]),
-        pd.Series([1, 2, 3], index=["a", "b", "c"], name="value").to_frame(),
-        lambda x: x["value"].to_numpy(),
-    ),
-    (
-        np.ones((2, 3)),
-        np.array([1, 2, 3]),
-        lambda x: x,
-    ),
-]
-
-
-@pytest.mark.parametrize("sample, x, to_internal", samples)
-def test_process_multistart_sample(sample, x, to_internal):
-    calculated = process_multistart_sample(sample, x, to_internal)
-    expeceted = np.ones((2, 3))
-    aaae(calculated, expeceted)
diff --git a/tests/optimagic/optimization/test_with_multistart.py b/tests/optimagic/optimization/test_with_multistart.py
index d83688fdf..a1153d9dd 100644
--- a/tests/optimagic/optimization/test_with_multistart.py
+++ b/tests/optimagic/optimization/test_with_multistart.py
@@ -1,6 +1,7 @@
 from itertools import product
 
 import numpy as np
+import optimagic as om
 import pandas as pd
 import pytest
 from numpy.testing import assert_array_almost_equal as aaae
@@ -52,63 +53,62 @@ def test_multistart_minimize_with_sum_of_squares_at_defaults(
 
     assert hasattr(res, "multistart_info")
     ms_info = res.multistart_info
-    assert len(ms_info["exploration_sample"]) == 40
-    assert len(ms_info["exploration_results"]) == 40
-    assert all(isinstance(entry, float) for entry in ms_info["exploration_results"])
-    assert all(isinstance(entry, OptimizeResult) for entry in ms_info["local_optima"])
-    assert all(isinstance(entry, pd.DataFrame) for entry in ms_info["start_parameters"])
+    assert len(ms_info.exploration_sample) == 400
+    assert len(ms_info.exploration_results) == 400
+    assert all(isinstance(entry, float) for entry in ms_info.exploration_results)
+    assert all(isinstance(entry, OptimizeResult) for entry in ms_info.local_optima)
+    assert all(isinstance(entry, pd.DataFrame) for entry in ms_info.start_parameters)
     assert np.allclose(res.fun, 0)
     aaae(res.params["value"], np.zeros(4))
 
 
 def test_multistart_with_existing_sample(params):
-    sample = pd.DataFrame(
-        np.arange(20).reshape(5, 4) / 10,
-        columns=params.index,
-    )
-    options = {"sample": sample}
+    sample = [params.assign(value=x) for x in np.arange(20).reshape(5, 4) / 10]
+    options = om.MultistartOptions(sample=sample)
 
     res = minimize(
         fun=sos_dict_criterion,
         params=params,
         algorithm="scipy_lbfgsb",
-        multistart=True,
-        multistart_options=options,
+        multistart=options,
     )
 
-    calc_sample = _params_list_to_aray(res.multistart_info["exploration_sample"])
-    aaae(calc_sample, options["sample"])
+    assert all(
+        got.equals(expected)
+        for expected, got in zip(
+            sample, res.multistart_info.exploration_sample, strict=False
+        )
+    )
 
 
 def test_convergence_via_max_discoveries_works(params):
-    options = {
-        "convergence_relative_params_tolerance": np.inf,
-        "convergence_max_discoveries": 2,
-    }
+    options = om.MultistartOptions(
+        convergence_xtol_rel=np.inf,
+        convergence_max_discoveries=2,
+    )
 
     res = maximize(
         fun=switch_sign(sos_dict_criterion),
         params=params,
         algorithm="scipy_lbfgsb",
-        multistart=True,
-        multistart_options=options,
+        multistart=options,
     )
 
-    assert len(res.multistart_info["local_optima"]) == 2
+    assert len(res.multistart_info.local_optima) == 2
 
 
 def test_steps_are_logged_as_skipped_if_convergence(params):
-    options = {
-        "convergence_relative_params_tolerance": np.inf,
-        "convergence_max_discoveries": 2,
-    }
+    options = om.MultistartOptions(
+        n_samples=10 * len(params),
+        convergence_xtol_rel=np.inf,
+        convergence_max_discoveries=2,
+    )
 
     minimize(
         fun=sos_dict_criterion,
         params=params,
         algorithm="scipy_lbfgsb",
-        multistart=True,
-        multistart_options=options,
+        multistart=options,
         logging="logging.db",
     )
 
@@ -118,14 +118,16 @@ def test_steps_are_logged_as_skipped_if_convergence(params):
 
 
 def test_all_steps_occur_in_optimization_iterations_if_no_convergence(params):
-    options = {"convergence_max_discoveries": np.inf}
+    options = om.MultistartOptions(
+        convergence_max_discoveries=np.inf,
+        n_samples=10 * len(params),
+    )
 
     minimize(
         fun=sos_dict_criterion,
         params=params,
         algorithm="scipy_lbfgsb",
-        multistart=True,
-        multistart_options=options,
+        multistart=options,
         logging="logging.db",
     )
 
@@ -165,11 +167,6 @@ def test_error_is_raised_with_transforming_constraints(params):
         )
 
 
-def _params_list_to_aray(params_list):
-    data = [params["value"].tolist() for params in params_list]
-    return np.array(data)
-
-
 def test_multistart_with_numpy_params():
     res = minimize(
         fun=lambda params: params @ params,
@@ -182,6 +179,20 @@ def test_multistart_with_numpy_params():
     aaae(res.params, np.zeros(5))
 
 
+def test_multistart_with_rng_seed():
+    rng = np.random.default_rng(12345)
+
+    res = minimize(
+        fun=lambda params: params @ params,
+        params=np.arange(5),
+        algorithm="scipy_lbfgsb",
+        bounds=Bounds(soft_lower=np.full(5, -10), soft_upper=np.full(5, 10)),
+        multistart=om.MultistartOptions(seed=rng),
+    )
+
+    aaae(res.params, np.zeros(5))
+
+
 def test_with_invalid_bounds():
     with pytest.raises(ValueError):
         minimize(
@@ -230,12 +241,11 @@ def ackley(x):
     minimize(
         **kwargs,
         algorithm="scipy_lbfgsb",
-        multistart=True,
-        multistart_options={
-            "n_samples": 200,
-            "share_optimizations": 0.1,
-            "convergence_max_discoveries": 10,
-        },
+        multistart=om.MultistartOptions(
+            n_samples=200,
+            stopping_maxopt=20,
+            convergence_max_discoveries=10,
+        ),
     )
 
 
@@ -245,8 +255,37 @@ def test_multistart_with_least_squares_optimizers():
         params=np.array([-1, 1.0]),
         bounds=Bounds(soft_lower=np.full(2, -10), soft_upper=np.full(2, 10)),
         algorithm="scipy_ls_trf",
-        multistart=True,
-        multistart_options={"n_samples": 3, "share_optimizations": 1.0},
+        multistart=om.MultistartOptions(n_samples=3, stopping_maxopt=3),
     )
 
     aaae(est.params, np.zeros(2))
+
+
+def test_with_ackley_using_dict_options():
+    def ackley(x):
+        out = (
+            -20 * np.exp(-0.2 * np.sqrt(np.mean(x**2)))
+            - np.exp(np.mean(np.cos(2 * np.pi * x)))
+            + 20
+            + np.exp(1)
+        )
+        return out
+
+    dim = 5
+
+    kwargs = {
+        "fun": ackley,
+        "params": np.full(dim, -10),
+        "bounds": Bounds(lower=np.full(dim, -32), upper=np.full(dim, 32)),
+        "algo_options": {"stopping.maxfun": 1000},
+    }
+
+    minimize(
+        **kwargs,
+        algorithm="scipy_lbfgsb",
+        multistart={
+            "n_samples": 200,
+            "stopping_maxopt": 20,
+            "convergence_max_discoveries": 10,
+        },
+    )
diff --git a/tests/optimagic/parameters/test_scaling.py b/tests/optimagic/parameters/test_scaling.py
index 56b34ac6c..2ae0395f0 100644
--- a/tests/optimagic/parameters/test_scaling.py
+++ b/tests/optimagic/parameters/test_scaling.py
@@ -1,6 +1,9 @@
 import pytest
 from optimagic.exceptions import InvalidScalingError
-from optimagic.parameters.scaling import ScalingOptions, pre_process_scaling
+from optimagic.parameters.scaling import (
+    ScalingOptions,
+    pre_process_scaling,
+)
 
 
 def test_pre_process_scaling_trivial_case():
@@ -39,25 +42,30 @@ def test_pre_process_scaling_invalid_type():
 
 
 def test_pre_process_scaling_invalid_dict_key():
-    with pytest.raises(InvalidScalingError, match="Invalid scaling options"):
+    with pytest.raises(InvalidScalingError, match="Invalid scaling options of type:"):
         pre_process_scaling(scaling={"wrong_key": "start_values"})
 
 
-def test_pre_process_scaling_invalid_method_value():
+def test_pre_process_scaling_invalid_dict_value():
+    with pytest.raises(InvalidScalingError, match="Invalid clipping value:"):
+        pre_process_scaling(scaling={"clipping_value": "invalid"})
+
+
+def test_scaling_options_invalid_method_value():
     with pytest.raises(InvalidScalingError, match="Invalid scaling method:"):
-        pre_process_scaling(scaling={"method": "invalid"})
+        ScalingOptions(method="invalid")
 
 
-def test_pre_process_scaling_invalid_clipping_value_type():
+def test_scaling_options_invalid_clipping_value_type():
     with pytest.raises(InvalidScalingError, match="Invalid clipping value:"):
-        pre_process_scaling(scaling={"clipping_value": "invalid"})
+        ScalingOptions(clipping_value="invalid")
 
 
-def test_pre_process_scaling_invalid_magnitude_value_type():
+def test_scaling_options_invalid_magnitude_value_type():
     with pytest.raises(InvalidScalingError, match="Invalid scaling magnitude:"):
-        pre_process_scaling(scaling={"magnitude": "invalid"})
+        ScalingOptions(magnitude="invalid")
 
 
-def test_pre_process_scaling_invalid_magnitude_value_range():
+def test_scaling_options_invalid_magnitude_value_range():
     with pytest.raises(InvalidScalingError, match="Invalid scaling magnitude:"):
-        pre_process_scaling(scaling={"magnitude": -1})
+        ScalingOptions(magnitude=-1)
diff --git a/tests/optimagic/test_deprecations.py b/tests/optimagic/test_deprecations.py
index 6242eeae6..f87062359 100644
--- a/tests/optimagic/test_deprecations.py
+++ b/tests/optimagic/test_deprecations.py
@@ -531,3 +531,73 @@ def test_old_scaling_options_are_deprecated_in_maximize():
             algorithm="scipy_lbfgsb",
             scaling_options={"method": "start_values", "magnitude": 1},
         )
+
+
+def test_old_multistart_options_are_deprecated_in_minimize():
+    msg = "Specifying multistart options via the argument `multistart_options` is"
+    with pytest.warns(FutureWarning, match=msg):
+        om.minimize(
+            lambda x: x @ x,
+            np.arange(3),
+            algorithm="scipy_lbfgsb",
+            multistart_options={"n_samples": 10},
+        )
+
+
+def test_old_multistart_options_are_deprecated_in_maximize():
+    msg = "Specifying multistart options via the argument `multistart_options` is"
+    with pytest.warns(FutureWarning, match=msg):
+        om.maximize(
+            lambda x: -x @ x,
+            np.arange(3),
+            algorithm="scipy_lbfgsb",
+            multistart_options={"n_samples": 10},
+        )
+
+
+def test_multistart_option_share_optimization_option_is_deprecated():
+    msg = "The share_optimization option is deprecated and will be removed in"
+    with pytest.warns(FutureWarning, match=msg):
+        om.minimize(
+            lambda x: x @ x,
+            np.arange(3),
+            algorithm="scipy_lbfgsb",
+            bounds=om.Bounds(lower=np.full(3, -1), upper=np.full(3, 2)),
+            multistart={"share_optimization": 0.1},
+        )
+
+
+def test_multistart_option_convergence_relative_params_tolerance_option_is_deprecated():
+    msg = "The convergence_relative_params_tolerance option is deprecated and will"
+    with pytest.warns(FutureWarning, match=msg):
+        om.minimize(
+            lambda x: x @ x,
+            np.arange(3),
+            algorithm="scipy_lbfgsb",
+            bounds=om.Bounds(lower=np.full(3, -1), upper=np.full(3, 2)),
+            multistart={"convergence_relative_params_tolerance": 0.01},
+        )
+
+
+def test_multistart_option_optimization_error_handling_option_is_deprecated():
+    msg = "The optimization_error_handling option is deprecated and will be removed"
+    with pytest.warns(FutureWarning, match=msg):
+        om.minimize(
+            lambda x: x @ x,
+            np.arange(3),
+            algorithm="scipy_lbfgsb",
+            bounds=om.Bounds(lower=np.full(3, -1), upper=np.full(3, 2)),
+            multistart={"optimization_error_handling": "continue"},
+        )
+
+
+def test_multistart_option_exploration_error_handling_option_is_deprecated():
+    msg = "The exploration_error_handling option is deprecated and will be removed"
+    with pytest.warns(FutureWarning, match=msg):
+        om.minimize(
+            lambda x: x @ x,
+            np.arange(3),
+            algorithm="scipy_lbfgsb",
+            bounds=om.Bounds(lower=np.full(3, -1), upper=np.full(3, 2)),
+            multistart={"exploration_error_handling": "continue"},
+        )
diff --git a/tests/optimagic/test_typed_dicts_consistency.py b/tests/optimagic/test_typed_dicts_consistency.py
new file mode 100644
index 000000000..5b66b10b8
--- /dev/null
+++ b/tests/optimagic/test_typed_dicts_consistency.py
@@ -0,0 +1,41 @@
+from typing import get_args, get_type_hints
+
+from optimagic.optimization.multistart_options import (
+    MultistartOptions,
+    MultistartOptionsDict,
+)
+from optimagic.parameters.scaling import ScalingOptions, ScalingOptionsDict
+
+
+def assert_attributes_and_type_hints_are_equal(dataclass, typed_dict):
+    """Test that dataclass and typed_dict have same attributes and types.
+
+    This assertion purposefully ignores that all type hints in the typed dict are
+    wrapped by typing.NotRequired.
+
+    As there is no easy way to *not* read the NotRequired types in 3.10, we need to
+    activate include_extras=True to get the NotRequired types in Python 3.11 and
+    above. Once we drop support for Python 3.10, we can remove the
+    include_extras=True argument and the removal of the NotRequired types.
+
+    Args:
+        dataclass: An instance of a dataclass
+        typed_dict: An instance of a typed dict
+
+    """
+    types_from_dataclass = get_type_hints(dataclass)
+    types_from_typed_dict = get_type_hints(typed_dict, include_extras=True)
+    types_from_typed_dict = {
+        # Remove typing.NotRequired from the types
+        k: get_args(v)[0]
+        for k, v in types_from_typed_dict.items()
+    }
+    assert types_from_dataclass == types_from_typed_dict
+
+
+def test_scaling_options_and_dict_have_same_attributes():
+    assert_attributes_and_type_hints_are_equal(ScalingOptions, ScalingOptionsDict)
+
+
+def test_multistart_options_and_dict_have_same_attributes():
+    assert_attributes_and_type_hints_are_equal(MultistartOptions, MultistartOptionsDict)
diff --git a/tests/optimagic/visualization/test_history_plots.py b/tests/optimagic/visualization/test_history_plots.py
index bfe4b03f6..7f79befc5 100644
--- a/tests/optimagic/visualization/test_history_plots.py
+++ b/tests/optimagic/visualization/test_history_plots.py
@@ -1,6 +1,7 @@
 import itertools
 
 import numpy as np
+import optimagic as om
 import pytest
 from optimagic.optimization.optimize import minimize
 from optimagic.parameters.bounds import Bounds
@@ -19,11 +20,9 @@ def minimize_result():
                 params=np.arange(5),
                 algorithm=algorithm,
                 bounds=bounds,
-                multistart=multistart,
-                multistart_options={
-                    "n_samples": 1000,
-                    "convergence.max_discoveries": 5,
-                },
+                multistart=om.MultistartOptions(
+                    n_samples=1000, convergence_max_discoveries=5
+                ),
             )
             res.append(_res)
         out[multistart] = res
@@ -102,8 +101,7 @@ def test_criterion_plot_different_input_types():
         params=np.arange(5),
         algorithm="scipy_lbfgsb",
         bounds=bounds,
-        multistart=True,
-        multistart_options={"n_samples": 1000, "convergence.max_discoveries": 5},
+        multistart=om.MultistartOptions(n_samples=1000, convergence_max_discoveries=5),
         log_options={"fast_logging": True},
         logging="test.db",
     )
@@ -113,8 +111,7 @@ def test_criterion_plot_different_input_types():
         params=np.arange(5),
         algorithm="scipy_lbfgsb",
         bounds=bounds,
-        multistart=True,
-        multistart_options={"n_samples": 1000, "convergence.max_discoveries": 5},
+        multistart=om.MultistartOptions(n_samples=1000, convergence_max_discoveries=5),
     )
 
     results = ["test.db", res]