diff --git a/src/scanpy/plotting/_tools/scatterplots.py b/src/scanpy/plotting/_tools/scatterplots.py
index 30356d8b04..a764dba57f 100644
--- a/src/scanpy/plotting/_tools/scatterplots.py
+++ b/src/scanpy/plotting/_tools/scatterplots.py
@@ -294,9 +294,8 @@ def embedding(  # noqa: PLR0912, PLR0913, PLR0915
         elif sort_order and color_type == "cat":
             # Null points go on bottom
             order = np.argsort(~pd.isnull(color_source_vector), kind="stable")
-        # Set orders
-        if isinstance(size, np.ndarray):
-            size = np.array(size)[order]
+        # `size` is not a loop variable, so don’t overwrite it here
+        size_plot = size[order] if isinstance(size, np.ndarray) else size
         color_source_vector = color_source_vector[order]
         color_vector = color_vector[order]
         coords = basis_values[:, dims][order, :]
@@ -348,11 +347,10 @@ def embedding(  # noqa: PLR0912, PLR0913, PLR0915
             )
         else:
             scatter = (
-                partial(ax.scatter, s=size, plotnonfinite=True)
+                partial(ax.scatter, s=size_plot, plotnonfinite=True)
                 if scale_factor is None
-                else partial(
-                    circles, s=size, ax=ax, scale_factor=scale_factor
-                )  # size in circles is radius
+                # size in circles is radius
+                else partial(circles, s=size_plot, ax=ax, scale_factor=scale_factor)
             )
 
             if add_outline:
@@ -366,7 +364,7 @@ def embedding(  # noqa: PLR0912, PLR0913, PLR0915
                 # with some transparency.
 
                 bg_width, gap_width = outline_width
-                point = np.sqrt(size)
+                point = np.sqrt(size_plot)
                 gap_size = (point + (point * gap_width) * 2) ** 2
                 bg_size = (np.sqrt(gap_size) + (point * bg_width) * 2) ** 2
                 # the default black and white colors can be changes using
diff --git a/tests/test_plotting.py b/tests/test_plotting.py
index af35236c6a..7f327e0b82 100644
--- a/tests/test_plotting.py
+++ b/tests/test_plotting.py
@@ -1875,6 +1875,38 @@ def test_umap_mask_no_modification():
     pd.testing.assert_series_equal(pbmc.obs["louvain"], data_copy)
 
 
+def test_scatter_size_not_mutated_across_panels():
+    """Per-point size array must not be cumulatively reordered across subplots.
+
+    Regression test for https://github.com/scverse/scanpy/issues/4024
+    Uses 3 panels (categorical + two continuous) to exercise cumulative reorder.
+    """
+    pbmc = pbmc3k_processed()
+    rng = np.random.default_rng(0)
+    sizes = rng.uniform(10, 200, size=pbmc.n_obs)
+    sizes_original = sizes.copy()
+
+    axes = sc.pl.umap(
+        pbmc,
+        color=["louvain", "n_genes", "n_counts"],
+        size=sizes,
+        show=False,
+    )
+
+    # The input array must not be modified
+    np.testing.assert_array_equal(sizes, sizes_original)
+
+    # Each panel must plot the correct per-point sizes (just reordered by
+    # z-order, not cumulatively scrambled).  Sorting makes the comparison
+    # independent of z-ordering.
+    expected_sorted = np.sort(sizes)
+    for ax in axes:
+        plotted = ax.collections[0].get_sizes()
+        np.testing.assert_allclose(np.sort(plotted), expected_sorted)
+
+    plt.close()
+
+
 def test_string_mask(tmp_path, check_same_image):
     """Check that the same mask given as string or bool array provides the same result."""
     pbmc = pbmc3k_processed()