diff --git a/contrib/pov-nodes/nodal_points.py b/contrib/pov-nodes/nodal_points.py
index 158e92f6..fa31b74f 100644
--- a/contrib/pov-nodes/nodal_points.py
+++ b/contrib/pov-nodes/nodal_points.py
@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import numpy as np
 from pov import (
     Background,
diff --git a/contrib/pov-nodes/pov.py b/contrib/pov-nodes/pov.py
index 93fd3a06..8f3de37c 100644
--- a/contrib/pov-nodes/pov.py
+++ b/contrib/pov-nodes/pov.py
@@ -1,4 +1,6 @@
 # http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/205451
+from __future__ import annotations
+
 import os
 
 import numpy as np
diff --git a/contrib/weights-from-festa-sommariva.py b/contrib/weights-from-festa-sommariva.py
index 9b197911..baab0f9e 100644
--- a/contrib/weights-from-festa-sommariva.py
+++ b/contrib/weights-from-festa-sommariva.py
@@ -11,6 +11,7 @@
 
 These rules are licensed under the GPL and are not included by default.
 """
+from __future__ import annotations
 
 import re
 
diff --git a/contrib/weights-from-fortran.py b/contrib/weights-from-fortran.py
index a26308b9..54105806 100644
--- a/contrib/weights-from-fortran.py
+++ b/contrib/weights-from-fortran.py
@@ -1,4 +1,5 @@
 #! /usr/bin/env python
+from __future__ import annotations
 
 import re
 
diff --git a/contrib/weights-from-jaskowiec-sukumar.py b/contrib/weights-from-jaskowiec-sukumar.py
index 028c5486..43253bc3 100644
--- a/contrib/weights-from-jaskowiec-sukumar.py
+++ b/contrib/weights-from-jaskowiec-sukumar.py
@@ -10,12 +10,12 @@
 The quadrature rules can be found in the supplemental material at the same
 URL and should be a file named ``nme6528-sup-0002-supinfo.tar``.
 """
+from __future__ import annotations
 
 import pathlib
 import re
 import sys
 import tarfile
-from typing import Optional
 
 
 FILE_ORDER_REGEX = re.compile(r"cubature_tet_sym_p(\d{1,2})_n\d{1,3}_expand.txt")
@@ -64,7 +64,7 @@ def extract_rules_if_not(path: pathlib.Path) -> pathlib.Path:
 
 
 def generate_jaskowiec_sukumar_quadrature_rules(
-    infolder: pathlib.Path, outfile: Optional[pathlib.Path] = None
+    infolder: pathlib.Path, outfile: pathlib.Path | None = None
 ) -> None:
     import numpy as np
 
diff --git a/contrib/weights-from-txt.py b/contrib/weights-from-txt.py
index 2399e0c6..daf1297e 100755
--- a/contrib/weights-from-txt.py
+++ b/contrib/weights-from-txt.py
@@ -1,4 +1,5 @@
 #! /usr/bin/env python
+from __future__ import annotations
 
 import sys
 
diff --git a/contrib/weights-from-witherden-vincent.py b/contrib/weights-from-witherden-vincent.py
index 1bae9d2a..f4ae78c0 100644
--- a/contrib/weights-from-witherden-vincent.py
+++ b/contrib/weights-from-witherden-vincent.py
@@ -9,6 +9,7 @@
 
 They are available under a Creative Commons license.
 """
+from __future__ import annotations
 
 import os
 
diff --git a/doc/conf.py b/doc/conf.py
index baea0041..f7043d0b 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 from importlib import metadata
 from urllib.request import urlopen
 
diff --git a/examples/derivative.py b/examples/derivative.py
index c64ca774..93eaf531 100644
--- a/examples/derivative.py
+++ b/examples/derivative.py
@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import numpy as np
 
 import modepy as mp
diff --git a/examples/plot-basis.py b/examples/plot-basis.py
index 25dc1410..0ce58065 100644
--- a/examples/plot-basis.py
+++ b/examples/plot-basis.py
@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import matplotlib.pyplot as pt
 import numpy as np
 
diff --git a/examples/plot-nodes.py b/examples/plot-nodes.py
index c6271f01..480f5a5f 100644
--- a/examples/plot-nodes.py
+++ b/examples/plot-nodes.py
@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import matplotlib.pyplot as pt
 import numpy as np
 
diff --git a/examples/plot-quadrature-nodes.py b/examples/plot-quadrature-nodes.py
index a6c801b5..4507c236 100644
--- a/examples/plot-quadrature-nodes.py
+++ b/examples/plot-quadrature-nodes.py
@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import matplotlib.pyplot as plt
 import numpy as np
 
diff --git a/examples/plot-warp.py b/examples/plot-warp.py
index e0634d49..3c4653c3 100644
--- a/examples/plot-warp.py
+++ b/examples/plot-warp.py
@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import matplotlib.pyplot as pt
 import numpy as np
 
diff --git a/modepy/__init__.py b/modepy/__init__.py
index fe800f8c..c591842a 100644
--- a/modepy/__init__.py
+++ b/modepy/__init__.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = """Copyright (C) 2009, 2010, 2013 Andreas Kloeckner, Tim Warburton,
     Jan Hesthaven, Xueyu Zhu"""
 
diff --git a/modepy/matrices.py b/modepy/matrices.py
index d6534d24..659f9ad1 100644
--- a/modepy/matrices.py
+++ b/modepy/matrices.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2013 Andreas Kloeckner"
 
 __license__ = """
@@ -21,7 +24,7 @@
 """
 
 
-from typing import Callable, Optional, Sequence, Tuple, Union
+from typing import Callable, Sequence
 from warnings import warn
 
 import numpy as np
@@ -117,7 +120,7 @@ def vandermonde(
 def multi_vandermonde(
             functions: Sequence[Callable[[np.ndarray], Sequence[np.ndarray]]],
             nodes: np.ndarray
-        ) -> Tuple[np.ndarray, ...]:
+        ) -> tuple[np.ndarray, ...]:
     """Evaluate multiple (generalized) Vandermonde matrices.
 
     The Vandermonde Matrix is given by :math:`V_{i,j} := f_j(x_i)`
@@ -213,8 +216,8 @@ def differentiation_matrices(
             basis: Sequence[Callable[[np.ndarray], np.ndarray]],
             grad_basis: Sequence[Callable[[np.ndarray], Sequence[np.ndarray]]],
             nodes: np.ndarray,
-            from_nodes: Optional[np.ndarray] = None
-        ) -> Tuple[np.ndarray, ...]:
+            from_nodes: np.ndarray | None = None
+        ) -> tuple[np.ndarray, ...]:
     """Return matrices carrying out differentiation on nodal values in the
     :math:`(r,s,t)` unit directions. (See :ref:`tri-coords` and
     :ref:`tet-coords`.)
@@ -253,7 +256,7 @@ def differentiation_matrices(
 def diff_matrices(
             basis: Basis,
             nodes: np.ndarray,
-            from_nodes: Optional[np.ndarray] = None
+            from_nodes: np.ndarray | None = None
         ):
     """Like :func:`differentiation_matrices`, but for a given :class:`~modepy.Basis`.
 
@@ -264,7 +267,7 @@ def diff_matrices(
 
 
 def diff_matrix_permutation(
-            node_tuples: Sequence[Tuple[int, ...]],
+            node_tuples: Sequence[tuple[int, ...]],
             ref_axis: int
         ) -> np.ndarray:
     """Return a :mod:`numpy` array *permutation* of integers so that::
@@ -290,7 +293,7 @@ def diff_matrix_permutation(
 
 
 def inverse_mass_matrix(
-            basis: Union[Basis, Sequence[Callable[[np.ndarray], np.ndarray]]],
+            basis: Basis | Sequence[Callable[[np.ndarray], np.ndarray]],
             nodes: np.ndarray
         ) -> np.ndarray:
     """Return a matrix :math:`A=M^{-1}`, which is the inverse of the one returned
@@ -325,7 +328,7 @@ def inverse_mass_matrix(
 
 
 def mass_matrix(
-            basis: Union[Basis, Sequence[Callable[[np.ndarray], np.ndarray]]],
+            basis: Basis | Sequence[Callable[[np.ndarray], np.ndarray]],
             nodes: np.ndarray
         ) -> np.ndarray:
     r"""Return a mass matrix :math:`M`, which obeys
@@ -370,7 +373,7 @@ def modal_quad_mass_matrix(
 def nodal_quad_mass_matrix(
             quadrature: Quadrature,
             test_functions: Sequence[Callable[[np.ndarray], np.ndarray]],
-            nodes: Optional[np.ndarray] = None,
+            nodes: np.ndarray | None = None,
         ) -> np.ndarray:
     r"""Using the *quadrature*, provide a matrix :math:`M` that
     satisfies:
diff --git a/modepy/modal_decay.py b/modepy/modal_decay.py
index 2adfb5ef..c8a44be0 100644
--- a/modepy/modal_decay.py
+++ b/modepy/modal_decay.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2010-2012 Andreas Kloeckner"
 
 __license__ = """
@@ -20,7 +23,6 @@
 THE SOFTWARE.
 """
 
-from typing import Tuple
 
 import numpy as np
 import numpy.linalg as la
@@ -88,7 +90,7 @@ def simplex_interp_error_coefficient_estimator_matrix(
 
 
 def make_mode_number_vector(
-        mode_order_tuples: Tuple[Tuple[int, ...], ...],
+        mode_order_tuples: tuple[tuple[int, ...], ...],
         ignored_modes: int) -> np.ndarray:
     node_cnt = len(mode_order_tuples)
 
@@ -152,7 +154,7 @@ def skyline_pessimize(modal_values: np.ndarray) -> np.ndarray:
 
 def fit_modal_decay(
         coeffs: np.ndarray, dims: int, n: int,
-        ignored_modes: int = 1) -> Tuple[np.ndarray, np.ndarray]:
+        ignored_modes: int = 1) -> tuple[np.ndarray, np.ndarray]:
     """Fit a curve to the modal decay on each element.
 
     :arg coeffs: an array of shape ``(num_elements, num_modes)`` containing modal
diff --git a/modepy/modes.py b/modepy/modes.py
index 385fea68..34535946 100644
--- a/modepy/modes.py
+++ b/modepy/modes.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = ("Copyright (C) 2009, 2010, 2013 "
 "Andreas Kloeckner, Tim Warburton, Jan Hesthaven, Xueyu Zhu")
 
@@ -29,12 +32,9 @@
     Callable,
     Hashable,
     Iterable,
-    List,
-    Optional,
     Sequence,
     Tuple,
     TypeVar,
-    Union,
 )
 
 import numpy as np
@@ -276,7 +276,7 @@ def scaled_jacobi(alpha: float, beta: float, n: int, x: RealValueT) -> RealValue
 
 def _rstoab(
         r: RealValueT, s: RealValueT,
-        tol: float = 1.0e-12) -> Tuple[RealValueT, RealValueT]:
+        tol: float = 1.0e-12) -> tuple[RealValueT, RealValueT]:
     """Transfer from (r, s) -> (a, b) coordinates in triangle."""
 
     # We may divide by zero below (or close to it), but we won't use the
@@ -287,7 +287,7 @@ def _rstoab(
     return a, b
 
 
-def pkdo_2d(order: Tuple[int, int], rs: np.ndarray) -> np.ndarray:
+def pkdo_2d(order: tuple[int, int], rs: np.ndarray) -> np.ndarray:
     """Evaluate a 2D orthonormal (with weight 1) polynomial on the unit simplex.
 
     :arg order: A tuple *(i, j)* representing the order of the polynomial.
@@ -311,8 +311,8 @@ def pkdo_2d(order: Tuple[int, int], rs: np.ndarray) -> np.ndarray:
 
 
 def grad_pkdo_2d(
-        order: Tuple[int, int],
-        rs: np.ndarray) -> Tuple[RealValueT, RealValueT]:
+        order: tuple[int, int],
+        rs: np.ndarray) -> tuple[RealValueT, RealValueT]:
     """Evaluate the derivatives of :func:`pkdo_2d`.
 
     :arg order: A tuple *(i, j)* representing the order of the polynomial.
@@ -368,7 +368,7 @@ def grad_pkdo_2d(
 
 def _rsttoabc(
         r: RealValueT, s: RealValueT, t: RealValueT,
-        tol: float = 1.0e-10) -> Tuple[RealValueT, RealValueT, RealValueT]:
+        tol: float = 1.0e-10) -> tuple[RealValueT, RealValueT, RealValueT]:
     # We may divide by zero below (or close to it), but we won't use the
     # results because of the conditional. Silence the resulting numpy warnings.
     with np.errstate(all="ignore"):
@@ -379,7 +379,7 @@ def _rsttoabc(
     return a, b, c
 
 
-def pkdo_3d(order: Tuple[int, int, int], rst: np.ndarray) -> np.ndarray:
+def pkdo_3d(order: tuple[int, int, int], rst: np.ndarray) -> np.ndarray:
     """Evaluate a 2D orthonormal (with weight 1) polynomial on the unit simplex.
 
     :arg order: A tuple *(i, j, k)* representing the order of the polynomial.
@@ -405,8 +405,8 @@ def pkdo_3d(order: Tuple[int, int, int], rst: np.ndarray) -> np.ndarray:
 
 
 def grad_pkdo_3d(
-        order: Tuple[int, int, int],
-        rst: np.ndarray) -> Tuple[RealValueT, RealValueT, RealValueT]:
+        order: tuple[int, int, int],
+        rst: np.ndarray) -> tuple[RealValueT, RealValueT, RealValueT]:
     """Evaluate the derivatives of :func:`pkdo_3d`.
 
     :arg order: A tuple *(i, j, k)* representing the order of the polynomial.
@@ -474,7 +474,7 @@ def grad_pkdo_3d(
 
 # {{{ monomials
 
-def monomial(order: Tuple[int, ...], rst: np.ndarray) -> np.ndarray:
+def monomial(order: tuple[int, ...], rst: np.ndarray) -> np.ndarray:
     """Evaluate the monomial of order *order* at the points *rst*.
 
     :arg order: A tuple *(i, j,...)* representing the order of the polynomial.
@@ -488,7 +488,7 @@ def monomial(order: Tuple[int, ...], rst: np.ndarray) -> np.ndarray:
     return product(rst[i] ** order[i] for i in range(dim))
 
 
-def grad_monomial(order: Tuple[int, ...], rst: np.ndarray) -> Tuple[RealValueT, ...]:
+def grad_monomial(order: tuple[int, ...], rst: np.ndarray) -> tuple[RealValueT, ...]:
     """Evaluate the derivative of the monomial of order *order* at the points *rst*.
 
     :arg order: A tuple *(i, j,...)* representing the order of the polynomial.
@@ -551,9 +551,9 @@ class _TensorProductBasisFunction:
     """
 
     def __init__(self,
-            multi_index: Tuple[Hashable, ...],
-            functions: Tuple[Callable[[np.ndarray], np.ndarray], ...], *,
-            dims_per_function: Tuple[int, ...]) -> None:
+            multi_index: tuple[Hashable, ...],
+            functions: tuple[Callable[[np.ndarray], np.ndarray], ...], *,
+            dims_per_function: tuple[int, ...]) -> None:
         assert len(dims_per_function) == len(functions)
 
         self.multi_index = multi_index
@@ -634,11 +634,11 @@ class _TensorProductGradientBasisFunction:
     """
 
     def __init__(self,
-            multi_index: Tuple[int, ...],
-            derivatives: Tuple[Tuple[
-                Callable[[np.ndarray], Union[np.ndarray, Tuple[np.ndarray, ...]]],
+            multi_index: tuple[int, ...],
+            derivatives: tuple[tuple[
+                Callable[[np.ndarray], np.ndarray | tuple[np.ndarray, ...]],
                 ...], ...], *,
-            dims_per_function: Tuple[int, ...]) -> None:
+            dims_per_function: tuple[int, ...]) -> None:
         assert all(len(dims_per_function) == len(df) for df in derivatives)
 
         self.multi_index = multi_index
@@ -684,10 +684,10 @@ def __repr__(self):
 # {{{ conversion to symbolic
 
 def symbolicize_function(
-        f: Callable[[RealValueT], Union[RealValueT, Tuple[RealValueT, ...]]],
+        f: Callable[[RealValueT], RealValueT | tuple[RealValueT, ...]],
         dim: int,
         ref_coord_var_name: str = "r",
-        ) -> Union[RealValueT, Tuple[RealValueT, ...]]:
+        ) -> RealValueT | tuple[RealValueT, ...]:
     """For a function *f* (basis or gradient) returned by one of the functions in
     this module, return a :mod:`pymbolic` expression representing the
     same function.
@@ -743,7 +743,7 @@ def orthonormality_weight(self) -> float:
 
     @property
     @abstractmethod
-    def mode_ids(self) -> Tuple[Hashable, ...]:
+    def mode_ids(self) -> tuple[Hashable, ...]:
         """A tuple of of mode (basis function) identifiers, one for
         each basis function. Mode identifiers should generally be viewed
         as opaque. They are hashable. For some bases, they are tuples of
@@ -753,7 +753,7 @@ def mode_ids(self) -> Tuple[Hashable, ...]:
 
     @property
     @abstractmethod
-    def functions(self) -> Tuple[BasisFunctionType, ...]:
+    def functions(self) -> tuple[BasisFunctionType, ...]:
         """A tuple of (callable) basis functions of length matching
         :attr:`mode_ids`.  Each function takes a vector of :math:`(r, s, t)`
         reference coordinates (depending on dimension) as input.
@@ -761,7 +761,7 @@ def functions(self) -> Tuple[BasisFunctionType, ...]:
 
     @property
     @abstractmethod
-    def gradients(self) -> Tuple[BasisGradientType, ...]:
+    def gradients(self) -> tuple[BasisGradientType, ...]:
         """A tuple of (callable) basis functions of length matching
         :attr:`mode_ids`.  Each function takes a vector of :math:`(r, s, t)`
         reference coordinates (depending on dimension) as input.
@@ -798,13 +798,13 @@ def zerod_basis(x: np.ndarray) -> np.ndarray:
 
 # {{{ PN bases
 
-def _pkdo_1d(order: Tuple[int], r: np.ndarray) -> np.ndarray:
+def _pkdo_1d(order: tuple[int], r: np.ndarray) -> np.ndarray:
     i, = order
     r0, = r
     return jacobi(0, 0, i, r0)
 
 
-def _grad_pkdo_1d(order: Tuple[int], r: np.ndarray) -> Tuple[np.ndarray]:
+def _grad_pkdo_1d(order: tuple[int], r: np.ndarray) -> tuple[np.ndarray]:
     i, = order
     r0, = r
     return (grad_jacobi(0, 0, i, r0),)
@@ -906,7 +906,7 @@ class TensorProductBasis(Basis):
 
     def __init__(self,
             bases: Sequence[Basis],
-            dims_per_basis: Optional[Tuple[int, ...]] = None) -> None:
+            dims_per_basis: tuple[int, ...] | None = None) -> None:
         """
         :param bases: a sequence of 1D bases used to construct the tensor
             product approximation basis.
@@ -940,7 +940,7 @@ def _nbases(self):
         return len(self._bases)
 
     @property
-    def _mode_index_tuples(self) -> Tuple[Tuple[int, ...], ...]:
+    def _mode_index_tuples(self) -> tuple[tuple[int, ...], ...]:
         from pytools import generate_nonnegative_integer_tuples_below as gnitb
 
         # ensure that these start numbering (0,0), (1,0), (i.e. x-axis first)
@@ -949,7 +949,7 @@ def _mode_index_tuples(self) -> Tuple[Tuple[int, ...], ...]:
                                        for b in self._bases[::-1]]))
 
     @property
-    def mode_ids(self) -> Tuple[Hashable, ...]:
+    def mode_ids(self) -> tuple[Hashable, ...]:
         underlying_mode_id_lists = [basis.mode_ids for basis in self._bases]
         is_all_singletons_with_int = [
                 all(isinstance(mid, tuple) and len(mid) == 1
@@ -957,9 +957,9 @@ def mode_ids(self) -> Tuple[Hashable, ...]:
                     for mid in mode_id_list)
                 for mode_id_list in underlying_mode_id_lists]
 
-        def part_flat_tuple(iterable: Iterable[Tuple[bool, Hashable]]
-                            ) -> Tuple[Hashable, ...]:
-            result: List[Hashable] = []
+        def part_flat_tuple(iterable: Iterable[tuple[bool, Hashable]]
+                            ) -> tuple[Hashable, ...]:
+            result: list[Hashable] = []
             for flatten, item in iterable:
                 if flatten:
                     assert isinstance(item, tuple)
@@ -977,7 +977,7 @@ def part_flat_tuple(iterable: Iterable[Tuple[bool, Hashable]]
                 for mode_index_tuple in self._mode_index_tuples)
 
     @property
-    def functions(self) -> Tuple[BasisFunctionType, ...]:
+    def functions(self) -> tuple[BasisFunctionType, ...]:
         return tuple(
                 _TensorProductBasisFunction(
                     mid,
@@ -987,7 +987,7 @@ def functions(self) -> Tuple[BasisFunctionType, ...]:
                 for mid in self._mode_index_tuples)
 
     @property
-    def gradients(self) -> Tuple[BasisGradientType, ...]:
+    def gradients(self) -> tuple[BasisGradientType, ...]:
         from pytools import wandering_element
         bases = [b.functions for b in self._bases]
         grad_bases = [b.gradients for b in self._bases]
diff --git a/modepy/nodes.py b/modepy/nodes.py
index b9339f72..5548124c 100644
--- a/modepy/nodes.py
+++ b/modepy/nodes.py
@@ -30,6 +30,9 @@
 .. autofunction:: legendre_gauss_lobatto_tensor_product_nodes
 """
 
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2009, 2010, 2013 Andreas Kloeckner, " \
         "Tim Warburton, Jan Hesthaven, Xueyu Zhu"
 
@@ -56,7 +59,7 @@
 # }}}
 
 from functools import partial, singledispatch
-from typing import Optional, Sequence, Tuple, Union
+from typing import Sequence
 
 import numpy as np
 import numpy.linalg as la
@@ -69,7 +72,7 @@
 
 def equidistant_nodes(
         dims: int, n: int,
-        node_tuples: Optional[Sequence[Tuple[int, ...]]] = None
+        node_tuples: Sequence[tuple[int, ...]] | None = None
         ) -> np.ndarray:
     """
     :arg dims: dimensionality of desired simplex
@@ -288,7 +291,7 @@ def warp_and_blend_nodes_3d(n, node_tuples=None):
 
 def warp_and_blend_nodes(
         dims: int, n: int,
-        node_tuples: Optional[Sequence[Tuple[int, ...]]] = None) -> np.ndarray:
+        node_tuples: Sequence[tuple[int, ...]] | None = None) -> np.ndarray:
     """Return interpolation nodes as described in [warburton-nodes]_
 
     .. [warburton-nodes] Warburton, T.
@@ -351,8 +354,8 @@ def warp_and_blend_nodes(
 # {{{ tensor product nodes
 
 def tensor_product_nodes(
-        dims_or_nodes: Union[int, Sequence[np.ndarray]],
-        nodes_1d: Optional[np.ndarray] = None) -> np.ndarray:
+        dims_or_nodes: int | Sequence[np.ndarray],
+        nodes_1d: np.ndarray | None = None) -> np.ndarray:
     """
     :returns: an array of shape ``(dims, nnodes_1d**dims)``.
 
@@ -428,7 +431,7 @@ def legendre_gauss_lobatto_tensor_product_nodes(dims: int, n: int) -> np.ndarray
 # {{{ space-based interface
 
 @singledispatch
-def node_tuples_for_space(space: FunctionSpace) -> Sequence[Tuple[int]]:
+def node_tuples_for_space(space: FunctionSpace) -> Sequence[tuple[int]]:
     raise NotImplementedError(type(space).__name__)
 
 
@@ -445,7 +448,7 @@ def edge_clustered_nodes_for_space(space: FunctionSpace, shape: Shape) -> np.nda
 @singledispatch
 def random_nodes_for_shape(
         shape: Shape, nnodes: int,
-        rng: Optional[np.random.Generator] = None) -> np.ndarray:
+        rng: np.random.Generator | None = None) -> np.ndarray:
     """
     :arg rng: a :class:`numpy.random.Generator`.
 
diff --git a/modepy/quadrature/__init__.py b/modepy/quadrature/__init__.py
index db74db0e..29143123 100644
--- a/modepy/quadrature/__init__.py
+++ b/modepy/quadrature/__init__.py
@@ -51,7 +51,8 @@
 """
 
 from functools import singledispatch
-from typing import Callable, Iterable, Optional, Sequence
+from numbers import Number
+from typing import Callable, Iterable, Sequence
 
 import numpy as np
 
@@ -226,7 +227,7 @@ class LegendreGaussTensorProductQuadrature(TensorProductQuadrature):
 
     def __init__(self,
                  N: int, dims: int,  # noqa: N803
-                 backend: Optional[str] = None) -> None:
+                 backend: str | None = None) -> None:
         from modepy.quadrature.jacobi_gauss import LegendreGaussQuadrature
         super().__init__([
             LegendreGaussQuadrature(N, backend=backend, force_dim_axis=True)
@@ -240,7 +241,7 @@ class LegendreGaussLobattoTensorProductQuadrature(TensorProductQuadrature):
 
     def __init__(self,
                  N: int, dims: int,  # noqa: N803
-                 backend: Optional[str] = None) -> None:
+                 backend: str | None = None) -> None:
         from modepy.quadrature.jacobi_gauss import LegendreGaussLobattoQuadrature
         super().__init__([
             LegendreGaussLobattoQuadrature(N, backend=backend, force_dim_axis=True)
diff --git a/modepy/quadrature/clenshaw_curtis.py b/modepy/quadrature/clenshaw_curtis.py
index 1471afb7..80bdd118 100644
--- a/modepy/quadrature/clenshaw_curtis.py
+++ b/modepy/quadrature/clenshaw_curtis.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2019 Xiaoyu Wei"
 
 __license__ = """
@@ -20,7 +23,6 @@
 THE SOFTWARE.
 """
 
-from typing import Tuple
 
 import numpy as np
 
@@ -29,7 +31,7 @@
 
 # {{{ Clenshaw-Curtis
 
-def _make_clenshaw_curtis_nodes_and_weights(n: int) -> Tuple[np.ndarray, np.ndarray]:
+def _make_clenshaw_curtis_nodes_and_weights(n: int) -> tuple[np.ndarray, np.ndarray]:
     """Nodes and weights of the Clenshaw-Curtis quadrature."""
     if n < 1:
         raise ValueError(f"Clenshaw-Curtis order must be at least 1: n = {n}")
@@ -94,7 +96,7 @@ def __init__(self, N: int, force_dim_axis: bool = False) -> None:  # noqa: N803
 
 # {{{ Fejer
 
-def _make_fejer1_nodes_and_weights(n: int) -> Tuple[np.ndarray, np.ndarray]:
+def _make_fejer1_nodes_and_weights(n: int) -> tuple[np.ndarray, np.ndarray]:
     """Nodes and weights of the Fejer quadrature of the first kind."""
     if n < 1:
         raise ValueError(f"Fejer1 order must be at least 1: n = {n}")
@@ -118,7 +120,7 @@ def _make_fejer1_nodes_and_weights(n: int) -> Tuple[np.ndarray, np.ndarray]:
     return x, w.real
 
 
-def _make_fejer2_nodes_and_weights(n: int) -> Tuple[np.ndarray, np.ndarray]:
+def _make_fejer2_nodes_and_weights(n: int) -> tuple[np.ndarray, np.ndarray]:
     """Nodes and weights of the Fejer quadrature of the second kind."""
     if n < 2:
         raise ValueError(f"Fejer2 order must be at least 2: n = {n}")
diff --git a/modepy/quadrature/grundmann_moeller.py b/modepy/quadrature/grundmann_moeller.py
index fcf62e21..82e50b71 100644
--- a/modepy/quadrature/grundmann_moeller.py
+++ b/modepy/quadrature/grundmann_moeller.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2009-2013 Andreas Kloeckner"
 
 __license__ = """
@@ -21,14 +24,13 @@
 """
 
 from functools import reduce
-from typing import Dict, List, Tuple
 
 import numpy as np
 
 from modepy.quadrature import Quadrature, inf
 
 
-def _extended_euclidean(q: int, r: int) -> Tuple[int, int, int]:
+def _extended_euclidean(q: int, r: int) -> tuple[int, int, int]:
     """Return a tuple (p, a, b) such that p = aq + br,
     where p is the greatest common divisor.
     """
@@ -55,7 +57,7 @@ def _gcd(q: int, r: int) -> int:
     return _extended_euclidean(q, r)[0]
 
 
-def _simplify_fraction(xxx_todo_changeme: Tuple[int, int]) -> Tuple[int, int]:
+def _simplify_fraction(xxx_todo_changeme: tuple[int, int]) -> tuple[int, int]:
     (a, b) = xxx_todo_changeme
     gcd = _gcd(a, b)
     return (a//gcd, b//gcd)
@@ -103,7 +105,7 @@ def __init__(self, order: int, dimension: int) -> None:
             wandering_element,
         )
 
-        points_to_weights: Dict[Tuple[Tuple[int, int], ...], np.ndarray] = {}
+        points_to_weights: dict[tuple[tuple[int, int], ...], np.ndarray] = {}
 
         for i in range(s + 1):
             weight = (-1)**i * 2**(-2*s) \
@@ -127,8 +129,8 @@ def __init__(self, order: int, dimension: int) -> None:
                 + [np.array(x)
                     for x in wandering_element(n, landscape=-1, wanderer=1)])
 
-        nodes: List[np.ndarray] = []
-        weights: List[np.ndarray] = []
+        nodes: list[np.ndarray] = []
+        weights: list[np.ndarray] = []
 
         dim_factor = 2**n
         for p, w in points_to_weights.items():
diff --git a/modepy/quadrature/jacobi_gauss.py b/modepy/quadrature/jacobi_gauss.py
index 69ade7d8..8fd41597 100644
--- a/modepy/quadrature/jacobi_gauss.py
+++ b/modepy/quadrature/jacobi_gauss.py
@@ -25,6 +25,8 @@
 .. autoclass:: JacobiGaussLobattoQuadrature
 .. autoclass:: LegendreGaussLobattoQuadrature
 """
+from __future__ import annotations
+
 
 __copyright__ = "Copyright (C) 2007 Andreas Kloeckner"
 
@@ -48,7 +50,6 @@
 THE SOFTWARE.
 """
 
-from typing import Optional, Tuple
 
 import numpy as np
 import numpy.linalg as la
@@ -73,7 +74,7 @@ class JacobiGaussQuadrature(Quadrature):
 
     def __init__(self,
             alpha: float, beta: float, N: int,                  # noqa: N803
-            backend: Optional[str] = None,
+            backend: str | None = None,
             force_dim_axis: bool = False) -> None:
         r"""
         :arg backend: Either ``"builtin"`` or ``"scipy"``. When the
@@ -119,7 +120,7 @@ def __init__(self,
     @staticmethod
     def compute_weights_and_nodes(
             N: int, alpha: float, beta: float,  # noqa: N803
-            ) -> Tuple[np.ndarray, np.ndarray]:
+            ) -> tuple[np.ndarray, np.ndarray]:
         """
         :arg N: order of the Gauss quadrature (the order of exactly
             integrated polynomials is :math:`2 N + 1`).
@@ -203,7 +204,7 @@ class LegendreGaussQuadrature(JacobiGaussQuadrature):
 
     def __init__(self,
                  N: int, backend:  # noqa: N803
-                 Optional[str] = None,
+                 str | None = None,
                  force_dim_axis: bool = False) -> None:
         super().__init__(
                 0, 0, N,
@@ -225,7 +226,7 @@ class ChebyshevGaussQuadrature(JacobiGaussQuadrature):
     def __init__(self,
                  N: int,  # noqa: N803
                  kind: int = 1,
-                 backend: Optional[str] = None,
+                 backend: str | None = None,
                  force_dim_axis: bool = False) -> None:
         if kind == 1:
             alpha = beta = -0.5
@@ -248,7 +249,7 @@ class GaussGegenbauerQuadrature(JacobiGaussQuadrature):
 
     def __init__(self,
                  alpha: float, N: int,  # noqa: N803
-                 backend: Optional[str] = None,
+                 backend: str | None = None,
                  force_dim_axis: bool = False) -> None:
         super().__init__(
                 alpha, alpha, N,
@@ -257,7 +258,7 @@ def __init__(self,
 
 def jacobi_gauss_lobatto_nodes(
         alpha: float, beta: float, N: int,          # noqa: N803
-        backend: Optional[str] = None,
+        backend: str | None = None,
         force_dim_axis: bool = False) -> np.ndarray:
     """Compute the Gauss-Lobatto quadrature
     nodes corresponding to the :class:`~modepy.JacobiGaussQuadrature`
@@ -287,7 +288,7 @@ def jacobi_gauss_lobatto_nodes(
 
 def legendre_gauss_lobatto_nodes(
         N: int,                     # noqa: N803
-        backend: Optional[str] = None,
+        backend: str | None = None,
         force_dim_axis: bool = False) -> np.ndarray:
     """Compute the Legendre-Gauss-Lobatto quadrature nodes.
     *N+1* is the number of nodes.
@@ -314,7 +315,7 @@ class JacobiGaussLobattoQuadrature(Quadrature):
 
     def __init__(self,
             alpha: float, beta: float, N: int,  # noqa: N803
-            *, backend: Optional[str] = None,
+            *, backend: str | None = None,
             force_dim_axis: bool = False) -> None:
         nodes = jacobi_gauss_lobatto_nodes(alpha, beta, N, backend)
 
@@ -381,7 +382,7 @@ def __init__(self,
 
 class LegendreGaussLobattoQuadrature(JacobiGaussLobattoQuadrature):
     def __init__(
-                self, N, *, backend: Optional[str] = None,  # noqa: N803
+                self, N, *, backend: str | None = None,  # noqa: N803
                 force_dim_axis: bool = False
             ) -> None:
         super().__init__(0, 0, N, backend=backend, force_dim_axis=force_dim_axis)
diff --git a/modepy/quadrature/jaskowiec_sukumar.py b/modepy/quadrature/jaskowiec_sukumar.py
index 299617cf..144e827c 100644
--- a/modepy/quadrature/jaskowiec_sukumar.py
+++ b/modepy/quadrature/jaskowiec_sukumar.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2024 Alexandru Fikl"
 
 __license__ = """
diff --git a/modepy/quadrature/jaskowiec_sukumar_tet_data.py b/modepy/quadrature/jaskowiec_sukumar_tet_data.py
index 61a45f37..820f5b09 100644
--- a/modepy/quadrature/jaskowiec_sukumar_tet_data.py
+++ b/modepy/quadrature/jaskowiec_sukumar_tet_data.py
@@ -1,5 +1,7 @@
 # GENERATED by modepy/contrib/weights-from-jaskowiec-sukumar.py
 # DO NOT EDIT
+from __future__ import annotations
+
 
 __copyright__ = "(C) 2020 Jan Jaśkowiec, N. Sukumar"
 
diff --git a/modepy/quadrature/vioreanu_rokhlin.py b/modepy/quadrature/vioreanu_rokhlin.py
index 4d0c7700..1b7d2749 100644
--- a/modepy/quadrature/vioreanu_rokhlin.py
+++ b/modepy/quadrature/vioreanu_rokhlin.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2009-2013 Andreas Kloeckner"
 
 __license__ = """
diff --git a/modepy/quadrature/vr_quad_data_tet.py b/modepy/quadrature/vr_quad_data_tet.py
index b0966f73..8fe16393 100644
--- a/modepy/quadrature/vr_quad_data_tet.py
+++ b/modepy/quadrature/vr_quad_data_tet.py
@@ -1,5 +1,6 @@
 # GENERATED by modepy/contrib/weights-from-fortran.py
 # DO NOT EDIT
+from __future__ import annotations
 
 import numpy as np
 
diff --git a/modepy/quadrature/vr_quad_data_tri.py b/modepy/quadrature/vr_quad_data_tri.py
index 5f482705..6996c9ef 100644
--- a/modepy/quadrature/vr_quad_data_tri.py
+++ b/modepy/quadrature/vr_quad_data_tri.py
@@ -1,5 +1,6 @@
 # GENERATED by modepy/contrib/weights-from-fortran.py
 # DO NOT EDIT
+from __future__ import annotations
 
 import numpy as np
 
diff --git a/modepy/quadrature/witherden_vincent.py b/modepy/quadrature/witherden_vincent.py
index 1d9d2a88..add2d296 100644
--- a/modepy/quadrature/witherden_vincent.py
+++ b/modepy/quadrature/witherden_vincent.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2020 Alexandru Fikl"
 
 __license__ = """
diff --git a/modepy/quadrature/witherden_vincent_quad_data.py b/modepy/quadrature/witherden_vincent_quad_data.py
index acb439b9..79225fc4 100644
--- a/modepy/quadrature/witherden_vincent_quad_data.py
+++ b/modepy/quadrature/witherden_vincent_quad_data.py
@@ -1,5 +1,7 @@
 # GENERATED by modepy/contrib/weights-from-witherden-vincent.py
 # DO NOT EDIT
+from __future__ import annotations
+
 
 __copyright__ = "(C) 2015 F. D. Witherden, P. E. Vincent"
 
diff --git a/modepy/quadrature/xg_quad_data.py b/modepy/quadrature/xg_quad_data.py
index f052989f..d67c6e0a 100644
--- a/modepy/quadrature/xg_quad_data.py
+++ b/modepy/quadrature/xg_quad_data.py
@@ -1,6 +1,8 @@
 # GENERATED, DO NOT EDIT
 # Xiao-Gimbutas quadratures
 # http://dx.doi.org/10.1016/j.camwa.2009.10.027
+from __future__ import annotations
+
 import numpy
 
 
diff --git a/modepy/quadrature/xiao_gimbutas.py b/modepy/quadrature/xiao_gimbutas.py
index f4b81300..ac1bdad4 100644
--- a/modepy/quadrature/xiao_gimbutas.py
+++ b/modepy/quadrature/xiao_gimbutas.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2009-2013 Andreas Kloeckner"
 
 __license__ = """
diff --git a/modepy/shapes.py b/modepy/shapes.py
index 7d833fde..63fa37f9 100644
--- a/modepy/shapes.py
+++ b/modepy/shapes.py
@@ -187,6 +187,8 @@
 """
 
 # }}}
+from __future__ import annotations
+
 
 __copyright__ = """
 Copyright (c) 2013 Andreas Kloeckner
@@ -217,7 +219,7 @@
 from abc import ABC, abstractmethod
 from dataclasses import dataclass
 from functools import partial, singledispatch
-from typing import Any, Callable, Dict, Sequence, Tuple
+from typing import Any, Callable, Sequence
 
 import numpy as np
 
@@ -258,7 +260,7 @@ class Face:
     """The volume :class:`Shape` to which the face belongs."""
     face_index: int
     """The face index in :attr:`volume_shape` of this face."""
-    volume_vertex_indices: Tuple[int, ...]
+    volume_vertex_indices: tuple[int, ...]
     """A tuple of indices into the vertices returned by
     :func:`unit_vertices_for_shape` for the :attr:`volume_shape`.
     """
@@ -301,7 +303,7 @@ def face_normal(face: Face, normalize: bool = True) -> np.ndarray:
 
 
 @singledispatch
-def faces_for_shape(shape: Shape) -> Tuple[Face, ...]:
+def faces_for_shape(shape: Shape) -> tuple[Face, ...]:
     r"""
     :returns: a tuple of :class:`Face`\ s representing the faces of *shape*.
     """
@@ -320,17 +322,17 @@ class TensorProductShape(Shape):
     results in a prism shape. Special cases include the :class:`Hypercube`.
     """
 
-    bases: Tuple[Shape, ...]
+    bases: tuple[Shape, ...]
     """A :class:`tuple` of base shapes that form the tensor product."""
 
     # NOTE: https://github.com/python/mypy/issues/1020
-    def __new__(cls, bases: Tuple[Shape, ...]) -> Any:
+    def __new__(cls, bases: tuple[Shape, ...]) -> Any:
         if len(bases) == 1:
             return bases[0]
         else:
             return Shape.__new__(cls)
 
-    def __init__(self, bases: Tuple[Shape, ...]) -> None:
+    def __init__(self, bases: tuple[Shape, ...]) -> None:
         # flatten input shapes
         bases = sum((
             s.bases if isinstance(s, TensorProductShape) else (s,)
@@ -411,7 +413,7 @@ def _simplex_face_to_vol_map(face_vertices, p: np.ndarray):
     return origin + np.einsum("ij,jk->ik", face_basis, (1 + p) / 2)
 
 
-_SIMPLEX_FACES: Dict[int, Tuple[Tuple[int, ...], ...]] = {
+_SIMPLEX_FACES: dict[int, tuple[tuple[int, ...], ...]] = {
             1: ((0,), (1,)),
             2: ((0, 1), (2, 0), (1, 2)),
             3: ((0, 2, 1), (0, 1, 3), (0, 3, 2), (1, 2, 3))
@@ -488,7 +490,7 @@ def _hypercube_face_to_vol_map(face_vertices: np.ndarray, p: np.ndarray):
     return origin + np.einsum("ij,jk->ik", face_basis, (1 + p) / 2)
 
 
-_HYPERCUBE_FACES: Dict[int, Tuple[Tuple[int, ...], ...]] = {
+_HYPERCUBE_FACES: dict[int, tuple[tuple[int, ...], ...]] = {
         1: ((0b0,), (0b1,)),
         2: ((0b00, 0b01), (0b11, 0b10), (0b10, 0b00), (0b01, 0b11)),
         3: (
@@ -525,8 +527,8 @@ def _faces_for_hypercube(shape: Hypercube):
 
 @singledispatch
 def submesh_for_shape(
-        shape: Shape, node_tuples: Sequence[Tuple[int, ...]]
-        ) -> Sequence[Tuple[int, ...]]:
+        shape: Shape, node_tuples: Sequence[tuple[int, ...]]
+        ) -> Sequence[tuple[int, ...]]:
     """Return a list of tuples of indices into the node list that
     generate a tessellation of the reference element.
 
diff --git a/modepy/spaces.py b/modepy/spaces.py
index b87994ac..9adb4b27 100644
--- a/modepy/spaces.py
+++ b/modepy/spaces.py
@@ -21,6 +21,8 @@
 
 .. autofunction:: space_for_shape
 """
+from __future__ import annotations
+
 
 __copyright__ = "Copyright (C) 2020 Andreas Kloeckner"
 
@@ -47,7 +49,7 @@
 from abc import ABC, abstractmethod
 from functools import singledispatch
 from numbers import Number
-from typing import Any, Tuple, Union
+from typing import Any
 
 import numpy as np
 
@@ -81,7 +83,7 @@ def space_dim(self) -> int:
 
 @singledispatch
 def space_for_shape(
-        shape: Shape, order: Union[int, Tuple[int, ...]]
+        shape: Shape, order: int | tuple[int, ...]
         ) -> FunctionSpace:
     r"""Return an unspecified instance of :class:`FunctionSpace` suitable
     for approximation on *shape* attaining interpolation error of
@@ -108,17 +110,17 @@ class TensorProductSpace(FunctionSpace):
     :func:`~modepy.tools.reshape_array_for_tensor_product_space`.
     """
 
-    bases: Tuple[FunctionSpace, ...]
+    bases: tuple[FunctionSpace, ...]
     """A :class:`tuple` of the base spaces that take part in the tensor product."""
 
     # NOTE: https://github.com/python/mypy/issues/1020
-    def __new__(cls, bases: Tuple[FunctionSpace, ...]) -> Any:
+    def __new__(cls, bases: tuple[FunctionSpace, ...]) -> Any:
         if len(bases) == 1:
             return bases[0]
         else:
             return FunctionSpace.__new__(cls)
 
-    def __init__(self, bases: Tuple[FunctionSpace, ...]) -> None:
+    def __init__(self, bases: tuple[FunctionSpace, ...]) -> None:
         self.bases = sum((
             space.bases if isinstance(space, TensorProductSpace) else (space,)
             for space in bases
@@ -156,7 +158,7 @@ def __repr__(self) -> str:
 @space_for_shape.register(TensorProductShape)
 def _space_for_tensor_product_shape(
         shape: TensorProductShape,
-        order: Union[int, Tuple[int, ...]]) -> TensorProductSpace:
+        order: int | tuple[int, ...]) -> TensorProductSpace:
     nbases = len(shape.bases)
     if isinstance(order, Number):
         assert isinstance(order, int)
@@ -213,7 +215,7 @@ def __repr__(self) -> str:
 
 
 @space_for_shape.register(Simplex)
-def _space_for_simplex(shape: Simplex, order: Union[int, Tuple[int, ...]]) -> PN:
+def _space_for_simplex(shape: Simplex, order: int | tuple[int, ...]) -> PN:
     assert isinstance(order, int)
     return PN(shape.dim, order)
 
diff --git a/modepy/test/test_matrices.py b/modepy/test/test_matrices.py
index 6b7fcc6a..8d8a2e9c 100644
--- a/modepy/test/test_matrices.py
+++ b/modepy/test/test_matrices.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2024 University of Illinois Board of Trustees"
 
 __license__ = """
diff --git a/modepy/test/test_modes.py b/modepy/test/test_modes.py
index 6dc1c109..4d7b5a2b 100644
--- a/modepy/test/test_modes.py
+++ b/modepy/test/test_modes.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2009-2013 Andreas Kloeckner"
 
 __license__ = """
diff --git a/modepy/test/test_nodes.py b/modepy/test/test_nodes.py
index 1306e4d8..71de1cc3 100644
--- a/modepy/test/test_nodes.py
+++ b/modepy/test/test_nodes.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2009-2013 Andreas Kloeckner"
 
 __license__ = """
diff --git a/modepy/test/test_quadrature.py b/modepy/test/test_quadrature.py
index 287efa21..cf6643ef 100644
--- a/modepy/test/test_quadrature.py
+++ b/modepy/test/test_quadrature.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2009-2013 Andreas Kloeckner"
 
 __license__ = """
diff --git a/modepy/test/test_tools.py b/modepy/test/test_tools.py
index b73f216b..0f2d27b3 100644
--- a/modepy/test/test_tools.py
+++ b/modepy/test/test_tools.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2013 Andreas Kloeckner"
 
 __license__ = """
diff --git a/modepy/tools.py b/modepy/tools.py
index 986c24b5..9f8c5765 100644
--- a/modepy/tools.py
+++ b/modepy/tools.py
@@ -31,6 +31,8 @@
     .. method::: reshape
 
 """
+from __future__ import annotations
+
 
 __copyright__ = "Copyright (C) 2013 Andreas Kloeckner"
 
@@ -59,7 +61,7 @@
     gamma,  # noqa: F401
     sqrt,
 )
-from typing import Dict, Optional, Protocol, Tuple, TypeVar, runtime_checkable
+from typing import Protocol, TypeVar, runtime_checkable
 from warnings import warn
 
 import numpy as np
@@ -172,7 +174,7 @@ def inverse(self):
 
 # {{{ simplex coordinate mapping
 
-EQUILATERAL_TO_UNIT_MAP: Dict[int, AffineMap] = {
+EQUILATERAL_TO_UNIT_MAP: dict[int, AffineMap] = {
         1: AffineMap([[1]], [0]),
         2: AffineMap([
             [1, -1/sqrt(3)],
@@ -294,8 +296,8 @@ def plot_element_values(
         n: int,
         nodes: np.ndarray,
         values: np.ndarray,
-        resample_n: Optional[int] = None,
-        node_tuples: Optional[Tuple[int, ...]] = None,
+        resample_n: int | None = None,
+        node_tuples: tuple[int, ...] | None = None,
         show_nodes: bool = False) -> None:
     """
     :arg n: order of the polynomial basis.
@@ -384,7 +386,7 @@ def _evaluate_lebesgue_function(n, nodes, shape):
 def estimate_lebesgue_constant(
         n: int,
         nodes: np.ndarray,
-        shape: Optional[shp.Shape] = None,
+        shape: shp.Shape | None = None,
         *,
         visualize: bool = False) -> float:
     """Estimate the
@@ -461,11 +463,11 @@ def estimate_lebesgue_constant(
 
 @runtime_checkable
 class Reshapeable(Protocol):
-    shape: Tuple[int, ...]
+    shape: tuple[int, ...]
 
     def reshape(
-            self: "ReshapeableT", *newshape: Tuple[int, ...], order: str
-            ) -> "ReshapeableT":
+            self: ReshapeableT, *newshape: tuple[int, ...], order: str
+            ) -> ReshapeableT:
         ...
 
 
diff --git a/modepy/version.py b/modepy/version.py
index e7d6757c..2988bf95 100644
--- a/modepy/version.py
+++ b/modepy/version.py
@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import re
 from importlib import metadata
 
diff --git a/pyproject.toml b/pyproject.toml
index 64e0bc52..c9489dcc 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -115,6 +115,7 @@ known-local-folder = [
     "modepy",
 ]
 lines-after-imports = 2
+required-imports = ["from __future__ import annotations"]
 
 [tool.mypy]
 warn_unused_ignores = true