diff --git a/docs/source/run_autodoc.py b/docs/source/run_autodoc.py
index af51cd8..6f08c64 100755
--- a/docs/source/run_autodoc.py
+++ b/docs/source/run_autodoc.py
@@ -27,7 +27,7 @@ def replace_in_file(infile, replacements, outfile=None):
     if outfile is None:
         outfile = infile
 
-    with open(infile, "r") as fp:
+    with open(infile) as fp:
         content = fp.read()
 
     for key, value in replacements.items():
diff --git a/docs/source/sphinx_simplify_typehints.py b/docs/source/sphinx_simplify_typehints.py
index b9080e3..e59bc02 100644
--- a/docs/source/sphinx_simplify_typehints.py
+++ b/docs/source/sphinx_simplify_typehints.py
@@ -80,7 +80,7 @@
 # replacement rules based on regular expressions
 REPLACEMENTS_REGEX = {
     # remove full package path and only leave the module/class identifier
-    "pde\.(\w+\.)*": "",
+    r"pde\.(\w+\.)*": "",
 }
 
 
diff --git a/droplets/droplet_tracks.py b/droplets/droplet_tracks.py
index 165c7c9..4f7f20a 100644
--- a/droplets/droplet_tracks.py
+++ b/droplets/droplet_tracks.py
@@ -15,7 +15,7 @@
 
 import json
 import logging
-from typing import Callable, List, Literal, Optional
+from typing import Callable, Literal
 
 import numpy as np
 from numpy.lib import recfunctions as rfn
@@ -159,7 +159,7 @@ def last(self) -> SphericalDroplet:
         return self.droplets[-1]  # type: ignore
 
     @property
-    def dim(self) -> Optional[int]:
+    def dim(self) -> int | None:
         """return the space dimension of the droplets"""
         try:
             return self.last.dim
@@ -167,7 +167,7 @@ def dim(self) -> Optional[int]:
             return None
 
     @property
-    def data(self) -> Optional[np.ndarray]:
+    def data(self) -> np.ndarray | None:
         """:class:`~numpy.ndarray`: an array containing the data of the full track"""
         if len(self) == 0:
             return None
@@ -187,7 +187,7 @@ def items(self):
         """iterate over all times and droplets, returning them in pairs"""
         return zip(self.times, self.droplets)
 
-    def append(self, droplet: SphericalDroplet, time: Optional[float] = None) -> None:
+    def append(self, droplet: SphericalDroplet, time: float | None = None) -> None:
         """append a new droplet with a time code
 
         Args:
@@ -260,7 +260,7 @@ def get_volumes(self, smoothing: float = 0) -> np.ndarray:
         """
         return self.get_trajectory(smoothing, attribute="volume")
 
-    def time_overlaps(self, other: "DropletTrack") -> bool:
+    def time_overlaps(self, other: DropletTrack) -> bool:
         """determine whether two DropletTrack instances overlaps in time
 
         Args:
@@ -334,7 +334,7 @@ def _write_hdf_dataset(self, hdf_path, key: str = "droplet_track"):
 
         return dataset
 
-    def to_file(self, path: str, info: Optional[InfoDict] = None) -> None:
+    def to_file(self, path: str, info: InfoDict | None = None) -> None:
         """store data in hdf5 file
 
         The data can be read using the classmethod :meth:`DropletTrack.from_file`.
@@ -360,7 +360,7 @@ def plot(
         self,
         attribute: str = "radius",
         smoothing: float = 0,
-        t_max: Optional[float] = None,
+        t_max: float | None = None,
         ax=None,
         **kwargs,
     ) -> PlotReference:
@@ -409,7 +409,7 @@ def plot(
 
     @plot_on_axes()
     def plot_positions(
-        self, grid: Optional[GridBase] = None, arrow: bool = True, ax=None, **kwargs
+        self, grid: GridBase | None = None, arrow: bool = True, ax=None, **kwargs
     ) -> PlotReference:
         """plot the droplet track
 
@@ -495,7 +495,7 @@ def from_emulsion_time_course(
         time_course: EmulsionTimeCourse,
         *,
         method: Literal["distance", "overlap"] = "overlap",
-        grid: Optional[GridBase] = None,
+        grid: GridBase | None = None,
         progress: bool = False,
         **kwargs,
     ) -> DropletTrackList:
@@ -531,13 +531,13 @@ def from_emulsion_time_course(
             # track droplets by their physical overlap
 
             def match_tracks(
-                emulsion: Emulsion, tracks_alive: List[DropletTrack], time: float
+                emulsion: Emulsion, tracks_alive: list[DropletTrack], time: float
             ) -> None:
                 """helper function adding emulsions to the tracks"""
                 found_multiple_overlap = False
                 for droplet in emulsion:
                     # determine which old tracks could be extended
-                    overlaps: List[DropletTrack] = []
+                    overlaps: list[DropletTrack] = []
                     for track in tracks_alive:
                         if track.last.overlaps(droplet, grid=grid):
                             overlaps.append(track)
@@ -558,7 +558,7 @@ def match_tracks(
             max_dist = kwargs.pop("max_dist", np.inf)
 
             def match_tracks(
-                emulsion: Emulsion, tracks_alive: List[DropletTrack], time: float
+                emulsion: Emulsion, tracks_alive: list[DropletTrack], time: float
             ) -> None:
                 """helper function adding emulsions to the tracks"""
                 added = set()
@@ -619,7 +619,7 @@ def from_storage(
         method: Literal["distance", "overlap"] = "overlap",
         refine: bool = False,
         num_processes: int | Literal["auto"] = 1,
-        progress: Optional[bool] = None,
+        progress: bool | None = None,
     ) -> DropletTrackList:
         r"""obtain droplet tracks from stored scalar field data
 
@@ -675,7 +675,7 @@ def from_file(cls, path: str) -> DropletTrackList:
                 obj.append(DropletTrack._from_hdf_dataset(dataset))
         return obj
 
-    def to_file(self, path: str, info: Optional[InfoDict] = None) -> None:
+    def to_file(self, path: str, info: InfoDict | None = None) -> None:
         """store data in hdf5 file
 
         The data can be read using the classmethod :meth:`DropletTrackList.from_file`.
diff --git a/droplets/droplets.py b/droplets/droplets.py
index 604ad3a..0b836f9 100644
--- a/droplets/droplets.py
+++ b/droplets/droplets.py
@@ -32,7 +32,7 @@
 import warnings
 from abc import ABCMeta, abstractmethod
 from pathlib import Path
-from typing import Any, Callable, Dict, List, Optional, Tuple, Type, TypeVar, Union
+from typing import Any, Callable, List, Tuple, Type, TypeVar, Union
 
 import numpy as np
 from numba.extending import register_jitable
@@ -103,7 +103,7 @@ class DropletBase:
     `dtype` of the array determines what information the droplet class stores.
     """
 
-    _subclasses: Dict[str, DropletBase] = {}  # collect all inheriting classes
+    _subclasses: dict[str, DropletBase] = {}  # collect all inheriting classes
 
     __slots__ = ["data"]
 
@@ -238,7 +238,7 @@ def merge(self: TDroplet, other: TDroplet, *, inplace: bool = False) -> TDroplet
             return self.__class__.from_data(result)  # type: ignore
 
     @property
-    def data_bounds(self) -> Tuple[np.ndarray, np.ndarray]:
+    def data_bounds(self) -> tuple[np.ndarray, np.ndarray]:
         """tuple: lower and upper bounds on the parameters"""
         num = len(self._data_array)
         return np.full(num, -np.inf), np.full(num, np.inf)
@@ -290,7 +290,7 @@ def dim(self) -> int:
         return get_dtype_field_size(self.data.dtype, "position")
 
     @property
-    def data_bounds(self) -> Tuple[np.ndarray, np.ndarray]:
+    def data_bounds(self) -> tuple[np.ndarray, np.ndarray]:
         """tuple: lower and upper bounds on the parameters"""
         l, h = super().data_bounds
         l[self.dim] = 0  # radius must be non-negative
@@ -356,9 +356,7 @@ def bbox(self) -> Cuboid:
             self.position - self.radius, self.position + self.radius
         )
 
-    def overlaps(
-        self, other: SphericalDroplet, grid: Optional[GridBase] = None
-    ) -> bool:
+    def overlaps(self, other: SphericalDroplet, grid: GridBase | None = None) -> bool:
         """determine whether another droplet overlaps with this one
 
         Note that this function so far only compares the distances of the droplets to
@@ -475,7 +473,7 @@ def get_phase_field(
         *,
         vmin: float = 0,
         vmax: float = 1,
-        label: Optional[str] = None,
+        label: str | None = None,
     ) -> ScalarField:
         """Creates an image of the droplet on the `grid`
 
@@ -497,7 +495,7 @@ def get_phase_field(
         data = vmin + (vmax - vmin) * data  # scale data
         return ScalarField(grid, data=data, label=label)
 
-    def get_triangulation(self, resolution: float = 1) -> Dict[str, Any]:
+    def get_triangulation(self, resolution: float = 1) -> dict[str, Any]:
         """obtain a triangulated shape of the droplet surface
 
         Args:
@@ -567,7 +565,7 @@ def _get_mpl_patch(self, dim=None, **kwargs):
         return mpl.patches.Circle(position, self.radius, **kwargs)
 
     @plot_on_axes()
-    def plot(self, ax, value: Optional[Callable] = None, **kwargs) -> PlotReference:
+    def plot(self, ax, value: Callable | None = None, **kwargs) -> PlotReference:
         """Plot the droplet
 
         Args:
@@ -601,7 +599,7 @@ def __init__(
         self,
         position: np.ndarray,
         radius: float,
-        interface_width: Optional[float] = None,
+        interface_width: float | None = None,
     ):
         """
         Args:
@@ -617,7 +615,7 @@ def __init__(
         self.interface_width = interface_width
 
     @property
-    def data_bounds(self) -> Tuple[np.ndarray, np.ndarray]:
+    def data_bounds(self) -> tuple[np.ndarray, np.ndarray]:
         """tuple: lower and upper bounds on the parameters"""
         l, h = super().data_bounds
         l[self.dim + 1] = 0  # interface width must be non-negative
@@ -651,7 +649,7 @@ def merge_data(drop1: np.ndarray, drop2: np.ndarray, out: np.ndarray) -> None:
         return merge_data  # type: ignore
 
     @property
-    def interface_width(self) -> Optional[float]:
+    def interface_width(self) -> float | None:
         """float: the width of the interface of this droplet"""
         if np.isnan(self.data["interface_width"]):
             return None
@@ -659,7 +657,7 @@ def interface_width(self) -> Optional[float]:
             return float(self.data["interface_width"])
 
     @interface_width.setter
-    def interface_width(self, value: Optional[float]) -> None:
+    def interface_width(self, value: float | None) -> None:
         if value is None:
             self.data["interface_width"] = math.nan
         elif value < 0:
@@ -718,8 +716,8 @@ def __init__(
         self,
         position: np.ndarray,
         radius: float,
-        interface_width: Optional[float] = None,
-        amplitudes: Optional[np.ndarray] = None,
+        interface_width: float | None = None,
+        amplitudes: np.ndarray | None = None,
     ):
         """
         Args:
@@ -768,7 +766,7 @@ def get_dtype(cls, **kwargs) -> DTypeList:
         return dtype + [("amplitudes", float, (modes,))]
 
     @property
-    def data_bounds(self) -> Tuple[np.ndarray, np.ndarray]:
+    def data_bounds(self) -> tuple[np.ndarray, np.ndarray]:
         """tuple: lower and upper bounds on the parameters"""
         l, h = super().data_bounds
         n = self.dim + 2
@@ -789,7 +787,7 @@ def amplitudes(self) -> np.ndarray:
         return np.atleast_1d(self.data["amplitudes"])  # type: ignore
 
     @amplitudes.setter
-    def amplitudes(self, value: Optional[np.ndarray] = None) -> None:
+    def amplitudes(self, value: np.ndarray | None = None) -> None:
         if value is None:
             assert self.modes == 0
             self.data["amplitudes"] = np.broadcast_to(0.0, (0,))
@@ -888,8 +886,8 @@ def __init__(
         self,
         position: np.ndarray,
         radius: float,
-        interface_width: Optional[float] = None,
-        amplitudes: Optional[np.ndarray] = None,
+        interface_width: float | None = None,
+        amplitudes: np.ndarray | None = None,
     ):
         r"""
         Args:
@@ -1066,8 +1064,8 @@ def __init__(
         self,
         position: np.ndarray,
         radius: float,
-        interface_width: Optional[float] = None,
-        amplitudes: Optional[np.ndarray] = None,
+        interface_width: float | None = None,
+        amplitudes: np.ndarray | None = None,
     ):
         r"""
         Args:
@@ -1097,7 +1095,7 @@ def __init__(
 
     @preserve_scalars
     def interface_distance(  # type: ignore
-        self, θ: np.ndarray, φ: Optional[np.ndarray] = None
+        self, θ: np.ndarray, φ: np.ndarray | None = None
     ) -> np.ndarray:
         r"""calculates the distance of the droplet interface to the origin
 
@@ -1122,7 +1120,7 @@ def interface_distance(  # type: ignore
 
     @preserve_scalars
     def interface_position(
-        self, θ: np.ndarray, φ: Optional[np.ndarray] = None
+        self, θ: np.ndarray, φ: np.ndarray | None = None
     ) -> np.ndarray:
         r"""calculates the position of the interface of the droplet
 
@@ -1146,7 +1144,7 @@ def interface_position(
 
     @preserve_scalars
     def interface_curvature(  # type: ignore
-        self, θ: np.ndarray, φ: Optional[np.ndarray] = None
+        self, θ: np.ndarray, φ: np.ndarray | None = None
     ) -> np.ndarray:
         r"""calculates the mean curvature of the interface of the droplet
 
@@ -1296,7 +1294,7 @@ class _TriangulatedSpheres:
     def __init__(self) -> None:
         self.path = Path(__file__).resolve().parent / "resources" / "spheres_3d.hdf5"
         self.num_list = np.zeros((0,))
-        self.data: Optional[Dict[int, Dict[str, Any]]] = None
+        self.data: dict[int, dict[str, Any]] | None = None
 
     def _load(self):
         """load the stored resource"""
@@ -1316,7 +1314,7 @@ def _load(self):
                 }
                 self.data[num] = tri
 
-    def get_triangulation(self, num_est: int = 1) -> Dict[str, Any]:
+    def get_triangulation(self, num_est: int = 1) -> dict[str, Any]:
         """get a triangulation of a sphere
 
         Args:
diff --git a/droplets/emulsions.py b/droplets/emulsions.py
index a78e92c..a37a570 100644
--- a/droplets/emulsions.py
+++ b/droplets/emulsions.py
@@ -21,16 +21,10 @@
     TYPE_CHECKING,
     Any,
     Callable,
-    Dict,
     Iterable,
     Iterator,
-    List,
     Literal,
-    Optional,
     Sequence,
-    Tuple,
-    Type,
-    Union,
     overload,
 )
 
@@ -61,12 +55,12 @@ class Emulsion(list):
 
     def __init__(
         self,
-        droplets: Optional[Iterable[SphericalDroplet]] = None,
+        droplets: Iterable[SphericalDroplet] | None = None,
         *,
         copy: bool = True,
         dtype: np.typing.DTypeLike | np.ndarray | SphericalDroplet = None,
         force_consistency: bool = False,
-        grid: Optional[GridBase] = None,
+        grid: GridBase | None = None,
     ):
         """
         Args:
@@ -123,12 +117,12 @@ def empty(cls, droplet: SphericalDroplet) -> Emulsion:
     def from_random(
         cls,
         num: int,
-        grid_or_bounds: GridBase | Sequence[Tuple[float, float]],
-        radius: float | Tuple[float, float],
+        grid_or_bounds: GridBase | Sequence[tuple[float, float]],
+        radius: float | tuple[float, float],
         *,
         remove_overlapping: bool = True,
-        droplet_class: Type[SphericalDroplet] = SphericalDroplet,
-        rng: Optional[np.random.Generator] = None,
+        droplet_class: type[SphericalDroplet] = SphericalDroplet,
+        rng: np.random.Generator | None = None,
     ) -> Emulsion:
         """
         Create an emulsion with random droplets
@@ -188,7 +182,7 @@ def get_position():
         return emulsion
 
     @property
-    def dim(self) -> Optional[int]:
+    def dim(self) -> int | None:
         """int: dimensionality of space in which droplets are defined"""
         if self.dtype:
             return self.dtype["position"].shape[0]  # type: ignore
@@ -227,7 +221,7 @@ def copy(self, min_radius: float = -1) -> Emulsion:
                 exactly min_radius are removed, so `min_radius == 0` can be used to
                 filter vanished droplets.
         """
-        droplets: List[SphericalDroplet] = [
+        droplets: list[SphericalDroplet] = [
             droplet.copy() for droplet in self if droplet.radius > min_radius
         ]
         return self.__class__(droplets, copy=False)
@@ -303,7 +297,7 @@ def data(self) -> np.ndarray:
 
         else:
             # emulsion contains at least one droplet
-            classes = set(d.__class__ for d in self)
+            classes = {d.__class__ for d in self}
             if len(classes) > 1:
                 raise TypeError(
                     "Emulsion data cannot be stored contiguously if it contains a "
@@ -343,7 +337,7 @@ def _from_hdf_dataset(cls, dataset) -> Emulsion:
         # there are values, so the emulsion is not empty
         droplet_class = dataset.attrs["droplet_class"]
         if droplet_class == "None":
-            droplets: List[SphericalDroplet] = []
+            droplets: list[SphericalDroplet] = []
         else:
             droplets = [
                 droplet_from_data(droplet_class, data)  # type: ignore
@@ -416,7 +410,7 @@ def to_file(self, path: str) -> None:
             self._write_hdf_dataset(fp)
 
     @property
-    def interface_width(self) -> Optional[float]:
+    def interface_width(self) -> float | None:
         """float: the average interface width across all droplets
 
         This averages the interface widths of the individual droplets weighted by their
@@ -446,9 +440,7 @@ def bbox(self) -> Cuboid:
             raise RuntimeError("Bounding box of empty emulsion is undefined")
         return sum((droplet.bbox for droplet in self[1:]), self[0].bbox)
 
-    def get_phasefield(
-        self, grid: GridBase, label: Optional[str] = None
-    ) -> ScalarField:
+    def get_phasefield(self, grid: GridBase, label: str | None = None) -> ScalarField:
         """create a phase field representing a list of droplets
 
         Args:
@@ -487,7 +479,7 @@ def remove_small(self, min_radius: float = -np.inf) -> None:
                 self.pop(i)
 
     def get_pairwise_distances(
-        self, subtract_radius: bool = False, grid: Optional[GridBase] = None
+        self, subtract_radius: bool = False, grid: GridBase | None = None
     ) -> np.ndarray:
         """return the pairwise distance between droplets
 
@@ -566,7 +558,7 @@ def get_neighbor_distances(self, subtract_radius: bool = False) -> np.ndarray:
             return dist[:, 1]  # type: ignore
 
     def remove_overlapping(
-        self, min_distance: float = 0, grid: Optional[GridBase] = None
+        self, min_distance: float = 0, grid: GridBase | None = None
     ) -> None:
         """remove all droplets that are overlapping
 
@@ -607,7 +599,7 @@ def total_droplet_volume(self) -> float:
         """float: the total volume of all droplets"""
         return sum(droplet.volume for droplet in self)  # type: ignore
 
-    def get_size_statistics(self, incl_vanished: bool = True) -> Dict[str, float]:
+    def get_size_statistics(self, incl_vanished: bool = True) -> dict[str, float]:
         """determine size statistics of the current emulsion
 
         Args:
@@ -645,10 +637,10 @@ def get_size_statistics(self, incl_vanished: bool = True) -> Dict[str, float]:
     def plot(
         self,
         ax,
-        field: Optional[ScalarField] = None,
-        image_args: Optional[Dict[str, Any]] = None,
+        field: ScalarField | None = None,
+        image_args: dict[str, Any] | None = None,
         repeat_periodically: bool = True,
-        color_value: Optional[Callable] = None,
+        color_value: Callable | None = None,
         cmap=None,
         norm=None,
         colorbar: bool | str = True,
@@ -740,7 +732,7 @@ def plot(
 
             # and map them to colors
             mapper = plt.cm.ScalarMappable(norm=norm, cmap=cmap)
-            colors: Union[List, np.ndarray] = mapper.to_rgba(values)
+            colors: list | np.ndarray = mapper.to_rgba(values)
 
             if kwargs.pop("color", None) is not None:
                 logger = logging.getLogger(self.__class__.__name__)
@@ -795,7 +787,7 @@ class EmulsionTimeCourse:
 
     def __init__(
         self,
-        emulsions: Optional[Iterable[Emulsion]] = None,
+        emulsions: Iterable[Emulsion] | None = None,
         times: np.ndarray | Sequence[float] | None = None,
     ) -> None:
         """
@@ -808,8 +800,8 @@ def __init__(
             times = emulsions.times
             emulsions = emulsions.emulsions
 
-        self.emulsions: List[Emulsion] = []
-        self.times: List[float] = []
+        self.emulsions: list[Emulsion] = []
+        self.times: list[float] = []
 
         # add all emulsions
         if emulsions is not None:
@@ -826,7 +818,7 @@ def __init__(
             raise ValueError("Lists of emulsions and times must have same length")
 
     def append(
-        self, emulsion: Emulsion, time: Optional[float] = None, copy: bool = True
+        self, emulsion: Emulsion, time: float | None = None, copy: bool = True
     ) -> None:
         """add an emulsion to the list
 
@@ -875,7 +867,7 @@ def __iter__(self) -> Iterator[Emulsion]:
         """iterate over the emulsions"""
         return iter(self.emulsions)
 
-    def items(self) -> Iterator[Tuple[float, Emulsion]]:
+    def items(self) -> Iterator[tuple[float, Emulsion]]:
         """iterate over all times and emulsions, returning them in pairs"""
         return zip(self.times, self.emulsions)
 
@@ -890,9 +882,9 @@ def from_storage(
         *,
         num_processes: int | Literal["auto"] = 1,
         refine: bool = False,
-        progress: Optional[bool] = None,
+        progress: bool | None = None,
         **kwargs,
-    ) -> "EmulsionTimeCourse":
+    ) -> EmulsionTimeCourse:
         r"""create an emulsion time course from a stored phase field
 
         Args:
@@ -943,7 +935,7 @@ def from_storage(
         return cls(emulsions, times=storage.times)
 
     @classmethod
-    def from_file(cls, path: str, progress: bool = True) -> "EmulsionTimeCourse":
+    def from_file(cls, path: str, progress: bool = True) -> EmulsionTimeCourse:
         """create emulsion time course by reading file
 
         Args:
@@ -970,7 +962,7 @@ def from_file(cls, path: str, progress: bool = True) -> "EmulsionTimeCourse":
                 )
         return obj
 
-    def to_file(self, path: str, info: Optional[InfoDict] = None) -> None:
+    def to_file(self, path: str, info: InfoDict | None = None) -> None:
         """store data in hdf5 file
 
         The data can be read using the classmethod :meth:`EmulsionTimeCourse.from_file`.
@@ -1010,10 +1002,10 @@ def get_emulsion(self, time: float) -> Emulsion:
     def tracker(
         self,
         interrupts: InterruptData = 1,
-        filename: Optional[str] = None,
+        filename: str | None = None,
         *,
         interval=None,
-    ) -> "DropletTracker":
+    ) -> DropletTracker:
         """return a tracker that analyzes emulsions during simulations
 
         Args:
diff --git a/droplets/image_analysis.py b/droplets/image_analysis.py
index a7e372c..add8fe9 100644
--- a/droplets/image_analysis.py
+++ b/droplets/image_analysis.py
@@ -22,17 +22,7 @@
 import warnings
 from functools import reduce
 from itertools import product
-from typing import (
-    Any,
-    Callable,
-    Dict,
-    Iterable,
-    List,
-    Literal,
-    Optional,
-    Sequence,
-    Tuple,
-)
+from typing import Any, Callable, Iterable, Literal, Sequence
 
 import numpy as np
 from numpy.lib.recfunctions import (
@@ -134,8 +124,8 @@ def _locate_droplets_in_mask_cartesian(mask: ScalarField) -> Emulsion:
     # connect clusters linked viaperiodic boundary conditions
     for ax in np.flatnonzero(grid.periodic):  # look at all periodic axes
         # compile list of all boundary points connected along the current axis
-        low: List[List[int] | np.ndarray] = []
-        high: List[List[int] | np.ndarray] = []
+        low: list[list[int] | np.ndarray] = []
+        high: list[list[int] | np.ndarray] = []
         for a in range(grid.num_axes):
             if a == ax:
                 low.append([0])
@@ -362,9 +352,9 @@ def locate_droplets(
     *,
     minimal_radius: float = 0,
     modes: int = 0,
-    interface_width: Optional[float] = None,
+    interface_width: float | None = None,
     refine: bool = False,
-    refine_args: Optional[Dict[str, Any]] = None,
+    refine_args: dict[str, Any] | None = None,
     num_processes: int | Literal["auto"] = 1,
 ) -> Emulsion:
     """Locates droplets in the phase field
@@ -455,7 +445,7 @@ def locate_droplets(
     for droplet in candidates:
         # check whether we need to add the interface width
         droplet_class = droplet.__class__
-        args: Dict[str, NumberOrArray] = {}
+        args: dict[str, NumberOrArray] = {}
 
         # change droplet class when interface width is given
         if interface_width is not None:
@@ -500,7 +490,7 @@ def refine_droplets(
     *,
     num_processes: int | Literal["auto"] = 1,
     **kwargs,
-) -> List[DiffuseDroplet]:
+) -> list[DiffuseDroplet]:
     r"""Refines many droplets by fitting to phase field
 
     Args:
@@ -553,8 +543,8 @@ def refine_droplet(
     vmin: float = 0.0,
     vmax: float = 1.0,
     adjust_values: bool = False,
-    tolerance: Optional[float] = None,
-    least_squares_params: Optional[Dict[str, Any]] = None,
+    tolerance: float | None = None,
+    least_squares_params: dict[str, Any] | None = None,
 ) -> DiffuseDroplet:
     """Refines droplet parameters by fitting to phase field
 
@@ -682,7 +672,7 @@ def get_structure_factor(
     smoothing: None | float | Literal["auto", "none"] = "auto",
     wave_numbers: Sequence[float] | Literal["auto"] = "auto",
     add_zero: bool = False,
-) -> Tuple[np.ndarray, np.ndarray]:
+) -> tuple[np.ndarray, np.ndarray]:
     r"""Calculates the structure factor associated with a field
 
     Here, the structure factor is basically the power spectral density of the field
@@ -791,7 +781,7 @@ def get_length_scale(
         "structure_factor_mean", "structure_factor_maximum", "droplet_detection"
     ] = "structure_factor_maximum",
     **kwargs,
-) -> float | Tuple[float, Any]:
+) -> float | tuple[float, Any]:
     """Calculates a length scale associated with a phase field
 
     Args:
diff --git a/droplets/resources/make_spheres_3d.py b/droplets/resources/make_spheres_3d.py
index 66c8652..56c24d7 100755
--- a/droplets/resources/make_spheres_3d.py
+++ b/droplets/resources/make_spheres_3d.py
@@ -1,6 +1,8 @@
 #!/usr/bin/env python3
 """script that creates the resource with discretized surfaces of unit spheres """
 
+from __future__ import annotations
+
 import h5py
 import numpy as np
 import pygmsh
diff --git a/droplets/tools/spherical.py b/droplets/tools/spherical.py
index a4f3288..7b4f236 100644
--- a/droplets/tools/spherical.py
+++ b/droplets/tools/spherical.py
@@ -45,7 +45,7 @@
 We will use these indices interchangeably, although the mode :math:`k` is
 preferred internally. Note that we also consider axisymmetric spherical
 harmonics, where the order is always zero and the degree :math:`l` and the mode
-:math:`k` are thus identical.    
+:math:`k` are thus identical.
 
 
 .. autosummary::
@@ -68,11 +68,12 @@
    spherical_harmonic_real
    spherical_harmonic_real_k
 
-
 .. codeauthor:: David Zwicker <david.zwicker@ds.mpg.de>    
 """
 
-from typing import Callable, Tuple, TypeVar
+from __future__ import annotations
+
+from typing import Callable, TypeVar
 
 import numpy as np
 from numba.extending import overload, register_jitable
@@ -369,7 +370,7 @@ def spherical_index_k(degree: int, order: int = 0) -> int:
     return degree * (degree + 1) + order
 
 
-def spherical_index_lm(k: int) -> Tuple[int, int]:
+def spherical_index_lm(k: int) -> tuple[int, int]:
     """returns the degree `l` and the order `m` from the mode `k`
 
     Args:
diff --git a/droplets/trackers.py b/droplets/trackers.py
index 1eb7872..61ca018 100644
--- a/droplets/trackers.py
+++ b/droplets/trackers.py
@@ -13,7 +13,7 @@
 from __future__ import annotations
 
 import math
-from typing import Any, Callable, Dict, List, Literal, Optional
+from typing import Any, Callable, Literal
 
 from pde.fields.base import FieldBase
 from pde.tools.docstrings import fill_in_docstring
@@ -36,7 +36,7 @@ class LengthScaleTracker(TrackerBase):
     def __init__(
         self,
         interrupts: InterruptData = 1,
-        filename: Optional[str] = None,
+        filename: str | None = None,
         *,
         method: Literal[
             "structure_factor_mean", "structure_factor_maximum", "droplet_detection"
@@ -67,8 +67,8 @@ def __init__(
                 Determines whether errors in determining the length scales are logged.
         """
         super().__init__(interrupts=interrupts, interval=interval)
-        self.length_scales: List[float] = []
-        self.times: List[float] = []
+        self.length_scales: list[float] = []
+        self.times: list[float] = []
         self.filename = filename
         self.method = method
         self.source = source
@@ -101,7 +101,7 @@ def handle(self, field: FieldBase, t: float):
         self.times.append(t)
         self.length_scales.append(length)  # type: ignore
 
-    def finalize(self, info: Optional[InfoDict] = None) -> None:
+    def finalize(self, info: InfoDict | None = None) -> None:
         """finalize the tracker, supplying additional information
 
         Args:
@@ -132,14 +132,14 @@ class DropletTracker(TrackerBase):
     def __init__(
         self,
         interrupts: InterruptData = 1,
-        filename: Optional[str] = None,
+        filename: str | None = None,
         *,
         emulsion_timecourse=None,
         source: None | int | Callable = None,
         threshold: float | Literal["auto", "extrema", "mean", "otsu"] = 0.5,
         minimal_radius: float = 0,
         refine: bool = False,
-        refine_args: Optional[Dict[str, Any]] = None,
+        refine_args: dict[str, Any] | None = None,
         perturbation_modes: int = 0,
         interval=None,
     ):
@@ -234,7 +234,7 @@ def handle(self, field: FieldBase, t: float) -> None:
         )
         self.data.append(emulsion, t)
 
-    def finalize(self, info: Optional[InfoDict] = None) -> None:
+    def finalize(self, info: InfoDict | None = None) -> None:
         """finalize the tracker, supplying additional information
 
         Args:
diff --git a/examples/tutorial/Using the py-droplets package.ipynb b/examples/tutorial/Using the py-droplets package.ipynb
index 09fda92..ba15fe1 100644
--- a/examples/tutorial/Using the py-droplets package.ipynb	
+++ b/examples/tutorial/Using the py-droplets package.ipynb	
@@ -107,7 +107,7 @@
     }
    ],
    "source": [
-    "drop2 = droplets.DiffuseDroplet(position=[6, 8], radius=2, interface_width=.3)\n",
+    "drop2 = droplets.DiffuseDroplet(position=[6, 8], radius=2, interface_width=0.3)\n",
     "drop2.volume"
    ]
   },
@@ -325,7 +325,7 @@
     }
    ],
    "source": [
-    "emulsion.get_phasefield().plot(title='An emulsion');"
+    "emulsion.get_phasefield().plot(title=\"An emulsion\");"
    ]
   },
   {
@@ -351,10 +351,13 @@
     }
    ],
    "source": [
-    "data = [droplets.SphericalDroplet(position=np.random.uniform(0, 100, 2),\n",
-    "                                  radius=np.random.uniform(5, 10))\n",
-    "        for _ in range(30)]\n",
-    "    \n",
+    "data = [\n",
+    "    droplets.SphericalDroplet(\n",
+    "        position=np.random.uniform(0, 100, 2), radius=np.random.uniform(5, 10)\n",
+    "    )\n",
+    "    for _ in range(30)\n",
+    "]\n",
+    "\n",
     "emulsion = droplets.Emulsion(data)\n",
     "emulsion.remove_overlapping()\n",
     "emulsion.plot(title=f\"{len(emulsion)} droplets\");"
@@ -383,7 +386,7 @@
     }
    ],
    "source": [
-    "emulsion.get_phasefield(pde.UnitGrid([100, 100])).plot(title='An emulsion');"
+    "emulsion.get_phasefield(pde.UnitGrid([100, 100])).plot(title=\"An emulsion\");"
    ]
   },
   {
@@ -424,7 +427,7 @@
    "source": [
     "# run a numerical simulation\n",
     "grid = pde.UnitGrid([64, 64], periodic=True)\n",
-    "field = pde.ScalarField.random_uniform(grid, -1, 0)        \n",
+    "field = pde.ScalarField.random_uniform(grid, -1, 0)\n",
     "eq = pde.CahnHilliardPDE()\n",
     "\n",
     "final = eq.solve(field, t_range=1e2, dt=0.01);"
@@ -626,7 +629,7 @@
    ],
    "source": [
     "tracker = droplets.DropletTracker(interrupts=50)\n",
-    "final = eq.solve(field, t_range=1e3, dt=0.01, tracker=['progress', tracker]);"
+    "final = eq.solve(field, t_range=1e3, dt=0.01, tracker=[\"progress\", tracker]);"
    ]
   },
   {
diff --git a/scripts/format_code.sh b/scripts/format_code.sh
index 494082e..2d8651a 100755
--- a/scripts/format_code.sh
+++ b/scripts/format_code.sh
@@ -1,6 +1,11 @@
 #!/usr/bin/env bash
 # This script formats the code of this package
 
+echo "Upgrading python syntax..."
+pushd .. > /dev/null
+find . -name '*.py' -exec pyupgrade --py38-plus {} +
+popd > /dev/null
+
 echo "Formating import statements..."
 isort ..
 
diff --git a/tests/requirements.txt b/tests/requirements.txt
index 86b202f..283f9a7 100644
--- a/tests/requirements.txt
+++ b/tests/requirements.txt
@@ -2,6 +2,7 @@
 black>=19.*
 h5py
 isort>=5.1
+pyupgrade>=3
 pytest>=5.4
 pytest-cov>=2.8
 pytest-xdist>=1.30
diff --git a/tests/test_examples.py b/tests/test_examples.py
index d194ae5..11c83c0 100644
--- a/tests/test_examples.py
+++ b/tests/test_examples.py
@@ -40,7 +40,7 @@ def test_example(path):
 
     msg = "Script `%s` failed with following output:" % path
     if outs:
-        msg = "%s\nSTDOUT:\n%s" % (msg, outs)
+        msg = f"{msg}\nSTDOUT:\n{outs}"
     if errs:
-        msg = "%s\nSTDERR:\n%s" % (msg, errs)
+        msg = f"{msg}\nSTDERR:\n{errs}"
     assert proc.returncode <= 0, msg