diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 027e6533..8e5ab234 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -25,6 +25,23 @@ jobs:
         -   uses: actions/checkout@v4
         -   uses: crate-ci/typos@master
 
+    mypy:
+        name: Mypy
+        runs-on: ubuntu-latest
+        steps:
+        -   uses: actions/checkout@v4
+        -
+            uses: actions/setup-python@v5
+            with:
+                python-version: '3.x'
+        -   name: "Main Script"
+            run: |
+                curl -L -O https://tiker.net/ci-support-v0
+                . ./ci-support-v0
+                build_py_project_in_conda_env
+                python -m pip install mypy
+                mypy $(get_proj_name)
+
     pylint:
         name: Pylint
         runs-on: ubuntu-latest
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 2326e653..c494bc44 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -140,6 +140,18 @@ Pylint:
   except:
   - tags
 
+Mypy:
+  script: |
+    curl -L -O https://tiker.net/ci-support-v0
+    . ./ci-support-v0
+    build_py_project_in_venv
+    python -m pip install mypy
+    mypy $(get_proj_name)
+  tags:
+  - python3
+  except:
+  - tags
+
 Downstream:
   parallel:
     matrix:
diff --git a/benchmarks/bench_translations.py b/benchmarks/bench_translations.py
index 8cd52af4..8c0e1772 100644
--- a/benchmarks/bench_translations.py
+++ b/benchmarks/bench_translations.py
@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import logging
 
 import numpy as np
diff --git a/pyproject.toml b/pyproject.toml
index bbdbb5e7..d6753eff 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -42,6 +42,12 @@ known-local-folder = [
     "sumpy",
 ]
 lines-after-imports = 2
+required-imports = ["from __future__ import annotations"]
+
+[tool.ruff.lint.per-file-ignores]
+"doc/**/*.py" = ["I002"]
+"examples/**/*.py" = ["I002"]
+"setup.py" = ["I002"]
 
 [tool.typos.default]
 extend-ignore-re = [
@@ -61,3 +67,20 @@ extend-exclude = [
   "contrib/*/*.ipynb",
   "notes/*/*.eps",
 ]
+
+[tool.mypy]
+warn_unused_ignores = true
+
+[[tool.mypy.overrides]]
+module = [
+    "pyopencl.*",
+    "pymbolic.*",
+    "loopy.*",
+    "symengine.*",
+    "boxtree.*",
+    "pyvkfft.*",
+    "pyfmmlib.*",
+    "mayavi.*",
+    "scipy.*",
+]
+ignore_missing_imports = true
diff --git a/sumpy/__init__.py b/sumpy/__init__.py
index 29c5904b..ab98d104 100644
--- a/sumpy/__init__.py
+++ b/sumpy/__init__.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2013 Andreas Kloeckner"
 
 __license__ = """
@@ -21,7 +24,9 @@
 """
 
 import os
+from typing import Hashable
 
+import loopy as lp
 from pytools.persistent_dict import WriteOncePersistentDict
 
 from sumpy.e2e import (
@@ -56,8 +61,8 @@
 ]
 
 
-code_cache = WriteOncePersistentDict("sumpy-code-cache-v6-"+VERSION_TEXT,
-                                     safe_sync=False)
+code_cache: WriteOncePersistentDict[Hashable, lp.TranslationUnit] = \
+        WriteOncePersistentDict("sumpy-code-cache-v6-"+VERSION_TEXT, safe_sync=False)
 
 
 # {{{ optimization control
diff --git a/sumpy/array_context.py b/sumpy/array_context.py
index 21432af3..eddfb2cc 100644
--- a/sumpy/array_context.py
+++ b/sumpy/array_context.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2022 Alexandru Fikl"
 
 __license__ = """
diff --git a/sumpy/assignment_collection.py b/sumpy/assignment_collection.py
index d4906455..678dbbd4 100644
--- a/sumpy/assignment_collection.py
+++ b/sumpy/assignment_collection.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2012 Andreas Kloeckner"
 
 __license__ = """
diff --git a/sumpy/codegen.py b/sumpy/codegen.py
index 139a17b0..22ee9d64 100644
--- a/sumpy/codegen.py
+++ b/sumpy/codegen.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2012 Andreas Kloeckner"
 
 __license__ = """
diff --git a/sumpy/cse.py b/sumpy/cse.py
index 35b39980..6dbb4d47 100644
--- a/sumpy/cse.py
+++ b/sumpy/cse.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = """
 Copyright (C) 2017 Matt Wala
 Copyright (C) 2006-2016 SymPy Development Team
diff --git a/sumpy/derivative_taker.py b/sumpy/derivative_taker.py
index 7e619202..ffe23485 100644
--- a/sumpy/derivative_taker.py
+++ b/sumpy/derivative_taker.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = """
 Copyright (C) 2012 Andreas Kloeckner
 Copyright (C) 2020 Isuru Fernando
@@ -338,7 +341,7 @@ def diff(self, mi, q=0):
 
 # {{{ DifferentiatedExprDerivativeTaker
 
-DerivativeCoeffDict = Dict[Tuple[int], Any]
+DerivativeCoeffDict = Dict[Tuple[int, ...], Any]
 
 
 @tag_dataclass
@@ -390,7 +393,7 @@ def diff_derivative_coeff_dict(derivative_coeff_dict: DerivativeCoeffDict,
     and return a new derivative transformation dictionary.
     """
     from collections import defaultdict
-    new_derivative_coeff_dict = defaultdict(lambda: 0)
+    new_derivative_coeff_dict: DerivativeCoeffDict = defaultdict(lambda: 0)
 
     for mi, coeff in derivative_coeff_dict.items():
         # In the case where we have x * u.diff(x), the result should
diff --git a/sumpy/distributed.py b/sumpy/distributed.py
index e0b6ee3b..84f3c4e8 100644
--- a/sumpy/distributed.py
+++ b/sumpy/distributed.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2022 Hao Gao"
 
 __license__ = """
diff --git a/sumpy/e2e.py b/sumpy/e2e.py
index dde5d273..e6bedebf 100644
--- a/sumpy/e2e.py
+++ b/sumpy/e2e.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2013 Andreas Kloeckner"
 
 __license__ = """
diff --git a/sumpy/e2p.py b/sumpy/e2p.py
index 07066f39..87e2f139 100644
--- a/sumpy/e2p.py
+++ b/sumpy/e2p.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2013 Andreas Kloeckner"
 
 __license__ = """
diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py
index e2edd47c..6a4bf3a4 100644
--- a/sumpy/expansion/__init__.py
+++ b/sumpy/expansion/__init__.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2012 Andreas Kloeckner"
 
 __license__ = """
@@ -22,7 +25,7 @@
 
 import logging
 from abc import ABC, abstractmethod
-from typing import Any, ClassVar, Dict, Hashable, List, Optional, Sequence, Tuple, Type
+from typing import Any, ClassVar, Hashable, Sequence
 
 import loopy as lp
 import pymbolic.primitives as prim
@@ -77,12 +80,12 @@ class ExpansionBase(ABC):
     .. automethod:: __eq__
     .. automethod:: __ne__
     """
-    init_arg_names = ("kernel", "order", "use_rscale")
+    init_arg_names: tuple[str, ...] = ("kernel", "order", "use_rscale")
 
     def __init__(self,
             kernel: Kernel,
             order: int,
-            use_rscale: Optional[bool] = None) -> None:
+            use_rscale: bool | None = None) -> None:
         # Don't be tempted to remove target derivatives here.
         # Line Taylor QBX can't do without them, because it can't
         # apply those derivatives to the expanded quantity.
@@ -125,7 +128,7 @@ def get_source_args(self):
     # {{{ abstract interface
 
     @abstractmethod
-    def get_coefficient_identifiers(self) -> List[Hashable]:
+    def get_coefficient_identifiers(self) -> list[Hashable]:
         """
         :returns: the identifiers of the coefficients that actually get stored.
         """
@@ -195,10 +198,10 @@ def loopy_evaluator(self, kernels: Sequence[Kernel]) -> lp.TranslationUnit:
 
     # {{{ copy
 
-    def with_kernel(self, kernel: Kernel) -> "ExpansionBase":
+    def with_kernel(self, kernel: Kernel) -> ExpansionBase:
         return type(self)(kernel, self.order, self.use_rscale)
 
-    def copy(self, **kwargs) -> "ExpansionBase":
+    def copy(self, **kwargs) -> ExpansionBase:
         new_kwargs = {
                 name: getattr(self, name)
                 for name in self.init_arg_names}
@@ -250,12 +253,12 @@ class ExpansionTermsWrangler(ABC):
 
     .. automethod:: get_full_coefficient_identifiers
     """
-    init_arg_names = ("order", "dim", "max_mi")
+    init_arg_names: tuple[str, ...] = ("order", "dim", "max_mi")
 
     def __init__(self,
             order: int,
             dim: int,
-            max_mi: Optional[int] = None) -> None:
+            max_mi: tuple[int, ...] | None = None) -> None:
         self.order = order
         self.dim = dim
         self.max_mi = max_mi
@@ -263,7 +266,7 @@ def __init__(self,
     # {{{ abstract interface
 
     @abstractmethod
-    def get_coefficient_identifiers(self) -> List[Hashable]:
+    def get_coefficient_identifiers(self) -> list[tuple[int, ...]]:
         pass
 
     @abstractmethod
@@ -279,7 +282,7 @@ def get_stored_mpole_coefficients_from_full(self,
     # }}}
 
     @memoize_method
-    def get_full_coefficient_identifiers(self) -> List[Hashable]:
+    def get_full_coefficient_identifiers(self) -> list[Hashable]:
         """
         Returns identifiers for every coefficient in the complete expansion.
         """
@@ -311,7 +314,7 @@ def copy(self, **kwargs):
 
     # {{{ hyperplane helpers
 
-    def _get_mi_hyperpplanes(self) -> List[Tuple[int, int]]:
+    def _get_mi_hyperpplanes(self) -> list[tuple[int, int]]:
         r"""
         Coefficient storage is organized into "hyperplanes" in multi-index
         space. Potentially only a subset of these hyperplanes contain
@@ -332,7 +335,7 @@ def _get_mi_hyperpplanes(self) -> List[Tuple[int, int]]:
 
     @memoize_method
     def _split_coeffs_into_hyperplanes(
-            self) -> List[Tuple[int, List[Tuple[int, ...]]]]:
+            self) -> list[tuple[int, list[tuple[int, ...]]]]:
         r"""
         This splits the coefficients into :math:`O(p)` disjoint sets
         so that for each set, all the identifiers have the form,
@@ -511,7 +514,7 @@ class LinearPDEBasedExpansionTermsWrangler(ExpansionTermsWrangler):
 
     def __init__(self,
             order: int, dim: int, knl: Kernel,
-            max_mi: Optional[int] = None) -> None:
+            max_mi: tuple[int, ...] | None = None) -> None:
         r"""
         :param order: order of the expansion
         :param dim: number of dimensions
@@ -592,7 +595,7 @@ def mi_key(ident):
 
         return mi_key, axis_permutation
 
-    def _get_mi_hyperpplanes(self) -> List[Tuple[int, int]]:
+    def _get_mi_hyperpplanes(self) -> list[tuple[int, int]]:
         mis = self.get_full_coefficient_identifiers()
         mi_to_index = {mi: i for i, mi in enumerate(mis)}
 
@@ -804,8 +807,8 @@ def get_projection_matrix(self, rscale):
 # FIXME: This is called an expansion but doesn't inherit from ExpansionBase?
 
 class VolumeTaylorExpansionMixin:
-    expansion_terms_wrangler_class: ClassVar[Type[ExpansionTermsWrangler]]
-    expansion_terms_wrangler_cache: ClassVar[Dict[Hashable, Any]] = {}
+    expansion_terms_wrangler_class: ClassVar[type[ExpansionTermsWrangler]]
+    expansion_terms_wrangler_cache: ClassVar[dict[Hashable, Any]] = {}
 
     @classmethod
     def get_or_make_expansion_terms_wrangler(cls, *key):
diff --git a/sumpy/expansion/diff_op.py b/sumpy/expansion/diff_op.py
index bc439e51..63c234a0 100644
--- a/sumpy/expansion/diff_op.py
+++ b/sumpy/expansion/diff_op.py
@@ -1,3 +1,8 @@
+# mypy: disallow-untyped-defs
+
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2019 Isuru Fernando"
 
 __license__ = """
@@ -23,9 +28,11 @@
 import logging
 from dataclasses import dataclass
 from itertools import accumulate
-from typing import List, Mapping, Sequence
+from typing import Mapping, Sequence, Union
 
+import numpy as np
 import sympy as sp
+import sympy.polys.agca.modules as sp_modules
 from pyrsistent import pmap
 
 from pytools import memoize
@@ -67,6 +74,9 @@ class DerivativeIdentifier:
     """
 
 
+Number_ish = Union[int, float, complex, np.number]
+
+
 @dataclass(frozen=True, eq=True)
 class LinearPDESystemOperator:
     r"""
@@ -86,16 +96,20 @@ class LinearPDESystemOperator:
     """
 
     dim: int
-    eqs: Sequence[Mapping[DerivativeIdentifier, sp.Expr]]
+    eqs: tuple[Mapping[DerivativeIdentifier, sp.Expr], ...]
+
+    if __debug__:
+        def __post_init__(self) -> None:
+            hash(self)
 
     @property
-    def order(self):
+    def order(self) -> int:
         deg = 0
         for eq in self.eqs:
             deg = max(deg, max(sum(ident.mi) for ident in eq.keys()))
         return deg
 
-    def __mul__(self, param):
+    def __mul__(self, param: Number_ish) -> LinearPDESystemOperator:
         eqs = []
         for eq in self.eqs:
             deriv_ident_to_coeff = {}
@@ -106,7 +120,9 @@ def __mul__(self, param):
 
     __rmul__ = __mul__
 
-    def __add__(self, other_diff_op):
+    def __add__(
+                self, other_diff_op: LinearPDESystemOperator
+            ) -> LinearPDESystemOperator:
         assert self.dim == other_diff_op.dim
         assert len(self.eqs) == len(other_diff_op.eqs)
         eqs = []
@@ -122,27 +138,31 @@ def __add__(self, other_diff_op):
 
     __radd__ = __add__
 
-    def __sub__(self, other_diff_op):
+    def __sub__(
+                self, other_diff_op: LinearPDESystemOperator
+            ) -> LinearPDESystemOperator:
         return self + (-1)*other_diff_op
 
-    def __repr__(self):
+    def __repr__(self) -> str:
         return f"LinearPDESystemOperator({self.dim}, {self.eqs!r})"
 
-    def __getitem__(self, idx):
+    def __getitem__(self, idx: int | slice) -> LinearPDESystemOperator:
         item = self.eqs.__getitem__(idx)
-        if not isinstance(item, tuple):
-            item = (item,)
-        return LinearPDESystemOperator(self.dim, tuple(item))
+        if isinstance(item, tuple):
+            eqs = item
+        else:
+            eqs = (item,)
+        return LinearPDESystemOperator(self.dim, eqs)
 
     @property
-    def total_dims(self):
+    def total_dims(self) -> int:
         """
         Returns the total number of dimensions including time
         """
         did = next(iter(self.eqs[0].keys()))
         return len(did.mi)
 
-    def to_sym(self, fnames=None):
+    def to_sym(self, fnames: Sequence[str] | None = None) -> list[sp.Expr]:
         from sumpy.symbolic import Function, make_sym_vector
         x = list(make_sym_vector("x", self.dim))
         x += list(make_sym_vector("t", self.total_dims - self.dim))
@@ -158,7 +178,7 @@ def to_sym(self, fnames=None):
 
         res = []
         for eq in self.eqs:
-            sym_eq = 0
+            sym_eq: sp.Expr = sp.sympify(0)
             for deriv_ident, coeff in eq.items():
                 expr = funcs[deriv_ident.vec_idx]
                 for i, val in enumerate(deriv_ident.mi):
@@ -169,7 +189,10 @@ def to_sym(self, fnames=None):
         return res
 
 
-def convert_module_to_matrix(module, generators):
+def convert_module_to_matrix(
+            module: Sequence[sp_modules.FreeModuleElement],
+            generators: Sequence[sp.Expr]
+        ) -> sp.Matrix:
     import sympy
     # poly is a sympy DMP (dense multi-variate polynomial)
     # type and we convert it to a sympy expression because
@@ -180,8 +203,7 @@ def convert_module_to_matrix(module, generators):
 
 
 @memoize
-def _get_all_scalar_pdes(pde: LinearPDESystemOperator) -> \
-        List[LinearPDESystemOperator]:
+def _get_all_scalar_pdes(pde: LinearPDESystemOperator) -> list[LinearPDESystemOperator]:
     import sympy
     from sympy.polys.orderings import grevlex
     gens = [sympy.symbols(f"_x{i}") for i in range(pde.dim)]
@@ -214,7 +236,10 @@ def _get_all_scalar_pdes(pde: LinearPDESystemOperator) -> \
     # for each column we calculate the intersection of the left modules and the
     # right modules. This requires only $3*(n-2)$ work.
 
-    def intersect(a, b):
+    def intersect(
+                a: sp_modules.SubModulePolyRing,
+                b: sp_modules.SubModulePolyRing,
+            ) -> sp_modules.SubModulePolyRing:
         return a.intersect(b)
 
     left_intersections = list(accumulate(column_syzygy_modules, func=intersect))
@@ -247,7 +272,7 @@ def intersect(a, b):
             DerivativeIdentifier(mi, 0): sym.sympify(coeff.as_expr().simplify()) for
             (mi, coeff) in zip(scalar_pde.monoms(), scalar_pde.coeffs())
         }
-        results.append(LinearPDESystemOperator(pde.dim, pmap(pde_dict)))
+        results.append(LinearPDESystemOperator(pde.dim, (pmap(pde_dict),)))
 
     return results
 
@@ -323,9 +348,10 @@ def as_scalar_pde(pde: LinearPDESystemOperator, comp_idx: int) \
     return _get_all_scalar_pdes(pde)[comp_idx]
 
 
-def laplacian(diff_op):
+def laplacian(diff_op: LinearPDESystemOperator) -> LinearPDESystemOperator:
     dim = diff_op.dim
-    empty = [pmap()] * len(diff_op.eqs)
+    empty: tuple[Mapping[DerivativeIdentifier, sp.Expr], ...] = \
+        (pmap(),) * len(diff_op.eqs)
     res = LinearPDESystemOperator(dim, empty)
     for j in range(dim):
         mi = [0]*diff_op.total_dims
@@ -334,7 +360,9 @@ def laplacian(diff_op):
     return res
 
 
-def diff(diff_op, mi):
+def diff(
+            diff_op: LinearPDESystemOperator, mi: tuple[int, ...]
+        ) -> LinearPDESystemOperator:
     eqs = []
     for eq in diff_op.eqs:
         res = {}
@@ -345,9 +373,9 @@ def diff(diff_op, mi):
     return LinearPDESystemOperator(diff_op.dim, tuple(eqs))
 
 
-def divergence(diff_op):
+def divergence(diff_op: LinearPDESystemOperator) -> LinearPDESystemOperator:
     assert len(diff_op.eqs) == diff_op.dim
-    res = LinearPDESystemOperator(diff_op.dim, pmap())
+    res = LinearPDESystemOperator(diff_op.dim, (pmap(),))
     for i in range(diff_op.dim):
         mi = [0]*diff_op.total_dims
         mi[i] = 1
@@ -355,7 +383,7 @@ def divergence(diff_op):
     return res
 
 
-def gradient(diff_op):
+def gradient(diff_op: LinearPDESystemOperator) -> LinearPDESystemOperator:
     assert len(diff_op.eqs) == 1
     eqs = []
     dim = diff_op.dim
@@ -366,26 +394,25 @@ def gradient(diff_op):
     return LinearPDESystemOperator(dim, tuple(eqs))
 
 
-def curl(pde):
-    assert len(pde.eqs) == 3
-    assert pde.dim == 3
+def curl(diff_op: LinearPDESystemOperator) -> LinearPDESystemOperator:
+    assert len(diff_op.eqs) == 3
+    assert diff_op.dim == 3
     eqs = []
     mis = []
     for i in range(3):
-        mi = [0]*pde.total_dims
+        mi = [0]*diff_op.total_dims
         mi[i] = 1
         mis.append(tuple(mi))
 
     for i in range(3):
-        new_pde = diff(pde[(i+2) % 3], mis[(i+1) % 3]) - \
-            diff(pde[(i+1) % 3], mis[(i+2) % 3])
+        new_pde = diff(diff_op[(i+2) % 3], mis[(i+1) % 3]) - \
+            diff(diff_op[(i+1) % 3], mis[(i+2) % 3])
         eqs.append(new_pde.eqs[0])
 
-    return LinearPDESystemOperator(pde.dim, tuple(eqs))
+    return LinearPDESystemOperator(diff_op.dim, tuple(eqs))
 
 
-def concat(*ops):
-    ops = list(ops)
+def concat(*ops: LinearPDESystemOperator) -> LinearPDESystemOperator:
     assert len(ops) >= 1
     dim = ops[0].dim
     for op in ops:
@@ -393,10 +420,12 @@ def concat(*ops):
     eqs = list(ops[0].eqs)
     for op in ops[1:]:
         eqs.extend(list(op.eqs))
-    return LinearPDESystemOperator(dim, eqs)
+    return LinearPDESystemOperator(dim, tuple(eqs))
 
 
-def make_identity_diff_op(ninput, noutput=1, time_dependent=False):
+def make_identity_diff_op(
+            ninput: int, noutput: int = 1, time_dependent: bool = False
+        ) -> LinearPDESystemOperator:
     """
     Returns the identity as a linear PDE system operator.
     if *include_time* is true, then the last dimension of the
@@ -410,5 +439,7 @@ def make_identity_diff_op(ninput, noutput=1, time_dependent=False):
         mi = tuple([0]*(ninput + 1))
     else:
         mi = tuple([0]*ninput)
-    eqs = [pmap({DerivativeIdentifier(mi, i): 1}) for i in range(noutput)]
+    eqs = tuple(pmap(
+                    {DerivativeIdentifier(mi, i): sp.sympify(1)})
+                    for i in range(noutput))
     return LinearPDESystemOperator(ninput, eqs)
diff --git a/sumpy/expansion/level_to_order.py b/sumpy/expansion/level_to_order.py
index 37a729b9..c81a4bdc 100644
--- a/sumpy/expansion/level_to_order.py
+++ b/sumpy/expansion/level_to_order.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2016 Matt Wala"
 
 __license__ = """
diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py
index 512ba782..712b26e3 100644
--- a/sumpy/expansion/local.py
+++ b/sumpy/expansion/local.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2012 Andreas Kloeckner"
 
 __license__ = """
diff --git a/sumpy/expansion/loopy.py b/sumpy/expansion/loopy.py
index 6c4ad370..22ee5d6a 100644
--- a/sumpy/expansion/loopy.py
+++ b/sumpy/expansion/loopy.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2022 Isuru Fernando"
 
 __license__ = """
diff --git a/sumpy/expansion/m2l.py b/sumpy/expansion/m2l.py
index 2dc5f8e9..134ed0c0 100644
--- a/sumpy/expansion/m2l.py
+++ b/sumpy/expansion/m2l.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2022 Isuru Fernando"
 
 __license__ = """
@@ -22,7 +25,7 @@
 
 import logging
 from abc import ABC, abstractmethod
-from typing import Any, ClassVar, Tuple
+from typing import Any, ClassVar
 
 import numpy as np
 
@@ -164,7 +167,7 @@ def loopy_translate(self, tgt_expansion, src_expansion):
             f"{src_expansion} to {tgt_expansion} using {self} is not implemented.")
 
     def translation_classes_dependent_data(self, tgt_expansion, src_expansion,
-            src_rscale, dvec, sac) -> Tuple[Any]:
+            src_rscale, dvec, sac) -> tuple[Any, ...]:
         """Return an iterable of expressions that needs to be precomputed
         for multipole-to-local translations that depend only on the
         distance between the multipole center and the local center which
diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py
index eed92b3b..d9a81764 100644
--- a/sumpy/expansion/multipole.py
+++ b/sumpy/expansion/multipole.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2012 Andreas Kloeckner"
 
 __license__ = """
diff --git a/sumpy/fmm.py b/sumpy/fmm.py
index 836182f3..aa2ff4c7 100644
--- a/sumpy/fmm.py
+++ b/sumpy/fmm.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2013 Andreas Kloeckner"
 
 __license__ = """
@@ -27,7 +30,7 @@
 """
 
 
-from typing import List, TypeVar, Union
+from typing import Protocol
 
 from boxtree.fmm import ExpansionWranglerInterface, TreeIndependentDataForWrangler
 
@@ -203,11 +206,12 @@ def opencl_fft_app(self, shape, dtype, inverse):
 _SECONDS_PER_NANOSECOND = 1e-9
 
 
-"""
-EventLike objects have an attribute native_event that returns
-a cl.Event that indicates the end of the event.
-"""
-EventLike = TypeVar("CLEventLike")
+class CLEventLike(Protocol):
+    """
+    EventLike objects have an attribute native_event that returns
+    a cl.Event that indicates the end of the event.
+    """
+    native_event: cl.Event
 
 
 class UnableToCollectTimingData(UserWarning):
@@ -216,12 +220,12 @@ class UnableToCollectTimingData(UserWarning):
 
 class SumpyTimingFuture:
 
-    def __init__(self, queue, events: List[Union[cl.Event, EventLike]]):
+    def __init__(self, queue, events: list[cl.Event | CLEventLike]):
         self.queue = queue
         self.events = events
 
     @property
-    def native_events(self) -> List[cl.Event]:
+    def native_events(self) -> list[cl.Event]:
         return [evt if isinstance(evt, cl.Event) else evt.native_event
                 for evt in self.events]
 
diff --git a/sumpy/kernel.py b/sumpy/kernel.py
index 1cf2c29e..f4a3cbe4 100644
--- a/sumpy/kernel.py
+++ b/sumpy/kernel.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2012 Andreas Kloeckner"
 
 __license__ = """
@@ -20,10 +23,12 @@
 THE SOFTWARE.
 """
 
+from abc import abstractmethod
 from collections import defaultdict
-from typing import ClassVar, Tuple
+from typing import TYPE_CHECKING, ClassVar
 
 import numpy as np
+import sympy as sp
 
 import loopy as lp
 from pymbolic import var
@@ -35,6 +40,10 @@
 from sumpy.symbolic import SpatialConstant, pymbolic_real_norm_2
 
 
+if TYPE_CHECKING:
+    from sumpy.expansion.diff_op import LinearPDESystemOperator
+
+
 __doc__ = """
 Kernel interface
 ----------------
@@ -140,9 +149,16 @@ class Kernel:
     .. automethod:: get_source_args
     """
 
-    init_arg_names: ClassVar[Tuple[str, ...]]
+    if TYPE_CHECKING:
+        @property
+        def is_complex_valued(self) -> bool:
+            ...
 
-    def __init__(self, dim):
+    dim: int
+    init_arg_names: ClassVar[tuple[str, ...]]
+    _hash_value: int
+
+    def __init__(self, dim: int) -> None:
         self.dim = dim
 
     # {{{ hashing/pickling/equality
@@ -159,12 +175,12 @@ def __eq__(self, other):
     def __ne__(self, other):
         return not self.__eq__(other)
 
-    def __hash__(self):
+    def __hash__(self) -> int:
         try:
-            return self.hash_value
+            return self._hash_value
         except AttributeError:
-            self.hash_value = hash((type(self), *self.__getinitargs__()))
-            return self.hash_value
+            self._hash_value = hash((type(self), *self.__getinitargs__()))
+            return self._hash_value
 
     def update_persistent_hash(self, key_hash, key_builder):
         key_hash.update(type(self).__name__.encode("utf8"))
@@ -183,13 +199,13 @@ def __setstate__(self, state):
 
     # }}}
 
-    def get_base_kernel(self):
+    def get_base_kernel(self) -> Kernel:
         """Return the kernel being wrapped by this one, or else
         *self*.
         """
         return self
 
-    def replace_base_kernel(self, new_base_kernel):
+    def replace_base_kernel(self, new_base_kernel: Kernel) -> Kernel:
         """Return the base kernel being wrapped by this one, or else
         *new_base_kernel*.
         """
@@ -208,9 +224,10 @@ def get_code_transformer(self):
         """
         return lambda expr: expr
 
-    def get_expression(self, dist_vec):
+    @abstractmethod
+    def get_expression(self, dist_vec: np.ndarray) -> sp.Expr:
         r"""Return a :mod:`sympy` expression for the kernel."""
-        raise NotImplementedError
+        ...
 
     def _diff(self, expr, vec, mi):
         """Take the derivative of an expression
@@ -268,6 +285,7 @@ def get_derivative_coeff_dict_at_source(self, expr_dict):
         """
         return expr_dict
 
+    @abstractmethod
     def get_global_scaling_const(self):
         r"""Return a global scaling constant of the kernel.
         Typically, this ensures that the kernel is scaled so that
@@ -276,21 +294,25 @@ def get_global_scaling_const(self):
         Not to be confused with *rscale*, which keeps expansion
         coefficients benignly scaled.
         """
-        raise NotImplementedError
+        ...
 
-    def get_args(self):
+    def get_args(self) -> list[KernelArgument]:
         """Return list of :class:`KernelArgument` instances describing
         extra arguments used by the kernel.
         """
         return []
 
-    def get_source_args(self):
+    def get_source_args(self) -> list[KernelArgument]:
         """Return list of :class:`KernelArgument` instances describing
         extra arguments used by kernel in picking up contributions
         from point sources.
         """
         return []
 
+    @abstractmethod
+    def get_pde_as_diff_op(self) -> LinearPDESystemOperator:
+        ...
+
     # TODO: Allow kernels that are not translation invariant
     is_translation_invariant = True
 
@@ -324,8 +346,8 @@ class ExpressionKernel(Kernel):
     .. automethod:: get_expression
     """
 
-    init_arg_names = ("dim", "expression", "global_scaling_const",
-            "is_complex_valued")
+    init_arg_names: ClassVar[tuple[str, ...]] = (
+            "dim", "expression", "global_scaling_const", "is_complex_valued")
 
     def __init__(self, dim, expression, global_scaling_const,
             is_complex_valued):
diff --git a/sumpy/p2e.py b/sumpy/p2e.py
index 4a481609..28a64f42 100644
--- a/sumpy/p2e.py
+++ b/sumpy/p2e.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2013 Andreas Kloeckner"
 
 __license__ = """
diff --git a/sumpy/p2p.py b/sumpy/p2p.py
index 6298da34..0db70d5c 100644
--- a/sumpy/p2p.py
+++ b/sumpy/p2p.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = """
 Copyright (C) 2012 Andreas Kloeckner
 Copyright (C) 2018 Alexandru Fikl
diff --git a/sumpy/point_calculus.py b/sumpy/point_calculus.py
index 238f0341..3694f89e 100644
--- a/sumpy/point_calculus.py
+++ b/sumpy/point_calculus.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2017 Andreas Kloeckner"
 
 __license__ = """
diff --git a/sumpy/qbx.py b/sumpy/qbx.py
index f587905f..ea3cded7 100644
--- a/sumpy/qbx.py
+++ b/sumpy/qbx.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = """
 Copyright (C) 2012 Andreas Kloeckner
 Copyright (C) 2018 Alexandru Fikl
diff --git a/sumpy/symbolic.py b/sumpy/symbolic.py
index 4e1b5276..43755e53 100644
--- a/sumpy/symbolic.py
+++ b/sumpy/symbolic.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2012 Andreas Kloeckner"
 
 __license__ = """
@@ -36,6 +39,7 @@
 
 import logging
 import math
+from typing import TYPE_CHECKING
 
 import numpy as np
 
@@ -135,7 +139,7 @@ def _coeff_isneg(a):
     return a.is_Number and a.is_negative
 
 
-if USE_SYMENGINE:
+if TYPE_CHECKING or USE_SYMENGINE:
     def UnevaluatedExpr(x):  # noqa: N802
         return x
 else:
@@ -379,11 +383,12 @@ class Hankel1(_BesselOrHankel):
 _SympyBesselJ = BesselJ
 _SympyHankel1 = Hankel1
 
-if USE_SYMENGINE:
-    def BesselJ(*args):   # noqa: N802  # pylint: disable=function-redefined
-        return sym.sympify(_SympyBesselJ(*args))
+if not TYPE_CHECKING:
+    if USE_SYMENGINE:
+        def BesselJ(*args):   # noqa: N802  # pylint: disable=function-redefined
+            return sym.sympify(_SympyBesselJ(*args))
 
-    def Hankel1(*args):   # noqa: N802  # pylint: disable=function-redefined
-        return sym.sympify(_SympyHankel1(*args))
+        def Hankel1(*args):   # noqa: N802  # pylint: disable=function-redefined
+            return sym.sympify(_SympyHankel1(*args))
 
 # vim: fdm=marker
diff --git a/sumpy/tools.py b/sumpy/tools.py
index c9d9d0a4..9ceffe80 100644
--- a/sumpy/tools.py
+++ b/sumpy/tools.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = """
 Copyright (C) 2012 Andreas Kloeckner
 Copyright (C) 2018 Alexandru Fikl
@@ -29,11 +32,12 @@
 import warnings
 from abc import ABC, abstractmethod
 from dataclasses import dataclass
-from typing import TYPE_CHECKING, Any, List, Optional, Sequence, Tuple
+from typing import TYPE_CHECKING, Any, Hashable, Sequence
 
 import numpy as np
 
 import loopy as lp
+import pyopencl as cl
 from pymbolic.mapper import WalkMapper
 from pytools import memoize_method
 from pytools.tag import Tag, tag_dataclass
@@ -111,7 +115,7 @@
 
 # {{{ multi_index helpers
 
-def add_mi(mi1: Sequence[int], mi2: Sequence[int]) -> Tuple[int, ...]:
+def add_mi(mi1: Sequence[int], mi2: Sequence[int]) -> tuple[int, ...]:
     return tuple([mi1i + mi2i for mi1i, mi2i in zip(mi1, mi2)])
 
 
@@ -125,13 +129,13 @@ def mi_factorial(mi: Sequence[int]) -> int:
 
 def mi_increment_axis(
         mi: Sequence[int], axis: int, increment: int
-        ) -> Tuple[int, ...]:
+        ) -> tuple[int, ...]:
     new_mi = list(mi)
     new_mi[axis] += increment
     return tuple(new_mi)
 
 
-def mi_set_axis(mi: Sequence[int], axis: int, value: int) -> Tuple[int, ...]:
+def mi_set_axis(mi: Sequence[int], axis: int, value: int) -> tuple[int, ...]:
     new_mi = list(mi)
     new_mi[axis] = value
     return tuple(new_mi)
@@ -284,12 +288,12 @@ class KernelComputation(ABC):
     """
 
     def __init__(self, ctx: Any,
-            target_kernels: List["Kernel"],
-            source_kernels: List["Kernel"],
-            strength_usage: Optional[List[int]] = None,
-            value_dtypes: Optional[List["numpy.dtype[Any]"]] = None,
-            name: Optional[str] = None,
-            device: Optional[Any] = None) -> None:
+            target_kernels: list[Kernel],
+            source_kernels: list[Kernel],
+            strength_usage: list[int] | None = None,
+            value_dtypes: list[numpy.dtype[Any]] | None = None,
+            name: str | None = None,
+            device: Any | None = None) -> None:
         """
         :arg target_kernels: list of :class:`~sumpy.kernel.Kernel` instances,
             with :class:`sumpy.kernel.DirectionalTargetDerivative` as
@@ -309,9 +313,9 @@ def __init__(self, ctx: Any,
             value_dtypes = []
             for knl in target_kernels:
                 if knl.is_complex_valued:
-                    value_dtypes.append(np.complex128)
+                    value_dtypes.append(np.dtype(np.complex128))
                 else:
-                    value_dtypes.append(np.float64)
+                    value_dtypes.append(np.dtype(np.float64))
 
         if not isinstance(value_dtypes, (list, tuple)):
             value_dtypes = [np.dtype(value_dtypes)] * len(target_kernels)
@@ -440,15 +444,21 @@ def __eq__(self, other):
 # }}}
 
 
-class KernelCacheMixin:
-    def get_cached_optimized_kernel(self, **kwargs):
-        from warnings import warn
-        warn("get_cached_optimized_kernel is deprecated. "
-             "Use get_cached_kernel_executor instead. "
-             "This will stop working in October 2023.",
-             DeprecationWarning, stacklevel=2)
+class KernelCacheMixin(ABC):
+    context: cl.Context
+    name: str
+
+    @abstractmethod
+    def get_cache_key(self) -> tuple[Hashable, ...]:
+        ...
 
-        return self.get_cached_kernel_executor(**kwargs)
+    @abstractmethod
+    def get_kernel(self, **kwargs: Any) -> lp.TranslationUnit:
+        ...
+
+    @abstractmethod
+    def get_optimized_kernel(self, **kwargs: Any) -> lp.TranslationUnit:
+        ...
 
     @memoize_method
     def get_cached_kernel_executor(self, **kwargs) -> lp.ExecutorBase:
@@ -890,7 +900,7 @@ class FFTBackend(enum.Enum):
     loopy = 2
 
 
-def _get_fft_backend(queue: "pyopencl.CommandQueue") -> FFTBackend:
+def _get_fft_backend(queue: pyopencl.CommandQueue) -> FFTBackend:
     import os
 
     env_val = os.environ.get("SUMPY_FFT_BACKEND")
@@ -931,9 +941,9 @@ def _get_fft_backend(queue: "pyopencl.CommandQueue") -> FFTBackend:
 
 
 def get_opencl_fft_app(
-        queue: "pyopencl.CommandQueue",
-        shape: Tuple[int, ...],
-        dtype: "numpy.dtype[Any]",
+        queue: pyopencl.CommandQueue,
+        shape: tuple[int, ...],
+        dtype: numpy.dtype[Any],
         inverse: bool) -> Any:
     """Setup an object for out-of-place FFT on with given shape and dtype
     on given queue.
@@ -954,12 +964,12 @@ def get_opencl_fft_app(
 
 
 def run_opencl_fft(
-        fft_app: Tuple[Any, FFTBackend],
-        queue: "pyopencl.CommandQueue",
+        fft_app: tuple[Any, FFTBackend],
+        queue: pyopencl.CommandQueue,
         input_vec: Any,
         inverse: bool = False,
-        wait_for: Optional[List["pyopencl.Event"]] = None
-    ) -> Tuple["pyopencl.Event", Any]:
+        wait_for: list[pyopencl.Event] | None = None
+    ) -> tuple[pyopencl.Event, Any]:
     """Runs an FFT on input_vec and returns a :class:`MarkerBasedProfilingEvent`
     that indicate the end and start of the operations carried out and the output
     vector.
diff --git a/sumpy/toys.py b/sumpy/toys.py
index 7f0d170c..4f569ef3 100644
--- a/sumpy/toys.py
+++ b/sumpy/toys.py
@@ -1,5 +1,7 @@
 from __future__ import annotations
 
+from sumpy.expansion.m2l import M2LTranslationClassFactoryBase
+
 
 __copyright__ = """
 Copyright (C) 2017 Andreas Kloeckner
@@ -28,7 +30,7 @@
 
 from functools import partial
 from numbers import Number
-from typing import TYPE_CHECKING, Sequence
+from typing import TYPE_CHECKING, Sequence, Union
 
 from pytools import memoize_method
 
@@ -128,7 +130,8 @@ def __init__(self, cl_context: pyopencl.Context, kernel: Kernel,
                 NonFFTM2LTranslationClassFactory,
             )
             if m2l_use_fft:
-                m2l_translation_class_factory = FFTM2LTranslationClassFactory()
+                m2l_translation_class_factory: M2LTranslationClassFactoryBase = \
+                    FFTM2LTranslationClassFactory()
             else:
                 m2l_translation_class_factory = NonFFTM2LTranslationClassFactory()
             local_expn_class = \
@@ -478,6 +481,9 @@ def _m2l(psource, to_center, to_rscale, to_order, e2e, expn_class, expn_kwargs,
 
 # {{{ potential source classes
 
+Number_ish = Union[int, float, complex, np.number]
+
+
 class PotentialSource:
     """A base class for all classes representing potentials that can be
     evaluated anywhere in space.
@@ -508,7 +514,7 @@ def eval(self, targets: np.ndarray) -> np.ndarray:
     def __neg__(self) -> PotentialSource:
         return -1*self
 
-    def __add__(self, other: Number | np.number | PotentialSource
+    def __add__(self, other: Number_ish | PotentialSource
                 ) -> PotentialSource:
         if isinstance(other, (Number, np.number)):
             other = ConstantPotential(self.toy_ctx, other)
@@ -520,16 +526,18 @@ def __add__(self, other: Number | np.number | PotentialSource
     __radd__ = __add__
 
     def __sub__(self,
-                other: Number | np.number | PotentialSource) -> PotentialSource:
+                other: Number_ish | PotentialSource) -> PotentialSource:
         return self.__add__(-other)
 
-    def __rsub__(self,
-                 other: Number | np.number | PotentialSource
+    # FIXME: mypy says " Forward operator "__sub__" is not callable"
+    # I don't know what it means. -AK, 2024-07-18
+    def __rsub__(self,  # type:ignore[misc]
+                 other: Number_ish | PotentialSource
                  ) -> PotentialSource:
         return (-self).__add__(other)
 
     def __mul__(self,
-                other: Number | np.number | PotentialSource) -> PotentialSource:
+                other: Number_ish | PotentialSource) -> PotentialSource:
         if isinstance(other, (Number, np.number)):
             other = ConstantPotential(self.toy_ctx, other)
         elif not isinstance(other, PotentialSource):
@@ -705,7 +713,7 @@ def __init__(self, psources: Sequence[PotentialSource]) -> None:
     def center(self) -> np.ndarray:
         for psource in self.psources:
             try:
-                return psource.center
+                return psource.center  # type: ignore[attr-defined]
             except AttributeError:
                 pass
 
@@ -718,7 +726,7 @@ class Sum(PotentialExpressionNode):
     """
 
     def eval(self, targets: np.ndarray) -> np.ndarray:
-        result = 0
+        result = np.zeros(targets.shape[1])
         for psource in self.psources:
             result = result + psource.eval(targets)
 
@@ -731,7 +739,7 @@ class Product(PotentialExpressionNode):
     """
 
     def eval(self, targets: np.ndarray) -> np.ndarray:
-        result = 1
+        result = np.zeros(targets.shape[1])
         for psource in self.psources:
             result = result * psource.eval(targets)
 
@@ -764,7 +772,7 @@ def multipole_expand(psource: PotentialSource,
 
 def local_expand(
         psource: PotentialSource,
-        center: np.ndarray, order: int | None = None, rscale: Number = 1,
+        center: np.ndarray, order: int | None = None, rscale: Number_ish = 1,
         **expn_kwargs) -> LocalExpansion:
     if isinstance(psource, PointSources):
         if order is None:
@@ -823,8 +831,10 @@ def combine_inner_outer(
         radius: float | None,
         center: np.ndarray | None = None) -> PotentialSource:
     if center is None:
+        assert isinstance(psource_inner, ExpansionPotentialSource)
         center = psource_inner.center
     if radius is None:
+        assert isinstance(psource_inner, ExpansionPotentialSource)
         radius = psource_inner.radius
 
     ball_one = OneOnBallPotential(psource_inner.toy_ctx, center, radius)
@@ -837,6 +847,7 @@ def combine_halfspace(psource_pos: PotentialSource,
                       psource_neg: PotentialSource, axis: int,
                       center: np.ndarray | None = None) -> PotentialSource:
     if center is None:
+        assert isinstance(psource_pos, ExpansionPotentialSource)
         center = psource_pos.center
 
     halfspace_one = HalfspaceOnePotential(psource_pos.toy_ctx, center, axis)
@@ -849,12 +860,14 @@ def combine_halfspace_and_outer(
         psource_pos: PotentialSource,
         psource_neg: PotentialSource,
         psource_outer: PotentialSource,
-        axis: int, radius: Number | None = None,
+        axis: int, radius: float | None = None,
         center: np.ndarray | None = None) -> PotentialSource:
 
     if center is None:
+        assert isinstance(psource_pos, ExpansionPotentialSource)
         center = psource_pos.center
     if radius is None:
+        assert isinstance(psource_pos, ExpansionPotentialSource)
         center = psource_pos.radius
 
     return combine_inner_outer(
@@ -866,6 +879,7 @@ def l_inf(psource: PotentialSource, radius: float,
           center: np.ndarray | None = None, npoints: int = 100,
           debug: bool = False) -> np.number:
     if center is None:
+        assert isinstance(psource, ExpansionPotentialSource)
         center = psource.center
 
     restr = psource * OneOnBallPotential(psource.toy_ctx, center, radius)
diff --git a/sumpy/version.py b/sumpy/version.py
index a812f4d6..eb6989b9 100644
--- a/sumpy/version.py
+++ b/sumpy/version.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2014 Andreas Kloeckner"
 
 __license__ = """
diff --git a/sumpy/visualization.py b/sumpy/visualization.py
index 02803a8e..d5f2f6c0 100644
--- a/sumpy/visualization.py
+++ b/sumpy/visualization.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2012 Andreas Kloeckner"
 
 __license__ = """
diff --git a/test/test_codegen.py b/test/test_codegen.py
index ff73bc73..4eb45b09 100644
--- a/test/test_codegen.py
+++ b/test/test_codegen.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2017 Matt Wala"
 
 __license__ = """
diff --git a/test/test_cse.py b/test/test_cse.py
index e28bf8e4..15330e44 100644
--- a/test/test_cse.py
+++ b/test/test_cse.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = """
 Copyright (C) 2017 Matt Wala
 Copyright (C) 2006-2016 SymPy Development Team
diff --git a/test/test_distributed.py b/test/test_distributed.py
index 60e94fff..9181f4cf 100644
--- a/test/test_distributed.py
+++ b/test/test_distributed.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2022 Hao Gao"
 
 __license__ = """
diff --git a/test/test_fmm.py b/test/test_fmm.py
index 433fca8a..437ca9ea 100644
--- a/test/test_fmm.py
+++ b/test/test_fmm.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2013 Andreas Kloeckner"
 
 __license__ = """
diff --git a/test/test_kernels.py b/test/test_kernels.py
index 6edb07f3..812f6ac7 100644
--- a/test/test_kernels.py
+++ b/test/test_kernels.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2012 Andreas Kloeckner"
 
 __license__ = """
diff --git a/test/test_matrixgen.py b/test/test_matrixgen.py
index 516a9581..1dc0dd85 100644
--- a/test/test_matrixgen.py
+++ b/test/test_matrixgen.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2018 Alexandru Fikl"
 
 __license__ = """
diff --git a/test/test_misc.py b/test/test_misc.py
index 302f5b85..06b45f64 100644
--- a/test/test_misc.py
+++ b/test/test_misc.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2017 Andreas Kloeckner"
 
 __license__ = """
diff --git a/test/test_qbx.py b/test/test_qbx.py
index 88c226bc..d48c1bdb 100644
--- a/test/test_qbx.py
+++ b/test/test_qbx.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2017 Matt Wala"
 
 __license__ = """
diff --git a/test/test_tools.py b/test/test_tools.py
index 309121c5..ab4d5bf4 100644
--- a/test/test_tools.py
+++ b/test/test_tools.py
@@ -1,3 +1,6 @@
+from __future__ import annotations
+
+
 __copyright__ = "Copyright (C) 2020 Isuru Fernando"
 
 __license__ = """