diff --git a/pypesto/visualize/_style.py b/pypesto/visualize/_style.py
new file mode 100644
index 000000000..4475f7e53
--- /dev/null
+++ b/pypesto/visualize/_style.py
@@ -0,0 +1,230 @@
+"""
+Visual style for ``pypesto.visualize``.
+
+Default constants, the ``style_kwargs`` registry, and small cross-module
+helpers.
+
+Users override any default per call via ``style_kwargs``, validated against
+:data:`_DEFAULTS`::
+
+    waterfall(result, style_kwargs={"mle_color": "tab:purple"})
+
+Keys are named after the **visual element**, not the plot using them:
+``line_*`` (any line), ``dash_*`` (any tick / rug marker), ``rectangle_*``
+(fills), ``bound_*`` (parameter bounds). Plot-specific keys (e.g.
+``trace_linewidth``) only when a default genuinely diverges.
+
+TODO (capstone): trim the naming note above once the series has settled.
+"""
+
+from __future__ import annotations
+
+import warnings
+from typing import Literal
+
+import matplotlib as mpl
+import matplotlib.axes
+import numpy as np
+from matplotlib.lines import Line2D
+
+# Colors — semantic roles
+# -----------------------
+MLE_COLOR = "#d62728"  # tab:red — best cluster + MLE markers
+OUTLIER_COLOR = "#b3b3b3"  # mid-grey — singleton / outlier starts
+
+# Colormaps
+# ---------
+CMAP_DISCRETE = "tab10"  # qualitative: cluster + per-variable colours
+
+# Lines (KDE curves, simulation / model-fit lines, …)
+# ---------------------------------------------------
+LINE_COLOR = "#145685"
+LINEWIDTH = 1.5
+
+# Dash markers (rug ticks, CI endpoints, …)
+# -----------------------------------------
+DASH_COLOR = "#174261"
+DASH_LINEWIDTH = 1.2  # markeredgewidth
+DASH_MARKERSIZE = 10  # marker length
+DASH_ALPHA = 0.8
+
+# Rectangle / histogram fills
+# ---------------------------
+RECTANGLE_COLOR = "#3182bd"
+RECTANGLE_EDGECOLOR = "#000000"
+RECTANGLE_LINEWIDTH = 1.0
+RECTANGLE_ALPHA = 0.6
+
+# Grid sizing — per-panel inches for multi-panel grids (size=None default)
+# ------------------------------------------------------------------------
+GRID_SIZE_PER_COL = 3.5
+GRID_SIZE_PER_ROW = 2.5
+
+# Parameter bounds
+# ----------------
+BOUND_LINESTYLE = "--"
+BOUND_COLOR = "0.5"
+BOUND_LINEWIDTH = 1.4
+BOUND_ALPHA = 0.95
+BOUND_VIEW_MARGIN = (
+    0.03  # axis-limit padding so bound lines aren't flush with the spine
+)
+
+# Style registry
+# --------------
+
+_DEFAULTS: dict[str, object] = {
+    "mle_color": MLE_COLOR,
+    "outlier_color": OUTLIER_COLOR,
+    "cmap_discrete": CMAP_DISCRETE,
+    "line_color": LINE_COLOR,
+    "linewidth": LINEWIDTH,
+    "dash_color": DASH_COLOR,
+    "dash_linewidth": DASH_LINEWIDTH,
+    "dash_markersize": DASH_MARKERSIZE,
+    "dash_alpha": DASH_ALPHA,
+    "rectangle_color": RECTANGLE_COLOR,
+    "rectangle_edgecolor": RECTANGLE_EDGECOLOR,
+    "rectangle_linewidth": RECTANGLE_LINEWIDTH,
+    "rectangle_alpha": RECTANGLE_ALPHA,
+    "bound_color": BOUND_COLOR,
+    "bound_linestyle": BOUND_LINESTYLE,
+    "bound_linewidth": BOUND_LINEWIDTH,
+    "bound_alpha": BOUND_ALPHA,
+}
+
+
+def resolve_style(style_kwargs: dict | None = None) -> dict:
+    """Return the effective style dict, merging defaults with caller overrides.
+
+    Parameters
+    ----------
+    style_kwargs:
+        User-supplied overrides. Unknown keys emit a ``UserWarning`` so
+        typos surface immediately.
+
+    Returns
+    -------
+    dict
+        Merged style dict with all keys from :data:`_DEFAULTS`, with
+        caller overrides applied on top.
+    """
+    style = dict(_DEFAULTS)
+    if style_kwargs:
+        unknown = set(style_kwargs) - set(_DEFAULTS)
+        if unknown:
+            warnings.warn(
+                f"Unknown style_kwargs keys: {sorted(unknown)}. "
+                f"Valid keys: {sorted(_DEFAULTS)}.",
+                UserWarning,
+                stacklevel=3,
+            )
+        style.update(style_kwargs)
+    return style
+
+
+# rcParams preset, not default, opt-in via ``apply_style()``
+# ---------------
+
+
+def apply_style() -> None:
+    """Apply pyPESTO's recommended matplotlib rcParams.
+
+    Sets larger axis/tick labels, removes top/right spines globally,
+    styles legends (auto-placed, framed, lightly translucent fill), and enables
+    ``constrained_layout`` for sensible panel spacing.
+
+    Opt-in: not called automatically. Users (and pyPESTO's example
+    notebooks/docs) call this once at the top of a session.
+    """
+    mpl.rcParams.update(
+        {
+            "axes.labelsize": 13,
+            "axes.labelweight": mpl.rcParamsDefault["axes.labelweight"],
+            "axes.titlesize": 14,
+            "axes.titleweight": "bold",
+            "xtick.labelsize": 11,
+            "ytick.labelsize": 11,
+            "legend.fontsize": mpl.rcParamsDefault["legend.fontsize"],
+            # Legends: auto-placed, framed, and lightly translucent so text reads
+            # clearly without making the legend feel heavy.
+            "legend.loc": "best",
+            "legend.frameon": True,
+            "legend.framealpha": 0.6,
+            "legend.edgecolor": "0.7",
+            "axes.spines.top": False,
+            "axes.spines.right": False,
+            "axes.grid": False,
+            "figure.constrained_layout.use": True,
+        }
+    )
+
+
+# Bound-line helpers
+# ------------------
+
+
+def _bounds_legend_handle(
+    label: str = "Bounds", style: dict | None = None
+) -> Line2D:
+    """Return a Line2D matching the bound style suitable as a legend handle."""
+    s = style or {}
+    return Line2D(
+        [0],
+        [0],
+        color=s.get("bound_color", BOUND_COLOR),
+        linestyle=s.get("bound_linestyle", BOUND_LINESTYLE),
+        linewidth=s.get("bound_linewidth", BOUND_LINEWIDTH),
+        alpha=s.get("bound_alpha", BOUND_ALPHA),
+        label=label,
+    )
+
+
+def draw_bounds_1d(
+    ax: matplotlib.axes.Axes,
+    lb: float,
+    ub: float,
+    *,
+    axis: Literal["x", "y"] = "x",
+    view_margin: bool = True,
+    style: dict | None = None,
+) -> Line2D:
+    """Draw the canonical pyPESTO parameter-bound lines on *ax*.
+
+    ``axis="x"`` draws two vertical dashed lines (``axvline``) at *lb* and
+    *ub*; ``axis="y"`` draws two horizontal dashed lines (``axhline``).
+
+    When *view_margin* is true the corresponding axis limits are extended by
+    :data:`BOUND_VIEW_MARGIN` * (ub - lb) so the bound lines are visible
+    rather than flush with the spine.
+
+    Returns a :class:`~matplotlib.lines.Line2D` that can be passed as a
+    legend handle (the lines drawn on the axis are not labeled to keep the
+    automatic legend clean).
+    """
+    if axis not in ("x", "y"):
+        raise ValueError(f"axis must be 'x' or 'y', got {axis!r}")
+    s = style or {}
+    color = s.get("bound_color", BOUND_COLOR)
+    linestyle = s.get("bound_linestyle", BOUND_LINESTYLE)
+    linewidth = s.get("bound_linewidth", BOUND_LINEWIDTH)
+    alpha = s.get("bound_alpha", BOUND_ALPHA)
+    drawer = ax.axvline if axis == "x" else ax.axhline
+    for bound in (lb, ub):
+        drawer(
+            bound,
+            color=color,
+            linestyle=linestyle,
+            linewidth=linewidth,
+            alpha=alpha,
+            zorder=1,
+        )
+    if view_margin and np.isfinite(lb) and np.isfinite(ub) and ub > lb:
+        margin = BOUND_VIEW_MARGIN * (ub - lb)
+        if axis == "x":
+            cur_lo, cur_hi = ax.get_xlim()
+            ax.set_xlim(min(cur_lo, lb - margin), max(cur_hi, ub + margin))
+        else:
+            cur_lo, cur_hi = ax.get_ylim()
+            ax.set_ylim(min(cur_lo, lb - margin), max(cur_hi, ub + margin))
+    return _bounds_legend_handle(style=s)
diff --git a/pypesto/visualize/clust_color.py b/pypesto/visualize/clust_color.py
index e73a5087b..98c5c5484 100644
--- a/pypesto/visualize/clust_color.py
+++ b/pypesto/visualize/clust_color.py
@@ -1,4 +1,7 @@
-import matplotlib.cm as cm
+from __future__ import annotations
+
+import matplotlib.colors as mcolors
+import matplotlib.pyplot as plt
 import numpy as np
 from matplotlib.colors import is_color_like
 
@@ -6,10 +9,56 @@
 
 # for typehints
 from ..C import COLOR
+from ._style import resolve_style
+
+
+def _build_cluster_palette(style: dict) -> np.ndarray:
+    """Sample non-best cluster colors from ``cmap_discrete``.
+
+    Colors close to ``mle_color`` or ``outlier_color`` are filtered out so
+    those reserved roles remain visually distinct from cycled cluster colors.
+    """
+    cmap = plt.get_cmap(style["cmap_discrete"])
+    reserved = [
+        np.array(mcolors.to_rgb(style["mle_color"])),
+        np.array(mcolors.to_rgb(style["outlier_color"])),
+    ]
+
+    # We remove colors from the cluster palette that are too close to the
+    # reserved colors (MLE red and outlier grey). The distance threshold is
+    # just a reasonable heuristic.
+    _RESERVED_COLOR_DISTANCE = 0.2
+
+    # Number of evenly-spaced samples taken from a continuous cmap (e.g. viridis)
+    # when ``cmap_discrete`` is set to one. Categorical cmaps (e.g. tab10) use
+    # all their listed colors and ignore this.
+    _CMAP_DISCRETE_SAMPLES = 10
+
+    if hasattr(cmap, "colors"):
+        candidates = [mcolors.to_rgba(c) for c in cmap.colors]
+    else:
+        candidates = [
+            cmap(i / (_CMAP_DISCRETE_SAMPLES - 1))
+            for i in range(_CMAP_DISCRETE_SAMPLES)
+        ]
+    palette = [
+        c
+        for c in candidates
+        if all(
+            np.linalg.norm(np.array(c[:3]) - r) > _RESERVED_COLOR_DISTANCE
+            for r in reserved
+        )
+    ]
+    if not palette:
+        palette = candidates
+    return np.array(palette)
 
 
 def assign_clustered_colors(
-    vals: np.ndarray, balance_alpha: bool = True, highlight_global: bool = True
+    vals: np.ndarray,
+    balance_alpha: bool = True,
+    highlight_global: bool = True,
+    style: dict | None = None,
 ):
     """
     Cluster and assign colors.
@@ -23,6 +72,10 @@ def assign_clustered_colors(
         avoid overplotting
     highlight_global:
         flag indicating whether global optimum should be highlighted
+    style:
+        Pre-resolved visualization style dict, as returned by
+        :func:`pypesto.visualize._style.resolve_style`. When ``None``, defaults
+        are used.
 
     Returns
     -------
@@ -36,6 +89,12 @@ def assign_clustered_colors(
     # assign clusters
     clusters, cluster_size = assign_clusters(vals)
 
+    if style is None:
+        style = resolve_style(None)
+    palette = _build_cluster_palette(style)
+    mle_rgba = list(mcolors.to_rgba(style["mle_color"]))
+    outlier_rgb = list(mcolors.to_rgb(style["outlier_color"]))
+
     # create list of colors, which has the correct shape
     n_clusters = 1 + max(clusters) - sum(cluster_size == 1)
 
@@ -43,13 +102,12 @@ def assign_clustered_colors(
     if highlight_global and cluster_size[0] > 1:
         n_clusters -= 1
 
-    # fill color array from colormap
-    colormap = cm.ScalarMappable().to_rgba
-    color_list = colormap(np.linspace(0.0, 1.0, n_clusters))
+    # fill color array by cycling through the categorical cluster palette
+    color_list = palette[np.arange(n_clusters) % len(palette)].copy()
 
-    # best optimum should be colored in red
+    # best optimum should be colored in MLE red
     if highlight_global and cluster_size[0] > 1:
-        color_list = np.concatenate(([[1.0, 0.0, 0.0, 1.0]], color_list))
+        color_list = np.concatenate(([mle_rgba], color_list))
 
     # We have clustered the results. However, clusters may have size 1,
     # so we need to rearrange the regroup the results into "no_clusters",
@@ -64,8 +122,8 @@ def assign_clustered_colors(
     if balance_alpha:
         # set minimal alpha value to avoid non-visible colors
         min_alpha = 0.01
-        # assign neutral color, add 1 for avoiding division by zero
-        grey = [0.7, 0.7, 0.7, min(1.0, 5.0 / (no_clusters.size + 1.0))]
+        # alpha shrinks with the number of singletons to avoid overplotting
+        grey = [*outlier_rgb, min(1.0, 5.0 / (no_clusters.size + 1.0))]
 
         # reduce alpha level depend on size of each cluster
         n_cluster_size = np.delete(cluster_size, no_clusters)
@@ -74,10 +132,9 @@ def assign_clustered_colors(
                 1.0, max(5.0 / n_cluster_size[icluster], min_alpha)
             )
     else:
-        # assign neutral color
-        grey = [0.7, 0.7, 0.7, 1.0]
+        grey = [*outlier_rgb, 1.0]
 
-    # create a color list, prfilled with grey values
+    # create a color list, prefilled with grey values
     colors = np.array([grey] * clusters.size)
 
     # assign colors to real clusters
@@ -86,9 +143,9 @@ def assign_clustered_colors(
         ind_of_iclust = np.argwhere(clusters == iclust).flatten()
         colors[ind_of_iclust, :] = color_list[icol, :]
 
-    # if best value was found only once: replace it with red
+    # if best value was found only once: replace it with MLE red
     if highlight_global and cluster_size[0] == 1:
-        colors[0] = [1.0, 0.0, 0.0, 1.0]
+        colors[0] = mle_rgba
 
     return colors
 
@@ -98,6 +155,7 @@ def assign_colors(
     colors: COLOR | list[COLOR] | np.ndarray | None = None,
     balance_alpha: bool = True,
     highlight_global: bool = True,
+    style: dict | None = None,
 ) -> np.ndarray:
     """
     Assign colors or format user specified colors.
@@ -113,6 +171,10 @@ def assign_colors(
         avoid overplotting
     highlight_global:
         flag indicating whether global optimum should be highlighted
+    style:
+        Pre-resolved visualization style dict, as returned by
+        :func:`pypesto.visualize._style.resolve_style`. When ``None``, defaults
+        are used.
 
     Returns
     -------
@@ -129,6 +191,7 @@ def assign_colors(
             vals,
             balance_alpha=balance_alpha,
             highlight_global=highlight_global,
+            style=style,
         )
 
     # Get number of elements and use user assigned colors
@@ -160,6 +223,7 @@ def assign_colors(
 def assign_colors_for_list(
     num_entries: int,
     colors: COLOR | list[COLOR] | np.ndarray | None = None,
+    style: dict | None = None,
 ) -> list[list[float]] | np.ndarray:
     """
     Create a list of colors for a list of items.
@@ -173,6 +237,10 @@ def assign_colors_for_list(
         number of results in list
     colors:
         list of colors, or single color
+    style:
+        Pre-resolved visualization style dict, as returned by
+        :func:`pypesto.visualize._style.resolve_style`. When ``None``, defaults
+        are used.
 
     Returns
     -------
@@ -188,7 +256,10 @@ def assign_colors_for_list(
 
         # we don't want alpha levels for all plotting routines in this case...
         colors = assign_colors(
-            dummy_clusters, balance_alpha=False, highlight_global=False
+            dummy_clusters,
+            balance_alpha=False,
+            highlight_global=False,
+            style=style,
         )
 
         # dummy cluster had twice as many entries as really there. Reduce.
@@ -201,4 +272,5 @@ def assign_colors_for_list(
         colors=colors,
         balance_alpha=False,
         highlight_global=False,
+        style=style,
     )
diff --git a/pypesto/visualize/misc.py b/pypesto/visualize/misc.py
index b74f2b55f..3eebf493f 100644
--- a/pypesto/visualize/misc.py
+++ b/pypesto/visualize/misc.py
@@ -26,6 +26,12 @@
 )
 from ..result import Result
 from ..util import assign_clusters, delete_nan_inf
+from ._style import (
+    GRID_SIZE_PER_COL,
+    GRID_SIZE_PER_ROW,
+    draw_bounds_1d,
+    resolve_style,
+)
 from .clust_color import assign_colors_for_list
 
 logger = logging.getLogger(__name__)
@@ -35,6 +41,7 @@ def process_result_list(
     results: Result | list[Result],
     colors: COLOR | list[COLOR] | np.ndarray | None = None,
     legends: str | list[str] | None = None,
+    style: dict | None = None,
 ) -> tuple[list[Result], list[COLOR], list[str]]:
     """
     Assign colors and legends to a list of results, check user provided lists.
@@ -47,6 +54,10 @@ def process_result_list(
         list of colors recognized by matplotlib, or single color
     legends:
         labels for line plots
+    style:
+        Pre-resolved visualization style dict, as returned by
+        :func:`pypesto.visualize._style.resolve_style`. When ``None``, defaults
+        are used.
 
     Returns
     -------
@@ -84,7 +95,7 @@ def process_result_list(
             legend_type_error = True
     else:
         # if more than one result is passed, we use one color per result
-        colors = assign_colors_for_list(len(results), colors)
+        colors = assign_colors_for_list(len(results), colors, style=style)
 
         # check whether list of legends has the correct length
         if legends is None:
@@ -486,7 +497,9 @@ def get_axes_array(
         Expected grid shape.
     size:
         Figure size ``(width, height)`` in inches; only used when ``axes``
-        is None.
+        is None. When ``None`` a single panel uses matplotlib's default
+        figure size, and a multi-panel grid uses
+        ``(GRID_SIZE_PER_COL * ncols, GRID_SIZE_PER_ROW * nrows)``.
 
     Returns
     -------
@@ -494,6 +507,8 @@ def get_axes_array(
         A 2-D NumPy object array containing matplotlib Axes.
     """
     if axes is None:
+        if size is None and nrows * ncols > 1:
+            size = (GRID_SIZE_PER_COL * ncols, GRID_SIZE_PER_ROW * nrows)
         _, axes = plt.subplots(
             nrows,
             ncols,
@@ -625,6 +640,140 @@ def plot_diagonal_marginal(
     ax.set_ylabel("Count")
 
 
+def plot_density_panel(
+    ax: matplotlib.axes.Axes,
+    values: np.ndarray,
+    bins: int | str = "auto",
+    bw_method: str = "scott",
+    style: dict | None = None,
+    *,
+    show_hist: bool = True,
+    show_kde: bool = True,
+    show_rug: bool = True,
+    show_bounds: bool = False,
+    lb: float | None = None,
+    ub: float | None = None,
+):
+    """Draw a density panel: histogram, KDE overlay, rug marks, and bounds.
+
+    Element styling is read from the resolved *style* dict:
+
+    - histogram bars: ``rectangle_*``
+    - KDE line: ``line_color`` / ``linewidth``
+    - rug marks: ``dash_color`` / ``dash_linewidth`` / ``dash_markersize`` /
+      ``dash_alpha``
+    - parameter-bound lines: ``bound_*`` (drawn only when ``show_bounds`` is
+      ``True`` and ``lb`` / ``ub`` are finite)
+
+    Sets x-axis limits to the data range with a 5% margin; if bounds are
+    drawn, the limits are extended to include them (never shrunk).
+
+    Parameters
+    ----------
+    ax:
+        Axes to draw into.
+    values:
+        1-D array of data values.
+    bins:
+        Histogram bins — passed directly to :func:`matplotlib.axes.Axes.hist`.
+    bw_method:
+        Bandwidth method for :class:`scipy.stats.gaussian_kde`.
+    style:
+        Pre-resolved visualization style dict, as returned by
+        :func:`pypesto.visualize._style.resolve_style`. When ``None``, defaults
+        are used.
+    show_hist:
+        Whether to draw the histogram bars.
+    show_kde:
+        Whether to draw the KDE line overlay.
+    show_rug:
+        Whether to draw rug marks along the x-axis.
+    show_bounds:
+        Whether to draw parameter-bound lines via
+        :func:`pypesto.visualize._style.draw_bounds_1d`. Requires ``lb`` and
+        ``ub`` to be passed and finite; silently skipped otherwise.
+    lb, ub:
+        Lower and upper parameter bounds. Used only when ``show_bounds`` is
+        ``True``.
+
+    Returns
+    -------
+    The bound legend handle if bounds were drawn, otherwise ``None``.
+    Callers wire this into their per-panel legend.
+    """
+    from scipy.stats import gaussian_kde
+
+    style = style if style is not None else resolve_style(None)
+
+    values = np.asarray(values, dtype=float)
+    values = values[np.isfinite(values)]
+    if values.size == 0:
+        return
+
+    if show_hist:
+        ax.hist(
+            values,
+            bins=bins,
+            density=True,
+            color=style["rectangle_color"],
+            alpha=style["rectangle_alpha"],
+            edgecolor=style["rectangle_edgecolor"],
+            linewidth=style["rectangle_linewidth"],
+        )
+
+    if show_kde and len(values) > 1 and np.std(values) > 0:
+        try:
+            kde = gaussian_kde(values, bw_method=bw_method)
+            # Extend 3 bandwidths beyond the data so the curve tapers to zero.
+            bw = kde.factor * np.std(values, ddof=1)
+            x_grid = np.linspace(
+                values.min() - 3 * bw, values.max() + 3 * bw, 300
+            )
+            ax.plot(
+                x_grid,
+                kde(x_grid),
+                color=style["line_color"],
+                linewidth=style["linewidth"],
+            )
+        except np.linalg.LinAlgError:
+            pass
+
+    if show_rug:
+        ax.plot(
+            values,
+            np.zeros(len(values)),
+            marker="|",
+            linestyle="none",
+            color=style["dash_color"],
+            alpha=style["dash_alpha"],
+            markersize=style["dash_markersize"],
+            markeredgewidth=style["dash_linewidth"],
+            transform=ax.get_xaxis_transform(),
+            clip_on=False,
+            zorder=5,
+        )
+
+    # Frame the panel tightly to the data (5% margin around the range, or a
+    # sensible fallback when the data is constant).
+    data_min = float(values.min())
+    data_max = float(values.max())
+    spread = data_max - data_min
+    margin = spread * 0.05 if spread > 0 else max(abs(data_min) * 0.05, 1.0)
+    ax.set_xlim(data_min - margin, data_max + margin)
+
+    # Draw parameter bounds (xlim is extended by draw_bounds_1d when bounds
+    # fall outside the data frame; never shrunk).
+    if (
+        show_bounds
+        and lb is not None
+        and ub is not None
+        and np.isfinite(lb)
+        and np.isfinite(ub)
+    ):
+        return draw_bounds_1d(ax, lb, ub, axis="x", style=style)
+    return None
+
+
 #: Sentinel meaning "this kwarg was not passed at all."
 #: Use as the default for deprecated kwargs so that an explicit
 #: ``f(old_kwarg=None)`` can be detected and warned about.
@@ -632,29 +781,38 @@ def plot_diagonal_marginal(
 
 
 def process_deprecated_kwarg(
-    canonical_name: str,
+    canonical_name: str | None,
     canonical_value,
     deprecated_name: str,
     deprecated_value=_UNSET,
     stacklevel: int = 3,
+    note: str | None = None,
 ):
     """
-    Resolve a kwarg that has been renamed.
+    Resolve a kwarg that has been renamed or removed.
 
     The deprecated kwarg must use :data:`_UNSET` as its default in the
     calling function so that an explicit ``f(old_kwarg=None)`` is correctly
     detected and warned about.
 
-    Returns the canonical value if the deprecated kwarg was not passed,
-    the deprecated value (with a ``DeprecationWarning``) if only the old
-    name was used, or raises ``ValueError`` if both are given.
+    Two modes:
+
+    - **Rename** (``canonical_name`` is given): returns the canonical value
+      if the deprecated kwarg was not passed, the deprecated value (with a
+      ``DeprecationWarning``) if only the old name was used, or raises
+      ``ValueError`` if both are given.
+    - **Removal** (``canonical_name`` is ``None``): emits a
+      ``DeprecationWarning`` if the deprecated kwarg was passed and returns
+      ``None``. ``canonical_value`` is ignored.
 
     Parameters
     ----------
     canonical_name:
-        Name of the canonical (new) kwarg, used in messages.
+        Name of the canonical (new) kwarg, used in messages. Pass ``None``
+        to indicate the kwarg is removed without replacement.
     canonical_value:
-        Value passed under the canonical name (or ``None``).
+        Value passed under the canonical name (or ``None``). Ignored when
+        ``canonical_name`` is ``None``.
     deprecated_name:
         Name of the deprecated (old) kwarg, used in messages.
     deprecated_value:
@@ -663,12 +821,26 @@ def process_deprecated_kwarg(
         Forwarded to :func:`warnings.warn`. Default 3 attributes the
         warning to the caller of the public function that invoked this
         helper.
+    note:
+        Optional additional sentence appended to the deprecation message
+        (e.g. explaining where the behaviour moved to). Useful for the
+        removal case.
 
     Returns
     -------
     value:
-        The resolved value, or ``None`` if neither was given.
+        The resolved value, or ``None`` if neither was given (or in the
+        removal case).
     """
+    if canonical_name is None:
+        if deprecated_value is _UNSET:
+            return None
+        message = f"`{deprecated_name}` is deprecated and has no effect."
+        if note:
+            message = f"{message} {note}"
+        warnings.warn(message, DeprecationWarning, stacklevel=stacklevel)
+        return None
+
     if deprecated_value is _UNSET:
         return canonical_value
     if canonical_value is not None:
@@ -676,9 +848,10 @@ def process_deprecated_kwarg(
             f"Pass either `{canonical_name}` or the deprecated "
             f"`{deprecated_name}`, not both."
         )
-    warnings.warn(
-        f"`{deprecated_name}` is deprecated; use `{canonical_name}` instead.",
-        DeprecationWarning,
-        stacklevel=stacklevel,
+    message = (
+        f"`{deprecated_name}` is deprecated; use `{canonical_name}` instead."
     )
+    if note:
+        message = f"{message} {note}"
+    warnings.warn(message, DeprecationWarning, stacklevel=stacklevel)
     return deprecated_value
diff --git a/pypesto/visualize/optimization_stats.py b/pypesto/visualize/optimization_stats.py
index 5669627f8..e98e82e0a 100644
--- a/pypesto/visualize/optimization_stats.py
+++ b/pypesto/visualize/optimization_stats.py
@@ -8,6 +8,7 @@
 
 from ..C import COLOR
 from ..result import Result
+from ._style import resolve_style
 from .clust_color import assign_colors, assign_colors_for_list
 from .misc import (
     get_ax,
@@ -28,6 +29,7 @@ def optimization_run_properties_one_plot(
     legends: str | list[str] | None = None,
     plot_type: str = "line",
     ax: matplotlib.axes.Axes | None = None,
+    style_kwargs: dict | None = None,
 ) -> matplotlib.axes.Axes:
     """
     Plot stats for allproperties specified in properties_to_plot on one plot.
@@ -54,6 +56,16 @@ def optimization_run_properties_one_plot(
         Labels, one label per optimization property
     plot_type:
         Specifies plot type. Possible values: 'line' and 'hist'
+    style_kwargs:
+        Style overrides. Keys used by this function:
+
+        - ``cmap_discrete`` — the categorical palette from which
+          per-property line colours are sampled. Only consulted when
+          ``colors`` is ``None``; an explicit ``colors`` short-circuits
+          palette selection.
+
+        All valid keys and their defaults are listed in
+        :data:`pypesto.visualize._style._DEFAULTS`.
 
     Returns
     -------
@@ -76,6 +88,8 @@ def optimization_run_properties_one_plot(
             colors=[[.5, .9, .9, .3], [.2, .1, .9, .5]]
         )
     """
+    style = resolve_style(style_kwargs)
+
     if properties_to_plot is None:
         properties_to_plot = [
             "time",
@@ -87,7 +101,7 @@ def optimization_run_properties_one_plot(
         ]
 
     if colors is None:
-        colors = assign_colors_for_list(len(properties_to_plot))
+        colors = assign_colors_for_list(len(properties_to_plot), style=style)
     elif is_color_like(colors):
         colors = [colors]
 
@@ -136,6 +150,7 @@ def optimization_run_properties_per_multistart(
     legends: str | list[str] | None = None,
     plot_type: str = "line",
     axes: np.ndarray | None = None,
+    style_kwargs: dict | None = None,
 ) -> np.ndarray:
     """
     One plot per optimization property in properties_to_plot.
@@ -160,6 +175,19 @@ def optimization_run_properties_per_multistart(
         Labels for line plots, one label per result object
     plot_type:
         Specifies plot type. Possible values: 'line' and 'hist'
+    style_kwargs:
+        Style overrides forwarded to
+        :func:`optimization_run_property_per_multistart`. Keys used by
+        this function:
+
+        - ``cmap_discrete``, ``mle_color``, ``outlier_color`` — colours
+          of the per-start scatter when clustering is applied (best
+          cluster, secondary clusters, isolated starts respectively).
+          Only consulted when ``colors`` is ``None``; an explicit
+          ``colors`` short-circuits clustering.
+
+        All valid keys and their defaults are listed in
+        :data:`pypesto.visualize._style._DEFAULTS`.
 
     Returns
     -------
@@ -223,6 +251,7 @@ def optimization_run_properties_per_multistart(
             colors=colors,
             legends=legends,
             plot_type=plot_type,
+            style_kwargs=style_kwargs,
         )
     return axes
 
@@ -236,6 +265,7 @@ def optimization_run_property_per_multistart(
     colors: COLOR | list[COLOR] | np.ndarray | None = None,
     legends: str | list[str] | None = None,
     plot_type: str = "line",
+    style_kwargs: dict | None = None,
 ) -> np.ndarray:
     """
     Plot stats for an optimization run property specified by opt_run_property.
@@ -268,12 +298,25 @@ def optimization_run_property_per_multistart(
         Labels for line plots, one label per result object
     plot_type:
         Specifies plot type. Possible values: 'line', 'hist', 'both'
+    style_kwargs:
+        Style overrides. Keys used by this function:
+
+        - ``cmap_discrete``, ``mle_color``, ``outlier_color`` — colours
+          of the per-start scatter when clustering is applied
+          (single-result default; best cluster, secondary clusters,
+          isolated starts respectively). Only consulted when
+          ``colors`` is ``None``; an explicit ``colors`` short-circuits
+          clustering.
+
+        All valid keys and their defaults are listed in
+        :data:`pypesto.visualize._style._DEFAULTS`.
 
     Returns
     -------
     axes:
         2-D NumPy array containing one matplotlib Axes per panel.
     """
+    style = resolve_style(style_kwargs)
     supported_properties = {
         "time": "Wall-clock time (seconds)",
         "n_fval": "Number of function evaluations",
@@ -291,7 +334,9 @@ def optimization_run_property_per_multistart(
         )
 
     # parse input
-    (results, colors, legends) = process_result_list(results, colors, legends)
+    (results, colors, legends) = process_result_list(
+        results, colors, legends, style=style
+    )
 
     ncols = 2 if plot_type == "both" else 1
     axes = get_axes_array(axes=axes, nrows=1, ncols=ncols, size=size)
@@ -320,6 +365,7 @@ def optimization_run_property_per_multistart(
                 start_indices,
                 colors[j],
                 legends[j],
+                style=style,
             )
 
             stats_lowlevel(
@@ -331,6 +377,7 @@ def optimization_run_property_per_multistart(
                 colors[j],
                 legends[j],
                 plot_type="hist",
+                style=style,
             )
         else:
             stats_lowlevel(
@@ -342,6 +389,7 @@ def optimization_run_property_per_multistart(
                 colors[j],
                 legends[j],
                 plot_type,
+                style=style,
             )
 
     if sum(legend is not None for legend in legends) > 0:
@@ -363,6 +411,7 @@ def stats_lowlevel(
     color: COLOR | list[COLOR] | np.ndarray | None = "C0",
     legend: str | None = None,
     plot_type: str = "line",
+    style: dict | None = None,
 ):
     """
     Plot values of the optimization run property across different multistarts.
@@ -389,6 +438,10 @@ def stats_lowlevel(
         Label describing the result
     plot_type:
         Specifies plot type. Possible values: 'line' and 'hist'
+    style:
+        Pre-resolved visualization style dict, as returned by
+        :func:`pypesto.visualize._style.resolve_style`. When ``None``, defaults
+        are used.
 
     Returns
     -------
@@ -407,7 +460,9 @@ def stats_lowlevel(
     n_starts = len(values)
 
     # assign colors
-    colors = assign_colors(vals=fvals, colors=color, balance_alpha=False)
+    colors = assign_colors(
+        vals=fvals, colors=color, balance_alpha=False, style=style
+    )
 
     sorted_indices = sorted(range(n_starts), key=lambda j: fvals[j])
     values = values[sorted_indices]
diff --git a/pypesto/visualize/optimizer_history.py b/pypesto/visualize/optimizer_history.py
index 6fc624010..8baa6ad3f 100644
--- a/pypesto/visualize/optimizer_history.py
+++ b/pypesto/visualize/optimizer_history.py
@@ -15,6 +15,7 @@
 )
 from ..history import HistoryBase
 from ..result import Result
+from ._style import resolve_style
 from .clust_color import assign_colors
 from .misc import (
     get_ax,
@@ -44,6 +45,7 @@ def optimizer_history(
     | list[dict]
     | None = None,
     legends: str | list[str] | None = None,
+    style_kwargs: dict | None = None,
 ) -> matplotlib.axes.Axes:
     """
     Plot history of optimizer.
@@ -86,17 +88,32 @@ def optimizer_history(
         least a function value fval
     legends:
         Labels for line plots, one label per result object
+    style_kwargs:
+        Style overrides. Keys used by this function:
+
+        - ``cmap_discrete``, ``mle_color``, ``outlier_color`` — colours
+          of the per-start history traces when clustering is applied
+          (best cluster, secondary clusters, isolated starts respectively).
+          Only consulted when ``colors`` is ``None``; an explicit
+          ``colors`` short-circuits clustering.
+
+        All valid keys and their defaults are listed in
+        :data:`pypesto.visualize._style._DEFAULTS`.
 
     Returns
     -------
     ax:
         The plot axes.
     """
+    style = resolve_style(style_kwargs)
+
     if isinstance(start_indices, int):
         start_indices = list(range(start_indices))
 
     # parse input
-    (results, colors, legends) = process_result_list(results, colors, legends)
+    (results, colors, legends) = process_result_list(
+        results, colors, legends, style=style
+    )
 
     for j, result in enumerate(results):
         # extract cost function values from result
@@ -119,6 +136,7 @@ def optimizer_history(
             x_label=x_label,
             y_label=y_label,
             legend_text=legends[j],
+            style=style,
         )
 
     # parse and apply plotting options
@@ -139,6 +157,7 @@ def optimizer_history_lowlevel(
     x_label: str = "Optimizer steps",
     y_label: str = "Objective value",
     legend_text: str | None = None,
+    style: dict | None = None,
 ) -> matplotlib.axes.Axes:
     """
     Plot optimizer history using list of numpy arrays.
@@ -162,6 +181,10 @@ def optimizer_history_lowlevel(
         label for y-axis
     legend_text:
         Label for line plots
+    style:
+        Pre-resolved visualization style dict, as returned by
+        :func:`pypesto.visualize._style.resolve_style`. When ``None``, defaults
+        are used.
 
     Returns
     -------
@@ -194,7 +217,7 @@ def optimizer_history_lowlevel(
     # assign colors
     # note: this has to happen before sorting
     # to get the same colors in different plots
-    colors = assign_colors(fvals, colors)
+    colors = assign_colors(fvals, colors, style=style)
 
     # sort
     indices = sorted(range(n_fvals), key=lambda j: fvals[j])
diff --git a/pypesto/visualize/parameters.py b/pypesto/visualize/parameters.py
index 597f615f0..9c1d7e88e 100644
--- a/pypesto/visualize/parameters.py
+++ b/pypesto/visualize/parameters.py
@@ -5,6 +5,7 @@
 import numpy as np
 import pandas as pd
 from matplotlib.colors import Colormap
+from matplotlib.lines import Line2D
 from matplotlib.ticker import MaxNLocator
 
 from pypesto.util import delete_nan_inf
@@ -17,11 +18,15 @@
     InnerParameterType,
 )
 from ..result import Result
+from ._style import resolve_style
 from .clust_color import assign_colors
 from .misc import (
+    _UNSET,
     get_ax,
     get_axes_array,
+    plot_density_panel,
     plot_diagonal_marginal,
+    process_deprecated_kwarg,
     process_parameter_indices,
     process_result_list,
     process_start_indices,
@@ -53,6 +58,7 @@ def parameters(
     scale_to_interval: tuple[float, float] | None = None,
     plot_inner_parameters: bool = True,
     log10_scale_hier_sigma: bool = True,
+    style_kwargs: dict | None = None,
 ) -> matplotlib.axes.Axes:
     """
     Plot parameter values.
@@ -95,14 +101,29 @@ def parameters(
     log10_scale_hier_sigma:
         Flag indicating whether to scale inner parameters of type
         ``InnerParameterType.SIGMA`` to log10 (default: True).
+    style_kwargs:
+        Style overrides. Keys used by this function:
+
+        - ``cmap_discrete``, ``mle_color``, ``outlier_color`` — colours
+          of the per-start parameter traces when clustering is used
+          (best cluster, secondary clusters, isolated starts respectively).
+          Only consulted when ``colors`` is ``None``; an explicit
+          ``colors`` short-circuits clustering.
+
+        All valid keys and their defaults are listed in
+        :data:`pypesto.visualize._style._DEFAULTS`.
 
     Returns
     -------
     ax:
         The plot axes.
     """
+    style = resolve_style(style_kwargs)
+
     # parse input
-    (results, colors, legends) = process_result_list(results, colors, legends)
+    (results, colors, legends) = process_result_list(
+        results, colors, legends, style=style
+    )
 
     if isinstance(parameter_indices, str):
         if parameter_indices == "all":
@@ -161,6 +182,7 @@ def scale_parameters(x):
             colors=colors[j],
             legend_text=legends[j],
             balance_alpha=balance_alpha,
+            style=style,
         )
 
     # parse and apply plotting options
@@ -197,39 +219,78 @@ def scale_parameters(x):
 def parameter_hist(
     result: Result,
     parameter_name: str,
+    start_indices: int | list[int] | None = None,
+    plot_type: str = "both",
     bins: int | str = "auto",
+    bw_method: str = "scott",
+    show_bounds: bool = True,
+    title: str | None = "Parameter histogram",
+    size: tuple[float, float] | None = None,
     ax: matplotlib.axes.Axes | None = None,
-    size: tuple[float, float] | None = (18.5, 10.5),
-    color: COLOR | None = None,
-    start_indices: int | list[int] | None = None,
+    style_kwargs: dict | None = None,
+    color: COLOR = _UNSET,
 ) -> matplotlib.axes.Axes:
     """
-    Plot parameter values as a histogram.
+    Plot one parameter's values across starts as a histogram + KDE + rug.
 
     Parameters
     ----------
     result:
-        Optimization result obtained by 'optimize.py'
+        Optimization result obtained by 'optimize.py'.
     parameter_name:
-        The name of the parameter that should be plotted
+        Name of the parameter to plot.
+    start_indices:
+        Which optimization starts to include: a list of indices, or an int
+        ``n`` for the first ``n`` starts. Default: all starts.
+    plot_type: {'hist'|'kde'|'both'}
+        Histogram only, KDE line only, or both with rug marks (default).
     bins:
-        Specifies bins of the histogram
-    ax:
-        Axes object to use
+        Number of bins, or a matplotlib binning strategy (``'auto'``,
+        ``'sturges'``, …). Passed to ``ax.hist``.
+    bw_method: {'scott', 'silverman' | scalar | pair of scalars}
+        Kernel bandwidth method for the KDE overlay.
+    show_bounds:
+        If ``True`` (default) draw the parameter bound lines and frame the
+        x-axis to include them; if ``False`` frame tightly to the data.
+    title:
+        Axes title. Pass ``None`` to suppress.
     size:
-        Figure size (width, height) in inches. Is only applied when no ax
-        object is specified
+        Figure size in inches. Defaults to matplotlib's default.
+    ax:
+        Axes object to use.
+    style_kwargs:
+        Style overrides. Keys used by this function:
+
+        - ``rectangle_color``, ``rectangle_alpha``, ``rectangle_edgecolor``,
+          ``rectangle_linewidth`` — histogram bar styling.
+        - ``line_color``, ``linewidth`` — KDE curve styling.
+        - ``dash_color``, ``dash_linewidth``, ``dash_markersize``,
+          ``dash_alpha`` — rug-mark styling.
+        - ``bound_color``, ``bound_linestyle``, ``bound_linewidth``,
+          ``bound_alpha`` — parameter-bound line styling.
+
+        All valid keys and their defaults are listed in
+        :data:`pypesto.visualize._style._DEFAULTS`.
     color:
-        Color recognized by matplotlib.
-    start_indices:
-        List of integers specifying the multistarts to be plotted or
-        int specifying up to which start index should be plotted
+        Deprecated. Pass ``style_kwargs`` instead — see
+        ``rectangle_color`` / ``line_color`` / ``dash_color`` above.
 
     Returns
     -------
     ax:
         The plot axes.
     """
+    process_deprecated_kwarg(
+        canonical_name=None,
+        canonical_value=None,
+        deprecated_name="color",
+        deprecated_value=color,
+        note=(
+            "Pass style_kwargs={'rectangle_color': ..., 'line_color': ..., "
+            "'dash_color': ...} instead."
+        ),
+    )
+    style = resolve_style(style_kwargs)
     ax = get_ax(ax, size)
 
     xs = result.optimize_result.x
@@ -241,12 +302,68 @@ def parameter_hist(
         xs = [xs[ind] for ind in start_indices]
 
     parameter_index = result.problem.x_names.index(parameter_name)
-    parameter_values = [x[parameter_index] for x in xs]
+    parameter_values = np.array([x[parameter_index] for x in xs])
+
+    # bounds and scale for this parameter
+    lb_val = result.problem.lb_full[parameter_index]
+    ub_val = result.problem.ub_full[parameter_index]
+    x_scales = getattr(result.problem, "x_scales", None)
+    scale = x_scales[parameter_index] if x_scales is not None else None
+
+    bound_handle = plot_density_panel(
+        ax,
+        parameter_values,
+        bins=bins,
+        bw_method=bw_method,
+        style=style,
+        show_hist=(plot_type in ("hist", "both")),
+        show_kde=(plot_type in ("kde", "both")),
+        show_rug=(plot_type in ("hist", "both")),
+        show_bounds=show_bounds,
+        lb=lb_val,
+        ub=ub_val,
+    )
 
-    ax.hist(parameter_values, color=color, bins=bins, label=parameter_name)
-    ax.set_xlabel(parameter_name)
-    ax.set_ylabel("counts")
-    ax.set_title(f"{parameter_name}")
+    legend_handles, legend_labels = [], []
+    show_kde = plot_type in ("kde", "both")
+    show_rug = plot_type in ("hist", "both")
+    finite_vals = parameter_values[np.isfinite(parameter_values)]
+    if finite_vals.size > 0:
+        if show_kde:
+            legend_handles.append(
+                Line2D(
+                    [0], [0], color=style["line_color"], lw=style["linewidth"]
+                )
+            )
+            legend_labels.append("KDE")
+        if show_rug:
+            legend_handles.append(
+                Line2D(
+                    [0],
+                    [0],
+                    color=style["dash_color"],
+                    marker="|",
+                    lw=0,
+                    markersize=style["dash_markersize"],
+                    markeredgewidth=style["dash_linewidth"],
+                )
+            )
+            legend_labels.append("Starts")
+
+    if bound_handle is not None:
+        legend_handles.append(bound_handle)
+        legend_labels.append("Bounds")
+
+    if legend_handles:
+        ax.legend(handles=legend_handles, labels=legend_labels)
+
+    xlabel = (
+        f"{parameter_name} ({scale})" if scale is not None else parameter_name
+    )
+    if title is not None:
+        ax.set_title(title)
+    ax.set_xlabel(xlabel)
+    ax.set_ylabel("Density")
 
     return ax
 
@@ -264,6 +381,7 @@ def parameters_lowlevel(
     linestyle: str = "-",
     legend_text: str | None = None,
     balance_alpha: bool = True,
+    style: dict | None = None,
 ) -> matplotlib.axes.Axes:
     """
     Plot parameters plot using list of parameters.
@@ -292,13 +410,16 @@ def parameters_lowlevel(
     balance_alpha:
         Flag indicating whether alpha for large clusters should be reduced to
         avoid overplotting (default: True)
+    style:
+        Pre-resolved visualization style dict, as returned by
+        :func:`pypesto.visualize._style.resolve_style`. When ``None``, defaults
+        are used.
 
     Returns
     -------
     ax:
         The plot axes.
     """
-
     if size is None:
         # 0.5 inch height per parameter
         size = (18.5, max(xs.shape[1], 1) / 2)
@@ -307,7 +428,7 @@ def parameters_lowlevel(
 
     # assign colors
     colors = assign_colors(
-        vals=fvals, colors=colors, balance_alpha=balance_alpha
+        vals=fvals, colors=colors, balance_alpha=balance_alpha, style=style
     )
 
     # parameter indices
diff --git a/pypesto/visualize/profiles.py b/pypesto/visualize/profiles.py
index d61f5bd7e..c5fef9512 100644
--- a/pypesto/visualize/profiles.py
+++ b/pypesto/visualize/profiles.py
@@ -12,6 +12,7 @@
 from ..problem import Problem
 from ..profile import chi2_quantile_to_ratio
 from ..result import Result
+from ._style import resolve_style
 from .clust_color import assign_colors
 from .misc import get_ax, process_result_list
 from .reference_points import ReferencePoint, create_references
@@ -123,6 +124,7 @@ def profiles(
     show_bounds: bool = False,
     plot_objective_values: bool = False,
     quality_colors: bool = False,
+    style_kwargs: dict | None = None,
 ) -> matplotlib.axes.Axes:
     """
     Plot classical 1D profile plot.
@@ -173,12 +175,24 @@ def profiles(
         had to resample the parameter vector due to optimization failure of the previous two.
         Black indicates a step for which none of the above was necessary. This option is only
         available if there is only one result and one profile_list_id (one profile per plot).
+    style_kwargs:
+        Style overrides. Keys used by this function:
+
+        - ``cmap_discrete``, ``mle_color``, ``outlier_color`` — colours
+          of the per-result / per-profile-list profile lines (best
+          cluster, secondary clusters, isolated starts respectively).
+          Only consulted when ``colors`` is ``None``; an explicit
+          ``colors`` short-circuits clustering.
+
+        All valid keys and their defaults are listed in
+        :data:`pypesto.visualize._style._DEFAULTS`.
 
     Returns
     -------
     ax:
         The plot axes.
     """
+    style = resolve_style(style_kwargs)
 
     if colors is not None and quality_colors:
         raise ValueError(
@@ -195,7 +209,7 @@ def profiles(
 
     # parse input
     results, profile_list_ids, colors, legends = process_result_list_profiles(
-        results, profile_list_ids, legends, colors
+        results, profile_list_ids, legends, colors, style=style
     )
 
     # get the parameter ids to be plotted
@@ -621,6 +635,7 @@ def process_result_list_profiles(
     profile_list_ids: int | Sequence[int] | None,
     legends: str | list[str],
     colors: COLOR | list[COLOR] | np.ndarray | None = None,  # todo: check
+    style: dict | None = None,
 ) -> tuple[list[Result], list[int] | Sequence[int], list, list[str]]:
     """
     Assign colors and legends to a list of results.
@@ -637,6 +652,10 @@ def process_result_list_profiles(
         list of colors for plotting.
     legends:
         Legends for plotting
+    style:
+        Pre-resolved visualization style dict, as returned by
+        :func:`pypesto.visualize._style.resolve_style`. When ``None``, defaults
+        are used.
 
     Returns
     -------
@@ -653,7 +672,7 @@ def process_result_list_profiles(
         if len(results) != 1:
             # if we have no single result, then use the standard api
             results, colors, legends = process_result_list(
-                results, colors, legends
+                results, colors, legends, style=style
             )
             return results, profile_list_ids, colors, legends
     else:
@@ -662,7 +681,9 @@ def process_result_list_profiles(
 
     # If we have a single result, we may still have multiple profile_list_ids
     # which should be plotted separately: use profile_list_ids as results dummy
-    _, colors, legends = process_result_list(profile_list_ids, colors, legends)
+    _, colors, legends = process_result_list(
+        profile_list_ids, colors, legends, style=style
+    )
 
     return results, profile_list_ids, colors, legends
 
diff --git a/pypesto/visualize/sampling.py b/pypesto/visualize/sampling.py
index bb520d352..fb7ab3c37 100644
--- a/pypesto/visualize/sampling.py
+++ b/pypesto/visualize/sampling.py
@@ -25,12 +25,14 @@
 from ..ensemble import EnsemblePrediction, get_percentile_label
 from ..result import McmcPtResult, PredictionResult, Result
 from ..sample import calculate_ci_mcmc_sample
+from ._style import resolve_style
 from .misc import (
     _UNSET,
     get_ax,
     get_axes_array,
     hide_unused_axes,
     make_grid_shape,
+    plot_density_panel,
     plot_diagonal_marginal,
     process_deprecated_kwarg,
     rgba2rgb,
@@ -1308,48 +1310,78 @@ def sampling_scatter(
 def sampling_1d_marginals(
     result: Result,
     i_chain: int = 0,
-    parameter_indices: Sequence[int] = None,
+    parameter_indices: Sequence[int] | None = None,
     stepsize: int = 1,
     plot_type: str = "both",
+    bins: int | str = "auto",
     bw_method: str = "scott",
-    suptitle: str | None = None,
+    show_bounds: bool = True,
+    title: str | None = None,
     size: tuple[float, float] | None = None,
     axes: np.ndarray | None = None,
+    style_kwargs: dict | None = None,
     par_indices: Sequence[int] = _UNSET,
+    suptitle: str | None = _UNSET,
 ) -> np.ndarray:
     """
-    Plot marginals.
+    Plot 1-D marginals of the sampled parameters as histogram + KDE + rug.
 
     Parameters
     ----------
     result:
         The pyPESTO result object with filled sample result.
     i_chain:
-        Which chain to plot. Default: First chain.
-    parameter_indices: list of integer values
-        List of integer values specifying which parameters to plot.
-        Default: All parameters are shown.
+        Which chain to plot. Default: first chain.
+    parameter_indices:
+        Which parameters to plot, as a list of indices. Default: all parameters.
     stepsize:
-        Only one in `stepsize` values is plotted.
+        Thinning factor — plot every ``stepsize``-th sample (``1`` = all).
+        Reduces overplotting and speeds up rendering for long chains.
     plot_type: {'hist'|'kde'|'both'}
-        Specify whether to plot a histogram ('hist'), a kernel density estimate
-        ('kde'), or both ('both').
+        Histogram only, KDE line only, or both with rug marks (default).
+    bins:
+        Number of bins, or a matplotlib binning strategy (``'auto'``,
+        ``'sturges'``, …). Passed to ``ax.hist``.
     bw_method: {'scott', 'silverman' | scalar | pair of scalars}
-        Kernel bandwidth method.
-    suptitle:
-        Figure super title.
+        Kernel bandwidth method for the KDE overlay.
+    show_bounds:
+        If ``True`` (default) draw the parameter bound lines and frame each
+        panel's x-axis to include them; if ``False`` frame each panel tightly
+        to its data.
+    title:
+        Figure title. Default: none (grids omit a title by default).
     size:
-        Figure size in inches.
+        Figure size in inches. When ``None`` the grid uses
+        ``GRID_SIZE_PER_COL * num_col`` × ``GRID_SIZE_PER_ROW * num_row``
+        (defaults from :mod:`pypesto.visualize._style`).
     axes:
         Axes grid to use. Must match the computed subplot layout.
+    style_kwargs:
+        Style overrides. Keys used by this function:
+
+        - ``rectangle_color``, ``rectangle_alpha``, ``rectangle_edgecolor``,
+          ``rectangle_linewidth`` — histogram bar styling.
+        - ``line_color``, ``linewidth`` — KDE curve styling.
+        - ``dash_color``, ``dash_linewidth``, ``dash_markersize``,
+          ``dash_alpha`` — rug-mark styling.
+        - ``bound_color``, ``bound_linestyle``, ``bound_linewidth``,
+          ``bound_alpha`` — parameter-bound line styling.
+
+        All valid keys and their defaults are listed in
+        :data:`pypesto.visualize._style._DEFAULTS`.
     par_indices:
         Deprecated. Use ``parameter_indices`` instead.
+    suptitle:
+        Deprecated. Use ``title`` instead.
 
-    Return
-    --------
+    Returns
+    -------
     axes:
         2-D NumPy array containing one matplotlib Axes per panel.
     """
+    style = resolve_style(style_kwargs)
+    title = process_deprecated_kwarg("title", title, "suptitle", suptitle)
+
     parameter_indices = process_deprecated_kwarg(
         "parameter_indices",
         parameter_indices,
@@ -1357,8 +1389,6 @@ def sampling_1d_marginals(
         par_indices,
     )
 
-    import seaborn as sns
-
     # get data which should be plotted
     nr_params, params_fval, theta_lb, theta_ub, param_names = get_data_to_plot(
         result=result,
@@ -1368,43 +1398,93 @@ def sampling_1d_marginals(
     )
 
     num_row, num_col = make_grid_shape(nr_params)
-    if size is None and axes is None:
-        size = (3.5 * num_col, 2.5 * num_row)
     axes = get_axes_array(axes=axes, nrows=num_row, ncols=num_col, size=size)
     fig = axes.flat[0].figure
     axes = hide_unused_axes(axes=axes, n_used=nr_params, clear=True)
 
-    par_ax = dict(zip(param_names, axes.flat, strict=True))
+    par_ax = dict(zip(param_names, axes.flat[:nr_params], strict=True))
+
+    # Build name→index map for looking up per-parameter lb/ub/scale.
+    all_reduced_names = result.problem.get_reduced_vector(
+        result.problem.x_names
+    )
+    name_to_reduced_idx = {name: i for i, name in enumerate(all_reduced_names)}
+    x_scales_reduced = (
+        result.problem.get_reduced_vector(result.problem.x_scales)
+        if getattr(result.problem, "x_scales", None) is not None
+        else None
+    )
+
+    _show_kde = plot_type in ("kde", "both")
+    _show_rug = plot_type in ("hist", "both")
 
-    # fig, ax = plt.subplots(nr_params, figsize=size)[1]
     for idx, par_id in enumerate(param_names):
-        if plot_type == "kde":
-            # TODO: add bw_adjust as option?
-            sns.kdeplot(
-                params_fval[par_id], bw_method=bw_method, ax=par_ax[par_id]
-            )
-        elif plot_type == "hist":
-            # fixes usage of sns distplot which throws a future warning
-            sns.histplot(
-                x=params_fval[par_id], ax=par_ax[par_id], stat="density"
-            )
-            sns.rugplot(x=params_fval[par_id], ax=par_ax[par_id])
-        elif plot_type == "both":
-            sns.histplot(
-                x=params_fval[par_id],
-                kde=True,
-                ax=par_ax[par_id],
-                stat="density",
-            )
-            sns.rugplot(x=params_fval[par_id], ax=par_ax[par_id])
+        ax = par_ax[par_id]
+        vals = np.asarray(params_fval[par_id])
+        finite_vals = vals[np.isfinite(vals)]
+        par_reduced_idx = name_to_reduced_idx.get(par_id, idx)
+        lb_val = theta_lb[par_reduced_idx]
+        ub_val = theta_ub[par_reduced_idx]
+
+        bound_handle = plot_density_panel(
+            ax,
+            vals,
+            bins=bins,
+            bw_method=bw_method,
+            style=style,
+            show_hist=(plot_type in ("hist", "both")),
+            show_kde=_show_kde,
+            show_rug=_show_rug,
+            show_bounds=show_bounds,
+            lb=lb_val,
+            ub=ub_val,
+        )
 
-        par_ax[par_id].set_xlabel(param_names[idx])
-        par_ax[par_id].set_ylabel("Density")
+        legend_handles, legend_labels = [], []
+        if finite_vals.size > 0 and idx == 0:
+            if _show_kde:
+                legend_handles.append(
+                    Line2D(
+                        [0],
+                        [0],
+                        color=style["line_color"],
+                        lw=style["linewidth"],
+                    )
+                )
+                legend_labels.append("KDE")
+            if _show_rug:
+                legend_handles.append(
+                    Line2D(
+                        [0],
+                        [0],
+                        color=style["dash_color"],
+                        marker="|",
+                        lw=0,
+                        markersize=style["dash_markersize"],
+                        markeredgewidth=style["dash_linewidth"],
+                    )
+                )
+                legend_labels.append("Samples")
 
-    sns.despine()
+        if bound_handle is not None and idx == 0:
+            legend_handles.append(bound_handle)
+            legend_labels.append("Bounds")
 
-    if suptitle:
-        fig.suptitle(suptitle)
+        if legend_handles:
+            ax.legend(handles=legend_handles, labels=legend_labels)
+
+        scale = (
+            x_scales_reduced[par_reduced_idx]
+            if x_scales_reduced is not None
+            else None
+        )
+        xlabel = f"{par_id} ({scale})" if scale is not None else par_id
+        ax.set_xlabel(xlabel)
+        # y-label only on the leftmost column to avoid grid-wide repetition
+        ax.set_ylabel("Density" if idx % num_col == 0 else "")
+
+    if title is not None:
+        fig.suptitle(title)
 
     return axes
 
diff --git a/pypesto/visualize/waterfall.py b/pypesto/visualize/waterfall.py
index c25158363..13be87a62 100644
--- a/pypesto/visualize/waterfall.py
+++ b/pypesto/visualize/waterfall.py
@@ -9,6 +9,7 @@
 
 from ..C import ALL, COLOR, WATERFALL_MAX_VALUE
 from ..result import Result
+from ._style import resolve_style
 from .clust_color import assign_colors
 from .misc import (
     get_ax,
@@ -33,6 +34,7 @@ def waterfall(
     colors: COLOR | list[COLOR] | np.ndarray | None = None,
     legends: Sequence[str] | str | None = None,
     order_by_id: bool = False,
+    style_kwargs: dict | None = None,
 ) -> matplotlib.axes.Axes:
     """
     Plot waterfall plot.
@@ -71,6 +73,17 @@ def waterfall(
         the same x-axis position. Only applicable when a list of result
         objects are provided. Default behavior is to sort the function values
         of each result independently of other results.
+    style_kwargs:
+        Style overrides. Keys used by this function:
+
+        - ``cmap_discrete``, ``mle_color``, ``outlier_color`` — colours
+          of the per-start scatter dots when clustering is used (best
+          cluster, secondary clusters, isolated starts respectively).
+          Only consulted when ``colors`` is ``None``; an explicit
+          ``colors`` short-circuits clustering.
+
+        All valid keys and their defaults are listed in
+        :data:`pypesto.visualize._style._DEFAULTS`.
 
     Returns
     -------
@@ -78,6 +91,7 @@ def waterfall(
         The plot axes.
     """
     ax = get_ax(ax, size)
+    style = resolve_style(style_kwargs)
 
     if n_starts_to_zoom:
         # create zoom in
@@ -89,7 +103,9 @@ def waterfall(
         inset_axes = None
 
     # parse input
-    (results, colors, legends) = process_result_list(results, colors, legends)
+    (results, colors, legends) = process_result_list(
+        results, colors, legends, style=style
+    )
 
     # handle `order_by_id`
     if order_by_id:
@@ -153,7 +169,7 @@ def waterfall(
             fvals.sort()
 
         # assign colors
-        coloring = assign_colors(fvals, colors=colors[j])
+        coloring = assign_colors(fvals, colors=colors[j], style=style)
 
         # call lowlevel plot routine
         ax = waterfall_lowlevel(
@@ -164,6 +180,7 @@ def waterfall(
             size=size,
             colors=coloring,
             legend_text=legends[j],
+            style=style,
         )
 
         if inset_axes is not None:
@@ -172,6 +189,7 @@ def waterfall(
                 scale_y=scale_y,
                 ax=inset_axes,
                 colors=coloring[:n_starts_to_zoom],
+                style=style,
             )
             # remove the title and axes labels for the zoom in subplot
             inset_axes.set(title=None, xlabel=None, ylabel=None)
@@ -203,6 +221,7 @@ def waterfall_lowlevel(
     offset_y: float = 0.0,
     colors: COLOR | list[COLOR] | np.ndarray | None = None,
     legend_text: str | None = None,
+    style: dict | None = None,
 ) -> matplotlib.axes.Axes:
     """
     Plot waterfall plot using list of function values.
@@ -226,6 +245,10 @@ def waterfall_lowlevel(
         and colors are assigned automatically
     legend_text:
         Label for line plots
+    style:
+        Pre-resolved visualization style dict, as returned by
+        :func:`pypesto.visualize._style.resolve_style`. When ``None``, defaults
+        are used.
 
     Returns
     -------
@@ -240,7 +263,7 @@ def waterfall_lowlevel(
         colors = [colors[i] for i in start_indices]
 
     # assign colors
-    colors = assign_colors(fvals, colors=colors)
+    colors = assign_colors(fvals, colors=colors, style=style)
 
     # plot
     ax.xaxis.set_major_locator(MaxNLocator(integer=True))