diff --git a/.github/workflows/quality-checks.yaml b/.github/workflows/quality-checks.yaml
index c1053b44..7c787550 100644
--- a/.github/workflows/quality-checks.yaml
+++ b/.github/workflows/quality-checks.yaml
@@ -16,7 +16,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest, macos-latest, windows-latest]
-        python-version: ["3.9", "3.10", "3.11", "3.12"]
+        python-version: ["3.10", "3.11", "3.12"]
       fail-fast: false
 
     runs-on: ${{ matrix.os }}
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 399859e3..cd176bbb 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,12 @@
 # CHANGELOG
 
+## v0.4.2 - 11.02.2026
+Maintance release (https://github.com/geometric-kernels/GeometricKernels/pull/173).
+* Fix compatibility with modern `plum-dispatch`/`numpy`.
+* Replace legacy `plum` typing aliases with modern Python type syntax.
+* Drop Python 3.9 support.
+* Discontinue testing for `gpflow` due to incompatibility with the latest `setuptools`.
+
 ## v0.4.1 - 15.12.2025
 * HammingGraph space by @colmont in https://github.com/geometric-kernels/GeometricKernels/pull/170
 * Added a link to a pull request example in README by @vabor112 in https://github.com/geometric-kernels/GeometricKernels/pull/171
diff --git a/Makefile b/Makefile
index d1a21268..6c3c2f5e 100644
--- a/Makefile
+++ b/Makefile
@@ -43,10 +43,6 @@ lint: sync
 
 test: sync ## Run the tests, start with the failing ones and break on first fail.
 	@$(UV_RUN) pytest -v -x --ff -rN -Wignore -s --tb=short --durations=0 --cov --cov-report=xml --cov-report=html:coverage_html tests
-	@$(UV_RUN) pytest --nbmake --nbmake-kernel=python3 --durations=0 --nbmake-timeout=1000 --ignore=notebooks/frontends/GPJax.ipynb notebooks/
-	@if [ "$(UV_PYTHON)" = "python3.9" ]; then \
-		echo "Skipping GPJax notebook on python3.9"; \
-	else \
-		$(UV_RUN) pytest --nbmake --nbmake-kernel=python3 --durations=0 --nbmake-timeout=1000 notebooks/frontends/GPJax.ipynb; \
-	fi;
+	# gpflow is ignored due to incompatibility with the recent setuptools
+	@$(UV_RUN) pytest --nbmake --nbmake-kernel=python3 --durations=0 --nbmake-timeout=1000 --ignore=notebooks/frontends/GPflow.ipynb notebooks/
 	@echo -e "$(SUCCESS)Tests done$(RESET)"
diff --git a/docs/index.rst b/docs/index.rst
index d5698e32..3cb2824c 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -33,7 +33,7 @@ Before doing anything, you might want to create and activate a new virtual envir
 
    uv venv --python python[version] [venv_dir]
 
-where [env_dir] is the directory (default is `.venv`) of the environment and [version] is the version of Python you want to use, we currently support 3.9, 3.10, 3.11, 3.12.
+where [env_dir] is the directory (default is `.venv`) of the environment and [version] is the version of Python you want to use, we currently support 3.10, 3.11, 3.12.
 
 [Optional] activate the environment. However, this is not strictly necessary for `uv`. Instead, use tools like `uv run python` to run Python inside the environment. See `uv documentation <https://docs.astral.sh/uv/>` for more details.
 
@@ -65,7 +65,7 @@ where [env_dir] is the directory (default is `.venv`) of the environment and [ve
    conda create -n [env_name] python=[version]
    conda activate [env_name]
 
-where [env_name] is the name of the environment and [version] is the version of Python you want to use, we currently support 3.9, 3.10, 3.11, 3.12.
+where [env_name] is the name of the environment and [version] is the version of Python you want to use, we currently support 3.10, 3.11, 3.12.
 
 .. raw:: html
 
@@ -449,4 +449,3 @@ Please also consider citing the theoretical papers the library is based on. You
    API reference <autoapi/geometric_kernels/index>
    Bibliography <bibliography>
    GitHub <https://github.com/geometric-kernels/GeometricKernels>
-
diff --git a/geometric_kernels/frontends/gpflow.py b/geometric_kernels/frontends/gpflow.py
index 172f483d..87fadea4 100644
--- a/geometric_kernels/frontends/gpflow.py
+++ b/geometric_kernels/frontends/gpflow.py
@@ -6,13 +6,40 @@
 :doc:`frontends/GPflow.ipynb </examples/frontends/GPflow>` notebook.
 """
 
-import gpflow
 import numpy as np
 import tensorflow as tf
 from beartype.typing import List, Optional, Union
-from gpflow.base import TensorType
-from gpflow.kernels.base import ActiveDims
-from gpflow.utilities import positive
+
+try:
+    import gpflow
+    from gpflow.base import TensorType
+    from gpflow.kernels.base import ActiveDims
+    from gpflow.utilities import positive
+except ImportError as error:
+    seen = set()
+    current: BaseException | None = error
+    missing_pkg_resources = False
+    while current is not None and id(current) not in seen:
+        seen.add(id(current))
+        if getattr(current, "name", None) == "pkg_resources":
+            missing_pkg_resources = True
+            break
+        text = str(current)
+        if "pkg_resources" in text and (
+            "No module named" in text or "cannot import name" in text
+        ):
+            missing_pkg_resources = True
+            break
+        current = current.__cause__ or current.__context__
+
+    if missing_pkg_resources:
+        raise ImportError(
+            "Importing `gpflow` failed because it depends on `pkg_resources`. "
+            "`pkg_resources` was removed from `setuptools==82.0.0`. "
+            "You can try pinning an older `setuptools` version, but this setup "
+            "is no longer tested by GeometricKernels."
+        ) from error
+    raise
 
 from geometric_kernels.kernels import BaseGeometricKernel
 from geometric_kernels.spaces import Space
diff --git a/geometric_kernels/lab_extras/extras.py b/geometric_kernels/lab_extras/extras.py
index 8d75a4fd..15f673f7 100644
--- a/geometric_kernels/lab_extras/extras.py
+++ b/geometric_kernels/lab_extras/extras.py
@@ -1,8 +1,6 @@
 import lab as B
-from beartype.typing import List
 from lab import dispatch
 from lab.util import abstract
-from plum import Union
 from scipy.sparse import spmatrix
 
 
@@ -23,7 +21,7 @@ def take_along_axis(a: B.Numeric, index: B.Numeric, axis: int = 0):
 
 @dispatch
 @abstract()
-def from_numpy(_: B.Numeric, b: Union[List, B.Numeric]):
+def from_numpy(_: B.Numeric, b: list | B.Numeric):
     """
     Converts the array `b` to a tensor of the same backend as `_`.
 
@@ -290,7 +288,7 @@ def eigvalsh(x: B.Numeric):
 
 @dispatch
 @abstract()
-def reciprocal_no_nan(x: Union[B.Numeric, spmatrix]):
+def reciprocal_no_nan(x: B.Numeric | spmatrix):
     """
     Return element-wise reciprocal (1/x). Whenever x = 0 puts 1/x = 0.
 
@@ -357,9 +355,7 @@ def bool_like(reference: B.Numeric):
     """
 
 
-def smart_cast(
-    dtype: Union[B.Bool, B.Int, B.Float, B.Complex, B.Numeric], x: B.Numeric
-):
+def smart_cast(dtype: B.Bool | B.Int | B.Float | B.Complex | B.Numeric, x: B.Numeric):
     """
     Return `x` cast to the `dtype` abstract data type.
 
diff --git a/geometric_kernels/lab_extras/jax/extras.py b/geometric_kernels/lab_extras/jax/extras.py
index 974fae99..6455325d 100644
--- a/geometric_kernels/lab_extras/jax/extras.py
+++ b/geometric_kernels/lab_extras/jax/extras.py
@@ -1,14 +1,15 @@
+from typing import TypeAlias
+
 import jax.numpy as jnp
 import lab as B
-from beartype.typing import List
 from lab import dispatch
-from plum import Union, convert
+from plum import convert
 
-_Numeric = Union[B.Number, B.JAXNumeric]
+_Numeric: TypeAlias = B.Number | B.JAXNumeric
 
 
 @dispatch
-def take_along_axis(a: Union[_Numeric, B.Numeric], index: _Numeric, axis: int = 0) -> _Numeric:  # type: ignore
+def take_along_axis(a: _Numeric | B.Numeric, index: _Numeric, axis: int = 0) -> _Numeric:  # type: ignore
     """
     Gathers elements of `a` along `axis` at `index` locations.
     """
@@ -18,7 +19,7 @@ def take_along_axis(a: Union[_Numeric, B.Numeric], index: _Numeric, axis: int =
 
 
 @dispatch
-def from_numpy(_: B.JAXNumeric, b: Union[List, B.NPNumeric, B.Number, B.JAXNumeric]):  # type: ignore
+def from_numpy(_: B.JAXNumeric, b: list | B.NPNumeric | B.Number | B.JAXNumeric):  # type: ignore
     """
     Converts the array `b` to a tensor of the same backend as `a`
     """
diff --git a/geometric_kernels/lab_extras/numpy/extras.py b/geometric_kernels/lab_extras/numpy/extras.py
index 01c69168..5167823c 100644
--- a/geometric_kernels/lab_extras/numpy/extras.py
+++ b/geometric_kernels/lab_extras/numpy/extras.py
@@ -1,11 +1,12 @@
+from typing import TypeAlias
+
 import lab as B
 import numpy as np
-from beartype.typing import Any, List, Optional
+from beartype.typing import Any
 from lab import dispatch
-from plum import Union
 from scipy.sparse import spmatrix
 
-_Numeric = Union[B.Number, B.NPNumeric]
+_Numeric: TypeAlias = B.Number | B.NPNumeric
 
 
 @dispatch
@@ -17,7 +18,7 @@ def take_along_axis(a: _Numeric, index: _Numeric, axis: int = 0) -> _Numeric:  #
 
 
 @dispatch
-def from_numpy(_: B.NPNumeric, b: Union[List, B.NPNumeric, B.Number]):  # type: ignore
+def from_numpy(_: B.NPNumeric, b: list | B.NPNumeric | B.Number):  # type: ignore
     """
     Converts the array `b` to a tensor of the same backend as `a`
     """
@@ -33,7 +34,7 @@ def trapz(y: _Numeric, x: _Numeric, dx: _Numeric = 1.0, axis: int = -1):  # type
 
 
 @dispatch
-def norm(x: _Numeric, ord: Optional[Any] = None, axis: Optional[int] = None):  # type: ignore
+def norm(x: _Numeric, ord: Any | None = None, axis: int | None = None):  # type: ignore
     """
     Matrix or vector norm.
     """
diff --git a/geometric_kernels/lab_extras/tensorflow/extras.py b/geometric_kernels/lab_extras/tensorflow/extras.py
index 594be7f9..90606723 100644
--- a/geometric_kernels/lab_extras/tensorflow/extras.py
+++ b/geometric_kernels/lab_extras/tensorflow/extras.py
@@ -1,13 +1,12 @@
-import sys
+from typing import TypeAlias
 
 import lab as B
 import tensorflow as tf
 import tensorflow_probability as tfp
-from beartype.typing import Any, List, Optional
+from beartype.typing import Any
 from lab import dispatch
-from plum import Union
 
-_Numeric = Union[B.Number, B.TFNumeric, B.NPNumeric]
+_Numeric: TypeAlias = B.Number | B.TFNumeric | B.NPNumeric
 
 
 @dispatch
@@ -15,15 +14,11 @@ def take_along_axis(a: _Numeric, index: _Numeric, axis: int = 0) -> _Numeric:  #
     """
     Gathers elements of `a` along `axis` at `index` locations.
     """
-    if sys.version_info[:2] <= (3, 9):
-        index = tf.cast(index, tf.int32)
-    return tf.experimental.numpy.take_along_axis(
-        a, index, axis=axis
-    )  # the absence of explicit cast to int64 causes an error for Python 3.9 and below
+    return tf.experimental.numpy.take_along_axis(a, index, axis=axis)
 
 
 @dispatch
-def from_numpy(_: B.TFNumeric, b: Union[List, B.Numeric, B.NPNumeric, B.TFNumeric]):  # type: ignore
+def from_numpy(_: B.TFNumeric, b: list | B.Numeric | B.NPNumeric | B.TFNumeric):  # type: ignore
     """
     Converts the array `b` to a tensor of the same backend as `a`
     """
@@ -39,7 +34,7 @@ def trapz(y: _Numeric, x: _Numeric, dx=None, axis=-1):  # type: ignore
 
 
 @dispatch
-def norm(x: _Numeric, ord: Optional[Any] = None, axis: Optional[int] = None):  # type: ignore
+def norm(x: _Numeric, ord: Any | None = None, axis: int | None = None):  # type: ignore
     """
     Matrix or vector norm.
     """
diff --git a/geometric_kernels/lab_extras/torch/extras.py b/geometric_kernels/lab_extras/torch/extras.py
index 2f80d42c..771aa205 100644
--- a/geometric_kernels/lab_extras/torch/extras.py
+++ b/geometric_kernels/lab_extras/torch/extras.py
@@ -1,14 +1,15 @@
+from typing import TypeAlias
+
 import lab as B
 import torch
-from beartype.typing import Any, List, Optional
+from beartype.typing import Any
 from lab import dispatch
-from plum import Union
 
-_Numeric = Union[B.Number, B.TorchNumeric]
+_Numeric: TypeAlias = B.Number | B.TorchNumeric
 
 
 @dispatch
-def take_along_axis(a: Union[_Numeric, B.Numeric], index: _Numeric, axis: int = 0) -> _Numeric:  # type: ignore
+def take_along_axis(a: _Numeric | B.Numeric, index: B.TorchNumeric, axis: int = 0) -> _Numeric:  # type: ignore
     """
     Gathers elements of `a` along `axis` at `index` locations.
     """
@@ -19,7 +20,7 @@ def take_along_axis(a: Union[_Numeric, B.Numeric], index: _Numeric, axis: int =
 
 @dispatch
 def from_numpy(
-    a: B.TorchNumeric, b: Union[List, B.Number, B.NPNumeric, B.TorchNumeric]
+    a: B.TorchNumeric, b: list | B.Number | B.NPNumeric | B.TorchNumeric
 ):  # type: ignore
     """
     Converts the array `b` to a tensor of the same backend as `a`
@@ -38,7 +39,7 @@ def trapz(y: B.TorchNumeric, x: _Numeric, dx: _Numeric = 1.0, axis: int = -1):
 
 
 @dispatch
-def norm(x: _Numeric, ord: Optional[Any] = None, axis: Optional[int] = None):  # type: ignore
+def norm(x: _Numeric, ord: Any | None = None, axis: int | None = None):  # type: ignore
     """
     Matrix or vector norm.
     """
diff --git a/pyproject.toml b/pyproject.toml
index 9bfe07ed..a8432136 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -11,7 +11,6 @@ description="A Python Package offering geometric kernels in NumPy, TensorFlow, P
 readme = "README.md"
 classifiers = [
     "License :: OSI Approved :: Apache Software License",
-    "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
@@ -21,21 +20,21 @@ classifiers = [
 keywords=[
     "geometric-kernels",
 ]
-requires-python = ">=3.9"
+requires-python = ">=3.10"
 dependencies = [
-    "backends", # >=1.7",
+    "backends>=1.8.0",
     "einops",
     "geomstats",
-    "numpy>=2.0",
+    "numpy>=2.0,<2.4",   # constraining from above due to geomstats==2.8.0 incompatibility
     "opt-einsum",
-    "plum-dispatch",
+    "plum-dispatch>=2.6.0",
     "potpourri3d",
     "robust_laplacian",
     "scipy>=1.3",
     "spherical-harmonics-basis",
     "sympy~=1.13",
 ]
-version="0.4.1"
+version="0.4.2"
 
 [project.urls]
 Documentation = "https://geometric-kernels.github.io/"
@@ -58,7 +57,7 @@ allow_redefinition = true
 
 [tool.black]
 line-length = 88
-target-version = ['py39', 'py310', 'py311', 'py312']
+target-version = ['py310', 'py311', 'py312']
 
 
 [tool.uv]
@@ -106,5 +105,5 @@ dev = [
   'jaxlib',
   'jaxtyping',
   'optax',
-  'gpjax>=0.12.2; python_version >= "3.10"',  # gpjax is not supported on python-3.9 or older.
+  'gpjax>=0.12.2; python_version >= "3.10"',  # gpjax requires python >= 3.10.
 ]
diff --git a/tests/spaces/test_hamming_graph.py b/tests/spaces/test_hamming_graph.py
index 55cefa34..b78bbfd9 100644
--- a/tests/spaces/test_hamming_graph.py
+++ b/tests/spaces/test_hamming_graph.py
@@ -1,7 +1,6 @@
 import lab as B
 import numpy as np
 import pytest
-from plum import Tuple
 
 from geometric_kernels.kernels import MaternGeometricKernel
 from geometric_kernels.spaces import HammingGraph, HypercubeGraph
@@ -11,7 +10,7 @@
 
 
 @pytest.fixture(params=[(1, 2), (2, 2), (5, 2), (10, 2), (10, 4)])
-def inputs(request) -> Tuple[B.Numeric]:
+def inputs(request) -> tuple[B.Numeric]:
     """
     Returns a tuple (space, eigenfunctions, X, X2, weights) where:
     - space is a HammingGraph object with (dim, n_cat) equal to request.param,
diff --git a/tests/spaces/test_hypercube_graph.py b/tests/spaces/test_hypercube_graph.py
index cb5ac9c9..fa820137 100644
--- a/tests/spaces/test_hypercube_graph.py
+++ b/tests/spaces/test_hypercube_graph.py
@@ -1,7 +1,6 @@
 import lab as B
 import numpy as np
 import pytest
-from plum import Tuple
 
 from geometric_kernels.kernels import MaternGeometricKernel
 from geometric_kernels.spaces import HypercubeGraph
@@ -11,7 +10,7 @@
 
 
 @pytest.fixture(params=[1, 2, 3, 5, 10])
-def inputs(request) -> Tuple[B.Numeric]:
+def inputs(request) -> tuple[B.Numeric]:
     """
     Returns a tuple (space, eigenfunctions, X, X2) where:
     - space is a HypercubeGraph object with dimension equal to request.param,