diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0357fb39e..c5852a528 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,6 +10,7 @@
 
 - Implement `split` method for the `Solenoid` element (see #380) (@cr-xu)
 - Implement a more robust RPN parser, fixing a bug where short strings in an Elegant variable definition would cause parsing to fail. (see #387) (@amylizzle, @Hespe, @jank324)
+- Add a `Sextupole` element (see #406) (@jank324, @Hespe)
 
 ### 🐛 Bug fixes
 
diff --git a/cheetah/__init__.py b/cheetah/__init__.py
index 5a37aa789..02d936a87 100644
--- a/cheetah/__init__.py
+++ b/cheetah/__init__.py
@@ -13,6 +13,7 @@
     RBend,
     Screen,
     Segment,
+    Sextupole,
     Solenoid,
     SpaceChargeKick,
     TransverseDeflectingCavity,
diff --git a/cheetah/accelerator/__init__.py b/cheetah/accelerator/__init__.py
index bc783dcee..7a8320e1a 100644
--- a/cheetah/accelerator/__init__.py
+++ b/cheetah/accelerator/__init__.py
@@ -11,6 +11,7 @@
 from .rbend import RBend  # noqa: F401
 from .screen import Screen  # noqa: F401
 from .segment import Segment  # noqa: F401
+from .sextupole import Sextupole  # noqa: F401
 from .solenoid import Solenoid  # noqa: F401
 from .space_charge_kick import SpaceChargeKick  # noqa: F401
 from .transverse_deflecting_cavity import TransverseDeflectingCavity  # noqa: F401
diff --git a/cheetah/accelerator/cavity.py b/cheetah/accelerator/cavity.py
index 83e901de2..6e34893c5 100644
--- a/cheetah/accelerator/cavity.py
+++ b/cheetah/accelerator/cavity.py
@@ -111,12 +111,10 @@ def _track_beam(self, incoming: Beam) -> Beam:
 
         tm = self.transfer_map(incoming.energy, incoming.species)
         if isinstance(incoming, ParameterBeam):
-            outgoing_mu = torch.matmul(tm, incoming.mu.unsqueeze(-1)).squeeze(-1)
-            outgoing_cov = torch.matmul(
-                tm, torch.matmul(incoming.cov, tm.transpose(-2, -1))
-            )
+            outgoing_mu = (tm @ incoming.mu.unsqueeze(-1)).squeeze(-1)
+            outgoing_cov = tm @ incoming.cov @ tm.transpose(-2, -1)
         else:  # ParticleBeam
-            outgoing_particles = torch.matmul(incoming.particles, tm.transpose(-2, -1))
+            outgoing_particles = incoming.particles @ tm.transpose(-2, -1)
         delta_energy = (
             self.voltage * torch.cos(phi) * incoming.species.num_elementary_charges * -1
         )
diff --git a/cheetah/accelerator/custom_transfer_map.py b/cheetah/accelerator/custom_transfer_map.py
index 355662fd7..a861e273c 100644
--- a/cheetah/accelerator/custom_transfer_map.py
+++ b/cheetah/accelerator/custom_transfer_map.py
@@ -69,10 +69,7 @@ def from_merging_elements(
             (*incoming_beam.energy.shape, 1, 1)
         )
         for element in elements:
-            tm = torch.matmul(
-                element.transfer_map(incoming_beam.energy, incoming_beam.species),
-                tm,
-            )
+            tm = element.transfer_map(incoming_beam.energy, incoming_beam.species) @ tm
             incoming_beam = element.track(incoming_beam)
 
         combined_length = sum(element.length for element in elements)
diff --git a/cheetah/accelerator/dipole.py b/cheetah/accelerator/dipole.py
index d14af9876..2e89ae483 100644
--- a/cheetah/accelerator/dipole.py
+++ b/cheetah/accelerator/dipole.py
@@ -407,11 +407,11 @@ def transfer_map(self, energy: torch.Tensor, species: Species) -> torch.Tensor:
             R[..., 2, 3] = self.length
 
         # Apply fringe fields
-        R = torch.matmul(R_exit, torch.matmul(R, R_enter))
+        R = R_exit @ R @ R_enter
+
         # Apply rotation for tilted magnets
-        R = torch.matmul(
-            rotation_matrix(-self.tilt), torch.matmul(R, rotation_matrix(self.tilt))
-        )
+        R = rotation_matrix(-self.tilt) @ R @ rotation_matrix(self.tilt)
+
         return R
 
     def _transfer_map_enter(self) -> torch.Tensor:
diff --git a/cheetah/accelerator/drift.py b/cheetah/accelerator/drift.py
index cf72e1dc3..5bf701402 100644
--- a/cheetah/accelerator/drift.py
+++ b/cheetah/accelerator/drift.py
@@ -5,7 +5,12 @@
 
 from cheetah.accelerator.element import Element
 from cheetah.particles import Beam, ParticleBeam, Species
-from cheetah.utils import UniqueNameGenerator, bmadx, compute_relativistic_factors
+from cheetah.utils import (
+    UniqueNameGenerator,
+    bmadx,
+    compute_relativistic_factors,
+    verify_device_and_dtype,
+)
 
 generate_unique_name = UniqueNameGenerator(prefix="unnamed_element")
 
@@ -30,6 +35,7 @@ def __init__(
         device: torch.device | None = None,
         dtype: torch.dtype | None = None,
     ) -> None:
+        device, dtype = verify_device_and_dtype([length], device, dtype)
         factory_kwargs = {"device": device, "dtype": dtype}
         super().__init__(name=name, **factory_kwargs)
 
diff --git a/cheetah/accelerator/element.py b/cheetah/accelerator/element.py
index 6fb656dc4..e90a094a8 100644
--- a/cheetah/accelerator/element.py
+++ b/cheetah/accelerator/element.py
@@ -68,8 +68,8 @@ def track(self, incoming: Beam) -> Beam:
         """
         if isinstance(incoming, ParameterBeam):
             tm = self.transfer_map(incoming.energy, incoming.species)
-            mu = torch.matmul(tm, incoming.mu.unsqueeze(-1)).squeeze(-1)
-            cov = torch.matmul(tm, torch.matmul(incoming.cov, tm.transpose(-2, -1)))
+            mu = (tm @ incoming.mu.unsqueeze(-1)).squeeze(-1)
+            cov = tm @ incoming.cov @ tm.transpose(-2, -1)
             return ParameterBeam(
                 mu,
                 cov,
@@ -79,7 +79,7 @@ def track(self, incoming: Beam) -> Beam:
             )
         elif isinstance(incoming, ParticleBeam):
             tm = self.transfer_map(incoming.energy, incoming.species)
-            new_particles = torch.matmul(incoming.particles, tm.transpose(-2, -1))
+            new_particles = incoming.particles @ tm.transpose(-2, -1)
             return ParticleBeam(
                 new_particles,
                 incoming.energy,
diff --git a/cheetah/accelerator/quadrupole.py b/cheetah/accelerator/quadrupole.py
index 3c58a0557..413a3047e 100644
--- a/cheetah/accelerator/quadrupole.py
+++ b/cheetah/accelerator/quadrupole.py
@@ -19,7 +19,7 @@ class Quadrupole(Element):
     :param length: Length in meters.
     :param k1: Strength of the quadrupole in 1/m^-2.
     :param misalignment: Misalignment vector of the quadrupole in x- and y-directions.
-    :param tilt: Tilt angle of the quadrupole in x-y plane [rad]. pi/4 for
+    :param tilt: Tilt angle of the quadrupole in x-y plane in radians. pi/4 for
         skew-quadrupole.
     :param num_steps: Number of drift-kick-drift steps to use for tracking through the
         element when tracking method is set to `"bmadx"`.
diff --git a/cheetah/accelerator/segment.py b/cheetah/accelerator/segment.py
index 916bdb016..7877c0b1d 100644
--- a/cheetah/accelerator/segment.py
+++ b/cheetah/accelerator/segment.py
@@ -399,7 +399,7 @@ def transfer_map(self, energy: torch.Tensor, species: Species) -> torch.Tensor:
         if self.is_skippable:
             tm = torch.eye(7, device=energy.device, dtype=energy.dtype)
             for element in self.elements:
-                tm = torch.matmul(element.transfer_map(energy, species), tm)
+                tm = element.transfer_map(energy, species) @ tm
             return tm
         else:
             return None
diff --git a/cheetah/accelerator/sextupole.py b/cheetah/accelerator/sextupole.py
new file mode 100644
index 000000000..a3b08ec2d
--- /dev/null
+++ b/cheetah/accelerator/sextupole.py
@@ -0,0 +1,142 @@
+import matplotlib.pyplot as plt
+import torch
+
+from cheetah.accelerator.element import Element
+from cheetah.particles import Beam, ParameterBeam, ParticleBeam, Species
+from cheetah.track_methods import base_rmatrix, base_ttensor, misalignment_matrix
+from cheetah.utils import verify_device_and_dtype
+
+
+class Sextupole(Element):
+    """
+    A sextupole element in a particle accelerator.
+
+    :param length: Length in meters.
+    :param k2: Sextupole strength in 1/m^3.
+    :param misalignment: Transverse misalignment in x and y directions in meters.
+    :param tilt: Tilt angle of the quadrupole in x-y plane in radians.
+    :param name: Unique identifier of the element.
+    """
+
+    def __init__(
+        self,
+        length: torch.Tensor,
+        k2: torch.Tensor | None = None,
+        misalignment: torch.Tensor | None = None,
+        tilt: torch.Tensor | None = None,
+        name: str | None = None,
+        device: torch.device | None = None,
+        dtype: torch.dtype | None = None,
+    ) -> None:
+        device, dtype = verify_device_and_dtype(
+            [length, k2, misalignment, tilt], device, dtype
+        )
+        factory_kwargs = {"device": device, "dtype": dtype}
+        super().__init__(name=name, **factory_kwargs)
+
+        self.length = torch.as_tensor(length, **factory_kwargs)
+
+        self.register_buffer_or_parameter(
+            "k2", torch.as_tensor(k2 if k2 is not None else 0.0, **factory_kwargs)
+        )
+        self.register_buffer_or_parameter(
+            "misalignment",
+            torch.as_tensor(
+                misalignment if misalignment is not None else (0.0, 0.0),
+                **factory_kwargs,
+            ),
+        )
+        self.register_buffer_or_parameter(
+            "tilt", torch.as_tensor(tilt if tilt is not None else 0.0, **factory_kwargs)
+        )
+
+    def transfer_map(self, energy: torch.Tensor, species: Species) -> torch.Tensor:
+        R = base_rmatrix(
+            length=self.length,
+            k1=torch.zeros_like(self.length),
+            hx=torch.zeros_like(self.length),
+            species=species,
+            tilt=self.tilt,
+            energy=energy,
+        )
+
+        if torch.all(self.misalignment == 0):
+            return R
+        else:
+            R_entry, R_exit = misalignment_matrix(self.misalignment)
+            R = R_exit @ R @ R_entry
+            return R
+
+    def track(self, incoming: Beam) -> Beam:
+        """
+        Track the beam through the sextupole element.
+
+        :param incoming: Beam entering the element.
+        :return: Beam exiting the element.
+        """
+        first_order_tm = self.transfer_map(incoming.energy, incoming.species)
+        second_order_tm = base_ttensor(
+            length=self.length,
+            k1=torch.zeros_like(self.length),
+            k2=self.k2,
+            hx=torch.zeros_like(self.length),
+            species=incoming.species,
+            tilt=self.tilt,
+            energy=incoming.energy,
+        )
+
+        if isinstance(incoming, ParameterBeam):
+            # For ParameterBeam, only first-order effects are applied
+            return super().track(incoming)
+        elif isinstance(incoming, ParticleBeam):
+            # Apply the transfer map to the incoming particles
+            first_order_particles = incoming.particles @ first_order_tm.transpose(
+                -2, -1
+            )
+            second_order_particles = torch.einsum(
+                "...ijk,...j,...k->...i",
+                second_order_tm.unsqueeze(-4),  # Add broadcast dimension for particles
+                incoming.particles,
+                incoming.particles,
+            )
+            outgoing_particles = second_order_particles + first_order_particles
+
+            return ParticleBeam(
+                particles=outgoing_particles,
+                energy=incoming.energy,
+                particle_charges=incoming.particle_charges,
+                survival_probabilities=incoming.survival_probabilities,
+                species=incoming.species,
+            )
+        else:
+            raise TypeError(
+                f"Unsupported beam type: {type(incoming)}. Expected ParameterBeam or "
+                "ParticleBeam."
+            )
+
+    @property
+    def is_skippable(self) -> bool:
+        return False
+
+    @property
+    def is_active(self) -> bool:
+        return torch.any(self.k2 != 0.0).item()
+
+    def split(self, resolution: torch.Tensor) -> list[Element]:
+        raise NotImplementedError
+
+    def plot(self, ax: plt.Axes, s: float, vector_idx: tuple | None = None) -> None:
+        raise NotImplementedError
+
+    @property
+    def defining_features(self) -> list[str]:
+        return super().defining_features + ["length", "k2", "misalignment", "tilt"]
+
+    def __repr__(self) -> str:
+        return (
+            f"{self.__class__.__name__}(length={repr(self.length)}, "
+            f"k2={repr(self.k2)}, "
+            f"misalignment={repr(self.misalignment)}, "
+            f"tilt={repr(self.tilt)}, "
+            f"name={repr(self.name)})"
+        )
diff --git a/cheetah/track_methods.py b/cheetah/track_methods.py
index 1c687c4de..8956a540e 100644
--- a/cheetah/track_methods.py
+++ b/cheetah/track_methods.py
@@ -6,29 +6,6 @@
 from cheetah.utils import compute_relativistic_factors
 
 
-def rotation_matrix(angle: torch.Tensor) -> torch.Tensor:
-    """Rotate the transfer map in x-y plane.
-
-    :param angle: Rotation angle in rad, for example `angle = np.pi/2` for vertical =
-        dipole.
-    :return: Rotation matrix to be multiplied to the element's transfer matrix.
-    """
-    cs = torch.cos(angle)
-    sn = torch.sin(angle)
-
-    tm = torch.eye(7, dtype=angle.dtype, device=angle.device).repeat(*angle.shape, 1, 1)
-    tm[..., 0, 0] = cs
-    tm[..., 0, 2] = sn
-    tm[..., 1, 1] = cs
-    tm[..., 1, 3] = sn
-    tm[..., 2, 0] = -sn
-    tm[..., 2, 2] = cs
-    tm[..., 3, 1] = -sn
-    tm[..., 3, 3] = cs
-
-    return tm
-
-
 def base_rmatrix(
     length: torch.Tensor,
     k1: torch.Tensor,
@@ -38,15 +15,15 @@ def base_rmatrix(
     energy: torch.Tensor | None = None,
 ) -> torch.Tensor:
     """
-    Create a universal transfer matrix for a beamline element.
+    Create a first order universal transfer map for a beamline element.
 
     :param length: Length of the element in m.
     :param k1: Quadrupole strength in 1/m**2.
-    :param hx: Curvature (1/radius) of the element in 1/m**2.
+    :param hx: Curvature (1/radius) of the element in 1/m.
     :param species: Particle species of the beam.
     :param tilt: Roation of the element relative to the longitudinal axis in rad.
     :param energy: Beam energy in eV.
-    :return: Transfer matrix for the element.
+    :return: First order transfer map for the element.
     """
     device = length.device
     dtype = length.dtype
@@ -96,12 +73,240 @@ def base_rmatrix(
 
     # Rotate the R matrix for skew / vertical magnets
     if torch.any(tilt != 0):
-        R = torch.einsum(
-            "...ij,...jk,...kl->...il", rotation_matrix(-tilt), R, rotation_matrix(tilt)
-        )
+        R = rotation_matrix(-tilt) @ R @ rotation_matrix(tilt)
+
     return R
 
 
+def base_ttensor(
+    length: torch.Tensor,
+    k1: torch.Tensor,
+    k2: torch.Tensor,
+    hx: torch.Tensor,
+    species: Species,
+    tilt: torch.Tensor | None = None,
+    energy: torch.Tensor | None = None,
+) -> torch.Tensor:
+    """
+    Create a second order universal transfer map for a beamline element. Uses MAD
+    convention.
+
+    :param length: Length of the element in m.
+    :param k1: Quadrupole strength in 1/m**2.
+    :param k2: Sextupole strength in 1/m**3.
+    :param hx: Curvature (1/radius) of the element in 1/m.
+    :param species: Particle species of the beam.
+    :param tilt: Roation of the element relative to the longitudinal axis in rad.
+    :param energy: Beam energy in eV.
+    :return: Second order transfer map for the element.
+    """
+    device = length.device
+    dtype = length.dtype
+
+    tilt = tilt if tilt is not None else torch.tensor(0.0, device=device, dtype=dtype)
+    energy = (
+        energy if energy is not None else torch.tensor(0.0, device=device, dtype=dtype)
+    )
+
+    _, igamma2, beta = compute_relativistic_factors(energy, species.mass_eV)
+
+    kx2 = k1 + hx**2
+    ky2 = -k1
+    kx = torch.sqrt(torch.complex(kx2, torch.tensor(0.0, device=device, dtype=dtype)))
+    ky = torch.sqrt(torch.complex(ky2, torch.tensor(0.0, device=device, dtype=dtype)))
+    cx = torch.cos(kx * length).real
+    cy = torch.cos(ky * length).real
+    sx = torch.where(kx != 0, (torch.sin(kx * length) / kx).real, length)
+    sy = torch.where(ky != 0, (torch.sin(ky * length) / ky).real, length)
+    dx = torch.where(kx != 0, (1.0 - cx) / kx2, length**2 / 2.0)
+
+    d2y = 0.5 * sy**2
+    s2y = sy * cy
+    c2y = torch.cos(2 * ky * length).real
+    fx = torch.where(kx2 != 0, (length - sx) / kx2, length**3 / 6.0)
+    f2y = torch.where(ky2 != 0, (length - s2y) / ky2, length**3 / 6.0)
+
+    j1 = torch.where(kx2 != 0, (length - sx) / kx2, length**3 / 6.0)
+    j2 = torch.where(
+        kx2 != 0,
+        (3.0 * length - 4.0 * sx + sx * cx) / (2 * kx2**2),
+        length**5 / 20.0,
+    )
+    j3 = torch.where(
+        kx2 != 0,
+        (15.0 * length - 22.5 * sx + 9.0 * sx * cx - 1.5 * sx * cx**2 + kx2 * sx**3)
+        / (6.0 * kx2**3),
+        length**7 / 56.0,
+    )
+    j_denominator = kx2 - 4 * ky2
+    jc = torch.where(j_denominator != 0, (c2y - cx) / j_denominator, 0.5 * length**2)
+    js = torch.where(
+        j_denominator != 0, (cy * sy - sx) / j_denominator, length**3 / 6.0
+    )
+    jd = torch.where(j_denominator != 0, (d2y - dx) / j_denominator, length**4 / 24.0)
+    jf = torch.where(j_denominator != 0, (f2y - fx) / j_denominator, length**5 / 120.0)
+
+    khk = k2 + 2 * hx * k1
+
+    vector_shape = torch.broadcast_shapes(
+        length.shape, k1.shape, k2.shape, hx.shape, tilt.shape, energy.shape
+    )
+
+    T = torch.zeros((7, 7, 7), dtype=dtype, device=device).repeat(
+        *vector_shape, 1, 1, 1
+    )
+    T[..., 0, 0, 0] = -1 / 6 * khk * (sx**2 + dx) - 0.5 * hx * kx2 * sx**2
+    T[..., 0, 0, 1] = 2 * (-1 / 6 * khk * sx * dx + 0.5 * hx * sx * cx)
+    T[..., 0, 1, 1] = -1 / 6 * khk * dx**2 + 0.5 * hx * dx * cx
+    T[..., 0, 0, 5] = 2 * (
+        -hx / 12 / beta * khk * (3 * sx * j1 - dx**2)
+        + 0.5 * hx**2 / beta * sx**2
+        + 0.25 / beta * k1 * length * sx
+    )
+    T[..., 0, 1, 5] = 2 * (
+        -hx / 12 / beta * khk * (sx * dx**2 - 2 * cx * j2)
+        + 0.25 * hx**2 / beta * (sx * dx + cx * j1)
+        - 0.25 / beta * (sx + length * cx)
+    )
+    T[..., 0, 5, 5] = (
+        -(hx**2) / 6 / beta**2 * khk * (dx**2 * dx - 2 * sx * j2)
+        + 0.5 * hx**3 / beta**2 * sx * j1
+        - 0.5 * hx / beta**2 * length * sx
+        - 0.5 * hx / (beta**2) * igamma2 * dx
+    )
+    T[..., 0, 2, 2] = k1 * k2 * jd + 0.5 * (k2 + hx * k1) * dx
+    T[..., 0, 2, 3] = 2 * (0.5 * k2 * js)
+    T[..., 0, 3, 3] = k2 * jd - 0.5 * hx * dx
+    T[..., 1, 0, 0] = -1 / 6 * khk * sx * (1 + 2 * cx)
+    T[..., 1, 0, 1] = 2 * (-1 / 6 * khk * dx * (1 + 2 * cx))
+    T[..., 1, 1, 1] = -1 / 3 * khk * sx * dx - 0.5 * hx * sx
+    T[..., 1, 0, 5] = 2 * (
+        -hx / 12 / beta * khk * (3 * cx * j1 + sx * dx)
+        - 0.25 / beta * k1 * (sx - length * cx)
+    )
+    T[..., 1, 1, 5] = 2 * (
+        -hx / 12 / beta * khk * (3 * sx * j1 + dx**2) + 0.25 / beta * k1 * length * sx
+    )
+    T[..., 1, 5, 5] = (
+        -(hx**2) / 6 / beta**2 * khk * (sx * dx**2 - 2 * cx * j2)
+        - 0.5 * hx / beta**2 * k1 * (cx * j1 - sx * dx)
+        - 0.5 * hx / beta**2 * igamma2 * sx
+    )
+    T[..., 1, 2, 2] = k1 * k2 * js + 0.5 * (k2 + hx * k1) * sx
+    T[..., 1, 2, 3] = 2 * (0.5 * k2 * jc)
+    T[..., 1, 3, 3] = k2 * js - 0.5 * hx * sx
+    T[..., 2, 0, 2] = 2 * (
+        0.5 * k2 * (cy * jc - 2 * k1 * sy * js) + 0.5 * hx * k1 * sx * sy
+    )
+    T[..., 2, 0, 3] = 2 * (0.5 * k2 * (sy * jc - 2 * cy * js) + 0.5 * hx * sx * cy)
+    T[..., 2, 1, 2] = 2 * (
+        0.5 * k2 * (cy * js - 2 * k1 * sy * jd) + 0.5 * hx * k1 * dx * sy
+    )
+    T[..., 2, 1, 3] = 2 * (0.5 * k2 * (sy * js - 2 * cy * jd) + 0.5 * hx * dx * cy)
+    T[..., 2, 2, 5] = 2 * (
+        0.5 * hx / beta * k2 * (cy * jd - 2 * k1 * sy * jf)
+        + 0.5 * hx**2 / beta * k1 * j1 * sy
+        - 0.25 / beta * k1 * length * sy
+    )
+    T[..., 2, 3, 5] = 2 * (
+        0.5 * hx / beta * k2 * (sy * jd - 2 * cy * jf)
+        + 0.5 * hx**2 / beta * j1 * cy
+        - 0.25 / beta * (sy + length * cy)
+    )
+    T[..., 3, 0, 2] = 2 * (
+        0.5 * k1 * k2 * (2 * cy * js - sy * jc) + 0.5 * (k2 + hx * k1) * sx * cy
+    )
+    T[..., 3, 0, 3] = 2 * (
+        0.5 * k2 * (2 * k1 * sy * js - cy * jc) + 0.5 * (k2 + hx * k1) * sx * sy
+    )
+    T[..., 3, 1, 2] = 2 * (
+        0.5 * k1 * k2 * (2 * cy * jd - sy * js) + 0.5 * (k2 + hx * k1) * dx * cy
+    )
+    T[..., 3, 1, 3] = 2 * (
+        0.5 * k2 * (2 * k1 * sy * jd - cy * js) + 0.5 * (k2 + hx * k1) * dx * sy
+    )
+    T[..., 3, 2, 5] = 2 * (
+        0.5 * hx / beta * k1 * k2 * (2 * cy * jf - sy * jd)
+        + 0.5 * hx / beta * (k2 + hx * k1) * j1 * cy
+        + 0.25 / beta * k1 * (sy - length * cy)
+    )
+    T[..., 3, 3, 5] = 2 * (
+        0.5 * hx / beta * k2 * (2 * k1 * sy * jf - cy * jd)
+        + 0.5 * hx / beta * (k2 + hx * k1) * j1 * sy
+        - 0.25 / beta * k1 * length * sy
+    )
+    T[..., 4, 0, 0] = -1 * (
+        hx / 12 / beta * khk * (sx * dx + 3 * j1)
+        - 0.25 / beta * k1 * (length - sx * cx)
+    )
+    T[..., 4, 0, 1] = -2 * (hx / 12 / beta * khk * dx**2 + 0.25 / beta * k1 * sx**2)
+    T[..., 4, 1, 1] = -1 * (
+        hx / 6 / beta * khk * j2 - 0.5 / beta * sx - 0.25 / beta * k1 * (j1 - sx * dx)
+    )
+    T[..., 4, 0, 5] = -2 * (
+        hx**2 / 12 / beta**2 * khk * (3 * dx * j1 - 4 * j2)
+        + 0.25 * hx / beta**2 * k1 * j1 * (1 + cx)
+        + 0.5 * hx / beta**2 * igamma2 * sx
+    )
+    T[..., 4, 1, 5] = -2 * (
+        hx**2 / 12 / beta**2 * khk * (dx * dx**2 - 2 * sx * j2)
+        + 0.25 * hx / beta**2 * k1 * sx * j1
+        + 0.5 * hx / beta**2 * igamma2 * dx
+    )
+    T[..., 4, 5, 5] = -1 * (
+        hx**3 / 6 / beta**3 * khk * (3 * j3 - 2 * dx * j2)
+        + hx**2 / 6 / beta**3 * k1 * (sx * dx**2 - j2 * (1 + 2 * cx))
+        + 1.5 / beta**3 * igamma2 * (hx**2 * j1 - length)
+    )
+    T[..., 4, 2, 2] = -1 * (
+        -hx / beta * k1 * k2 * jf
+        - 0.5 * hx / beta * (k2 + hx * k1) * j1
+        + 0.25 / beta * k1 * (length - cy * sy)
+    )
+    T[..., 4, 2, 3] = -2 * (-0.5 * hx / beta * k2 * jd - 0.25 / beta * k1 * sy**2)
+    T[..., 4, 3, 3] = -1 * (
+        -hx / beta * k2 * jf
+        + 0.5 * hx**2 / beta * j1
+        - 0.25 / beta * (length + cy * sy)
+    )
+    T[..., 6, 6, 6] = 0.0  # Constant term currently handled by first order transfer map
+
+    # Rotate the R matrix for skew / vertical magnets
+    if torch.any(tilt != 0):
+        T = torch.einsum(
+            "...ij,...jkl,...kn,...lm->...inm",
+            rotation_matrix(-tilt),
+            T,
+            rotation_matrix(tilt),
+            rotation_matrix(tilt),
+        )
+    return T
+
+
+def rotation_matrix(angle: torch.Tensor) -> torch.Tensor:
+    """
+    Rotate the coordinate system in the x-y plane.
+
+    :param angle: Rotation angle in rad, for example `angle = np.pi/2` for vertical =
+        dipole.
+    :return: Rotation matrix to be multiplied to the element's transfer matrix.
+    """
+    cs = torch.cos(angle)
+    sn = torch.sin(angle)
+
+    tm = torch.eye(7, dtype=angle.dtype, device=angle.device).repeat(*angle.shape, 1, 1)
+    tm[..., 0, 0] = cs
+    tm[..., 0, 2] = sn
+    tm[..., 1, 1] = cs
+    tm[..., 1, 3] = sn
+    tm[..., 2, 0] = -sn
+    tm[..., 2, 2] = cs
+    tm[..., 3, 1] = -sn
+    tm[..., 3, 3] = cs
+
+    return tm
+
+
 def misalignment_matrix(
     misalignment: torch.Tensor,
 ) -> tuple[torch.Tensor, torch.Tensor]:
diff --git a/cheetah/utils/kde.py b/cheetah/utils/kde.py
index 7d5c82973..22904785b 100644
--- a/cheetah/utils/kde.py
+++ b/cheetah/utils/kde.py
@@ -101,7 +101,7 @@ def _kde_joint_pdf_2d(
             + f"Got {type(kernel_values2)}"
         )
 
-    joint_kernel_values = torch.matmul(kernel_values1.transpose(-2, -1), kernel_values2)
+    joint_kernel_values = kernel_values1.transpose(-2, -1) @ kernel_values2
     normalization = (
         torch.sum(joint_kernel_values, dim=(-2, -1)).unsqueeze(-1).unsqueeze(-1)
         + epsilon
diff --git a/docs/accelerator.rst b/docs/accelerator.rst
index 8cc0a737d..959728d07 100644
--- a/docs/accelerator.rst
+++ b/docs/accelerator.rst
@@ -55,6 +55,10 @@ Accelerator
     :members:
     :undoc-members:
 
+.. automodule:: accelerator.sextupole
+    :members:
+    :undoc-members:
+
 .. automodule:: accelerator.solenoid
     :members:
     :undoc-members:
diff --git a/setup.py b/setup.py
index 284fb8cad..b5c61a8e5 100644
--- a/setup.py
+++ b/setup.py
@@ -19,6 +19,6 @@
     long_description_content_type="text/markdown",
     packages=[package for package in find_packages() if package.startswith("cheetah")],
     python_requires=">=3.10",
-    install_requires=["matplotlib", "numpy", "scipy", "torch<=2.6"],
+    install_requires=["matplotlib", "numpy", "scipy", "torch[opt-einsum]<=2.6"],
     extras_require={"openpmd": ["openpmd-beamphysics"]},
 )
diff --git a/tests/test_clone.py b/tests/test_clone.py
index 7d12d0452..1c3382aec 100644
--- a/tests/test_clone.py
+++ b/tests/test_clone.py
@@ -13,6 +13,7 @@
         cheetah.HorizontalCorrector,
         cheetah.Quadrupole,
         cheetah.RBend,
+        cheetah.Sextupole,
         cheetah.Solenoid,
         cheetah.TransverseDeflectingCavity,
         cheetah.Undulator,
diff --git a/tests/test_device_dtype.py b/tests/test_device_dtype.py
index b293b6bc1..ba6868ccf 100644
--- a/tests/test_device_dtype.py
+++ b/tests/test_device_dtype.py
@@ -53,6 +53,7 @@ def test_move_quadrupole_to_device(target_device: torch.device):
         cheetah.HorizontalCorrector,
         cheetah.Quadrupole,
         cheetah.RBend,
+        cheetah.Sextupole,
         cheetah.Solenoid,
         cheetah.TransverseDeflectingCavity,
         cheetah.Undulator,
@@ -80,6 +81,7 @@ def test_forced_element_dtype(ElementClass):
         cheetah.HorizontalCorrector,
         cheetah.Quadrupole,
         cheetah.RBend,
+        cheetah.Sextupole,
         cheetah.Solenoid,
         cheetah.TransverseDeflectingCavity,
         cheetah.Undulator,
diff --git a/tests/test_elements.py b/tests/test_elements.py
index 6e5fb3992..771a8bb72 100644
--- a/tests/test_elements.py
+++ b/tests/test_elements.py
@@ -11,6 +11,7 @@
     cheetah.HorizontalCorrector: {"length": torch.tensor(1.0)},
     cheetah.Quadrupole: {"length": torch.tensor(1.0)},
     cheetah.Segment: {"elements": [cheetah.Drift(length=torch.tensor(1.0))]},
+    cheetah.Sextupole: {"length": torch.tensor(1.0)},
     cheetah.Solenoid: {"length": torch.tensor(1.0)},
     cheetah.SpaceChargeKick: {"effect_length": torch.tensor(1.0)},
     cheetah.TransverseDeflectingCavity: {"length": torch.tensor(1.0)},
diff --git a/tests/test_sextupole.py b/tests/test_sextupole.py
new file mode 100644
index 000000000..b8fb4d940
--- /dev/null
+++ b/tests/test_sextupole.py
@@ -0,0 +1,138 @@
+from copy import deepcopy
+
+import ocelot
+import torch
+
+import cheetah
+
+
+def test_compare_sextupole_to_ocelot_particle():
+    """
+    Compare the results of tracking through a sextupole in Cheetah and Ocelot. For a
+    `ParticleBeam` with second order effects in Ocelot.
+    """
+    length = 0.11
+    k2 = 87.0
+    tilt = torch.pi / 2
+
+    # Track through a sextupole in Cheetah
+    incoming = cheetah.ParticleBeam.from_astra(
+        "tests/resources/ACHIP_EA1_2021.1351.001"
+    )
+    cheetah_sextupole = cheetah.Sextupole(
+        length=torch.tensor(length), k2=torch.tensor(k2), tilt=torch.tensor(tilt)
+    )
+    outgoing_cheetah = cheetah_sextupole.track(incoming)
+
+    # Convert to Ocelot sextupole
+    incoming_p_array = ocelot.astraBeam2particleArray(
+        "tests/resources/ACHIP_EA1_2021.1351.001"
+    )
+    lattice = ocelot.MagneticLattice(
+        [ocelot.Sextupole(l=length, k2=k2, tilt=tilt)],
+        method={"global": ocelot.SecondTM},
+    )
+    navigator = ocelot.Navigator(lattice)
+    _, outgoing_p_array = ocelot.track(lattice, deepcopy(incoming_p_array), navigator)
+    outgoing_ocelot = cheetah.ParticleBeam.from_ocelot(outgoing_p_array)
+
+    # Compare the results
+    assert torch.allclose(
+        outgoing_cheetah.particles, outgoing_ocelot.particles, atol=1e-5, rtol=1e-6
+    )
+
+
+def test_compare_sextupole_to_ocelot_particle_vectorized():
+    """
+    Compare the results of tracking through a sextupole in Cheetah and Ocelot. For a
+    `ParticleBeam` with second order effects in Ocelot.
+
+    Vectorised version of the test.
+    """
+    length = 0.11
+    k2 = 87.0
+    tilt = torch.pi / 2
+
+    # Track through a sextupole in Cheetah
+    incoming = cheetah.ParticleBeam.from_astra(
+        "tests/resources/ACHIP_EA1_2021.1351.001"
+    )
+    cheetah_sextupole = cheetah.Sextupole(
+        length=torch.tensor(length),
+        k2=torch.tensor(k2).repeat([2]),
+        tilt=torch.tensor(tilt).repeat([3, 1]),
+    )
+    outgoing_cheetah = cheetah_sextupole.track(incoming)
+
+    # Convert to Ocelot sextupole
+    incoming_p_array = ocelot.astraBeam2particleArray(
+        "tests/resources/ACHIP_EA1_2021.1351.001"
+    )
+    lattice = ocelot.MagneticLattice(
+        [ocelot.Sextupole(l=length, k2=k2, tilt=tilt)],
+        method={"global": ocelot.SecondTM},
+    )
+    navigator = ocelot.Navigator(lattice)
+    _, outgoing_p_array = ocelot.track(lattice, deepcopy(incoming_p_array), navigator)
+    outgoing_ocelot = cheetah.ParticleBeam.from_ocelot(outgoing_p_array)
+
+    # Compare the results
+    assert torch.allclose(
+        outgoing_cheetah.particles, outgoing_ocelot.particles, atol=1e-5, rtol=1e-6
+    )
+
+
+def test_compare_sextupole_to_ocelot_parameter():
+    """
+    Compare the results of tracking through a sextupole in Cheetah and Ocelot for a
+    `ParameterBeam` with only first order effects in Ocelot.
+    """
+    length = 0.11
+    k2 = 87.0
+    tilt = torch.pi / 2
+
+    # Track through a sextupole in Cheetah
+    incoming = cheetah.ParameterBeam.from_astra(
+        "tests/resources/ACHIP_EA1_2021.1351.001"
+    )
+    cheetah_sextupole = cheetah.Sextupole(
+        length=torch.tensor(length), k2=torch.tensor(k2), tilt=torch.tensor(tilt)
+    )
+    outgoing_cheetah = cheetah_sextupole.track(incoming)
+
+    # Convert to Ocelot sextupole
+    incoming_p_array = ocelot.astraBeam2particleArray(
+        "tests/resources/ACHIP_EA1_2021.1351.001"
+    )
+    lattice = ocelot.MagneticLattice(
+        [ocelot.Sextupole(l=length, k2=k2, tilt=tilt)],
+        method={"global": ocelot.TransferMap},
+    )
+    navigator = ocelot.Navigator(lattice)
+    _, outgoing_p_array = ocelot.track(lattice, deepcopy(incoming_p_array), navigator)
+    outgoing_ocelot = cheetah.ParameterBeam.from_ocelot(outgoing_p_array)
+
+    # Compare the results
+    assert torch.allclose(outgoing_cheetah.mu, outgoing_ocelot.mu, atol=1e-5, rtol=1e-6)
+    assert torch.allclose(
+        outgoing_cheetah.cov, outgoing_ocelot.cov, atol=1e-5, rtol=1e-6
+    )
+
+
+def test_sextupole_as_drift():
+    """Test that a sextupole with k2=0 is equivalent to a drift."""
+    incoming = cheetah.ParticleBeam.from_astra(
+        "tests/resources/ACHIP_EA1_2021.1351.001"
+    )
+
+    sextupole = cheetah.Sextupole(length=torch.tensor(0.11), k2=torch.tensor(0.0))
+    drift = cheetah.Drift(length=torch.tensor(0.11))
+
+    # Track through the sextupole and drift
+    sextupole_outgoing = sextupole.track(incoming)
+    drift_outgoing = drift.track(incoming)
+
+    # Check that the results are the same
+    assert torch.allclose(
+        sextupole_outgoing.particles, drift_outgoing.particles, atol=1e-5, rtol=1e-6
+    )
diff --git a/tests/test_vectorized.py b/tests/test_vectorized.py
index 1458eae09..3b668ff44 100644
--- a/tests/test_vectorized.py
+++ b/tests/test_vectorized.py
@@ -342,6 +342,7 @@ def test_vectorized_screen_2d(BeamClass, method):
         cheetah.HorizontalCorrector,
         cheetah.Quadrupole,
         cheetah.RBend,
+        cheetah.Sextupole,
         cheetah.Solenoid,
         cheetah.TransverseDeflectingCavity,
         cheetah.Undulator,