diff --git a/scipy-stubs/cluster/_hierarchy.pyi b/scipy-stubs/cluster/_hierarchy.pyi
new file mode 100644
index 00000000..2dbf545f
--- /dev/null
+++ b/scipy-stubs/cluster/_hierarchy.pyi
@@ -0,0 +1,65 @@
+from typing import Any, Final, Literal, final, type_check_only
+
+import numpy as np
+import optype.numpy as onp
+
+@type_check_only
+class _CythonMixin:
+    def __reduce_cython__(self) -> tuple[Any, ...]: ...
+    def __setstate_cython__(self, /, state: tuple[Any, ...]) -> None: ...
+
+# defined in `cluster/_structures.pxi`
+@type_check_only
+@final
+class Pair:
+    key: Final[int]
+    value: Final[float]
+
+# defined in `cluster/_structures.pxi`
+class Heap(_CythonMixin):
+    index_by_key: Final[onp.Array1D[np.int32]]
+    keys_by_index: Final[onp.Array1D[np.int32]]
+    values: Final[onp.Array1D[np.float64]]
+    size: Final[int]
+
+    def __init__(self, /, values: onp.Array1D[np.float64]) -> None: ...
+    def get_min(self, /) -> Pair: ...
+    def remove_min(self, /) -> None: ...
+    def change_value(self, /, key: int, value: float) -> None: ...
+
+class LinkageUnionFind(_CythonMixin):
+    parent: Final[onp.Array1D[np.int32]]
+    size: Final[onp.Array1D[np.int32]]
+    next_label: int
+    def __init__(self, /, n: int) -> None: ...
+
+def calculate_cluster_sizes(Z: onp.Array2D[np.float64], cs: onp.Array1D[np.float64], n: int) -> None: ...
+def cluster_dist(Z: onp.Array2D[np.float64], T: onp.Array1D[np.int32], cutoff: float, n: int) -> None: ...
+def cluster_in(
+    Z: onp.Array2D[np.float64], R: onp.Array2D[np.float64], T: onp.Array1D[np.int32], cutoff: float, n: int
+) -> None: ...
+def cluster_maxclust_dist(Z: onp.Array2D[np.float64], T: onp.Array1D[np.int32], n: int, mc: int) -> None: ...
+def cluster_maxclust_monocrit(
+    Z: onp.Array2D[np.float64], MC: onp.Array1D[np.float64], T: onp.Array1D[np.int32], n: int, max_nc: int
+) -> None: ...
+def cluster_monocrit(
+    Z: onp.Array2D[np.float64], MC: onp.Array1D[np.float64], T: onp.Array1D[np.int32], cutoff: float, n: int
+) -> None: ...
+def cophenetic_distances(Z: onp.Array2D[np.float64], d: onp.Array1D[np.float64], n: int) -> None: ...
+def get_max_Rfield_for_each_cluster(
+    Z: onp.Array2D[np.float64], R: onp.Array2D[np.float64], max_rfs: onp.Array1D[np.float64], n: int, rf: int
+) -> None: ...
+def get_max_dist_for_each_cluster(Z: onp.Array2D[np.float64], MD: onp.Array1D[np.float64], n: int) -> None: ...
+def inconsistent(Z: onp.Array2D[np.float64], R: onp.Array2D[np.float64], n: int, d: float) -> None: ...
+def leaders(
+    Z: onp.Array2D[np.float64], T: onp.Array1D[np.int32], L: onp.Array1D[np.int32], M: onp.Array1D[np.int32], nc: int, n: float
+) -> None: ...
+
+#
+def linkage(
+    dists: onp.Array1D[np.float64], n: np.int64 | int, method: Literal[0, 1, 2, 3, 4, 5, 6]
+) -> onp.Array2D[np.float64]: ...
+def fast_linkage(dists: onp.Array1D[np.float64], n: int, method: Literal[0, 1, 2, 3, 4, 5, 6]) -> onp.Array2D[np.float64]: ...
+def nn_chain(dists: onp.Array1D[np.float64], n: int, method: Literal[0, 1, 2, 3, 4, 5, 6]) -> onp.Array2D[np.float64]: ...
+def mst_single_linkage(dists: onp.Array1D[np.float64], n: int) -> onp.Array2D[np.float64]: ...
+def prelist(Z: onp.Array2D[np.float64], members: onp.Array1D[np.int32], n: int) -> None: ...
diff --git a/scipy-stubs/cluster/_optimal_leaf_ordering.pyi b/scipy-stubs/cluster/_optimal_leaf_ordering.pyi
new file mode 100644
index 00000000..a7a01df2
--- /dev/null
+++ b/scipy-stubs/cluster/_optimal_leaf_ordering.pyi
@@ -0,0 +1,4 @@
+import numpy as np
+import optype.numpy as onp
+
+def optimal_leaf_ordering(Z: onp.Array2D[np.float64], D: onp.Array2D[np.float64]) -> onp.Array2D[np.float64]: ...
diff --git a/scipy-stubs/cluster/_vq.pyi b/scipy-stubs/cluster/_vq.pyi
new file mode 100644
index 00000000..37cfb105
--- /dev/null
+++ b/scipy-stubs/cluster/_vq.pyi
@@ -0,0 +1,14 @@
+from typing import TypeVar
+
+import numpy as np
+import optype.numpy as onp
+
+_FloatT = TypeVar("_FloatT", np.float32, np.float64)
+_ShapeT = TypeVar("_ShapeT", bound=tuple[int] | tuple[int, int])
+
+def update_cluster_means(
+    obs: onp.ArrayND[_FloatT, _ShapeT], labels: onp.Array1D[np.int32], nc: int
+) -> tuple[onp.ArrayND[_FloatT, _ShapeT], onp.Array1D[np.int32]]: ...
+def vq(
+    obs: onp.ArrayND[_FloatT, tuple[int] | tuple[int, int]], codes: onp.Array2D[_FloatT]
+) -> tuple[onp.Array1D[np.int32], onp.Array1D[_FloatT]]: ...