diff --git a/benchmarks/benchmark_multipole_pyat.py b/benchmarks/benchmark_multipole_pyat.py
new file mode 100644
index 000000000..5c096cb18
--- /dev/null
+++ b/benchmarks/benchmark_multipole_pyat.py
@@ -0,0 +1,383 @@
+#!/usr/bin/env python3
+# Import pyAT components
+import at
+import numpy as np
+import pytest
+import torch
+
+# Import Cheetah components
+from cheetah.accelerator.multipole import Multipole
+from cheetah.particles import ParticleBeam, Species
+
+# Configuration parameters
+ELEMENT_LENGTH = 0.5  # meters
+SHORT_LENGTH = 1e-9  # meters for short dipole test
+NUM_STEPS = 10  # Number of integration steps
+
+# Magnet strengths
+DRIFT_K1 = 0.0  # 1/meters^2, no focusing (drift)
+QUAD_K1 = 2.0  # 1/meters^2, quadrupole strength
+SEXT_K2 = 10.0  # 1/meters^3, sextupole strength
+DIPOLE_K0 = 1.0  # 1/meters, dipole strength
+
+# Fringe field flags
+FRINGE_ON = 1  # Enable fringe fields
+FRINGE_OFF = 0  # Disable fringe fields
+
+# Radiation flags
+RADIATION_ON = True  # Enable radiation effects
+RADIATION_OFF = False  # Disable radiation effects
+
+# Initial particle coordinates
+X0 = 0.001  # meters
+PX0 = 0.0  # radians
+Y0 = 0.0  # meters
+PY0 = 0.001  # radians
+TAU0 = 0.1  # seconds
+DELTA0 = -0.1  # relative energy deviation
+
+# Tolerance for coordinate comparisons
+TOLERANCE = 5e-8  # Maximum acceptable difference between PyAT and Cheetah results
+
+
+def create_pyat_quadrupole(length, k1, num_steps, fringe=0, radiation=False):
+    """Create a quadrupole element in pyAT"""
+    pass_method = (
+        "StrMPoleSymplectic4RadPass" if radiation else "StrMPoleSymplectic4Pass"
+    )
+    quad = at.Quadrupole(
+        "Q1",
+        length,
+        k1,
+        NumIntSteps=num_steps,
+        Energy=6.0e9,
+        FringeQuadEntrance=fringe,
+        FringeQuadExit=fringe,
+        PassMethod=pass_method,
+    )
+    return quad
+
+
+def create_pyat_sextupole(length, k2, num_steps, radiation=False):
+    """Create a sextupole element in pyAT"""
+    pass_method = (
+        "StrMPoleSymplectic4RadPass" if radiation else "StrMPoleSymplectic4Pass"
+    )
+    sext = at.Sextupole(
+        "S1", length, k2, NumIntSteps=num_steps, Energy=6.0e9, PassMethod=pass_method
+    )
+    return sext
+
+
+def create_pyat_corrector(kick_x, kick_y, length):
+    """Create a corrector element in pyAT"""
+    corrector = at.Corrector("C1", length, np.array([kick_x, kick_y]))
+    return corrector
+
+
+def create_pyat_dipole(length, k0, num_steps, radiation=False):
+    """Create a dipole element in pyAT"""
+    pass_method = (
+        "StrMPoleSymplectic4RadPass" if radiation else "StrMPoleSymplectic4Pass"
+    )
+
+    # For a dipole, use PolynomB[0] = k0
+    poly_a = []  # No skew components
+    poly_b = [k0]  # Dipole component only
+
+    dipole = at.Multipole(
+        "D1",
+        length,
+        poly_a,
+        poly_b,
+        MaxOrder=0,  # Only using dipole component
+        NumIntSteps=num_steps,
+        Energy=6.0e9,
+        PassMethod=pass_method,
+    )
+    return dipole
+
+
+def create_cheetah_dipole(length, k0, num_steps, radiation=False):
+    """Create a dipole element in Cheetah"""
+    polynom_b = torch.zeros(1, dtype=torch.float64)
+    polynom_b[0] = k0
+
+    tracking_method = "symplectic4_rad" if radiation else "symplectic4"
+
+    dipole = Multipole(
+        length=torch.tensor(length, dtype=torch.float64),
+        polynom_b=polynom_b,
+        max_order=0,
+        num_steps=num_steps,
+        tracking_method=tracking_method,
+        name="D1",
+    )
+    return dipole
+
+
+def create_cheetah_quadrupole(length, k1, num_steps, fringe=0, radiation=False):
+    """Create a quadrupole element in Cheetah"""
+    polynom_b = torch.zeros(2, dtype=torch.float64)
+    polynom_b[1] = k1
+
+    tracking_method = "symplectic4_rad" if radiation else "symplectic4"
+
+    quad = Multipole(
+        length=torch.tensor(length, dtype=torch.float64),
+        polynom_b=polynom_b,
+        max_order=1,
+        num_steps=num_steps,
+        fringe_quad_entrance=fringe,
+        fringe_quad_exit=fringe,
+        tracking_method=tracking_method,
+        name="Q1",
+    )
+    return quad
+
+
+def create_cheetah_sextupole(length, k2, num_steps, radiation=False):
+    """Create a sextupole element in Cheetah"""
+    polynom_b = torch.zeros(3, dtype=torch.float64)
+    polynom_b[2] = k2
+
+    tracking_method = "symplectic4_rad" if radiation else "symplectic4"
+
+    sext = Multipole(
+        length=torch.tensor(length, dtype=torch.float64),
+        polynom_b=polynom_b,
+        max_order=2,
+        num_steps=num_steps,
+        tracking_method=tracking_method,
+        name="S1",
+    )
+    return sext
+
+
+def track_pyat(element, particle):
+    """Track particle through pyAT element"""
+    # Make a copy to avoid modifying the original
+    particle_out = element.track(particle.copy())
+    return particle_out
+
+
+def track_cheetah(element, particle):
+    """Track particle through Cheetah element"""
+    electron = Species("electron")
+    beam = ParticleBeam(
+        particles=particle,
+        energy=torch.tensor(6.0e9, dtype=torch.float64),
+        species=electron,
+    )
+
+    beam_out = element.track(beam)
+    return beam_out
+
+
+def get_aligned_coordinates(pyat_final):
+    """
+    Align coordinates from pyAT and Cheetah for comparison
+
+    PyAT: [0]=x, [1]=px, [2]=y, [3]=py, [4]=delta, [5]=tau
+    Cheetah: [0]=x, [1]=px, [2]=y, [3]=py, [4]=tau, [5]=delta
+    """
+    aligned_pyat = np.array(
+        [
+            pyat_final[0],  # x
+            pyat_final[1],  # px
+            pyat_final[2],  # y
+            pyat_final[3],  # py
+            pyat_final[5],  # tau (at[5] -> cheetah[4])
+            pyat_final[4],  # delta (at[4] -> cheetah[5])
+        ]
+    )
+    return aligned_pyat
+
+
+def run_element_tracking(
+    element_type, length, strength, num_steps, fringe=0, radiation=False
+):
+    """Run tracking for a specific element type and return results for comparison"""
+    # Create elements based on type
+    if element_type == "drift":
+        pyat_element = create_pyat_quadrupole(
+            length, 0.0, num_steps, fringe, radiation
+        )  # Drift is quad with k1=0
+        cheetah_element = create_cheetah_quadrupole(
+            length, 0.0, num_steps, fringe, radiation
+        )
+    elif element_type == "quadrupole":
+        pyat_element = create_pyat_quadrupole(
+            length, strength, num_steps, fringe, radiation
+        )
+        cheetah_element = create_cheetah_quadrupole(
+            length, strength, num_steps, fringe, radiation
+        )
+    elif element_type == "sextupole":
+        pyat_element = create_pyat_sextupole(length, strength, num_steps, radiation)
+        cheetah_element = create_cheetah_sextupole(
+            length, strength, num_steps, radiation
+        )
+    elif element_type == "dipole":
+        pyat_element = create_pyat_dipole(length, strength, num_steps, radiation)
+        cheetah_element = create_cheetah_dipole(length, strength, num_steps, radiation)
+    elif element_type == "corrector":
+        # For corrector, convert integrated dipole strength to kick angle
+        # kick = B0 * L / (Bρ), where Bρ = p/q ≈ E/c for ultrarelativistic particles
+        # For E = 6 GeV, Bρ ≈ 20 T⋅m
+        brho = 20.0  # T⋅m for 6 GeV electrons
+        kick_angle = strength * length / brho  # Convert B0*L to angle in radians
+
+        pyat_element = create_pyat_corrector(kick_angle, 0.0, length)
+        # Corrector doesn't support radiation, so force it off for Cheetah comparison
+        cheetah_element = create_cheetah_dipole(
+            length, strength, num_steps, radiation=False
+        )
+    else:
+        raise ValueError(f"Unknown element type: {element_type}")
+
+    # Create particles
+    # PyAT coordinates: [0]=x, [1]=px, [2]=y, [3]=py, [4]=delta, [5]=tau
+    pyat_particle = np.array([[X0, PX0, Y0, PY0, DELTA0, TAU0]]).T
+
+    # Cheetah coordinates: [0]=x, [1]=px, [2]=y, [3]=py, [4]=tau, [5]=delta, [6]=1.0
+    # (particle flag)
+    cheetah_particle = torch.tensor(
+        [[X0, PX0, Y0, PY0, TAU0, DELTA0, 1.0]], dtype=torch.float64
+    )
+
+    # Track particles
+    pyat_out = track_pyat(pyat_element, pyat_particle)
+    cheetah_out = track_cheetah(cheetah_element, cheetah_particle)
+
+    # Get final coordinates
+    pyat_final = pyat_out.flatten()
+    cheetah_final = cheetah_out.particles[0, :6].cpu().numpy()
+
+    return pyat_final, cheetah_final
+
+
+# Helper function to assert coordinates are close (for pytest)
+def assert_coordinates_close(pyat_final, cheetah_final, test_name):
+    """Assert that pyAT and Cheetah coordinates are within tolerance"""
+    aligned_pyat = get_aligned_coordinates(pyat_final)
+    max_diff = np.max(np.abs(aligned_pyat - cheetah_final))
+
+    assert (
+        max_diff <= TOLERANCE
+    ), f"Coordinates differ by {max_diff} in {test_name} test"
+
+
+# Converted tests to pytest style
+def test_drift():
+    """Test drift space tracking"""
+    pyat_final, cheetah_final = run_element_tracking(
+        "drift", ELEMENT_LENGTH, DRIFT_K1, NUM_STEPS
+    )
+    assert_coordinates_close(pyat_final, cheetah_final, "drift")
+
+    # Verify that PyAT and Cheetah agree on the drift behavior by comparing their
+    # coordinate changes
+    aligned_pyat = get_aligned_coordinates(pyat_final)
+
+    # For drift, check that x and y changes match between PyAT and Cheetah
+    assert (
+        abs((aligned_pyat[0] - X0) - (cheetah_final[0] - X0)) <= TOLERANCE
+    ), "Drift x coordinate change doesn't match between PyAT and Cheetah"
+    assert (
+        abs((aligned_pyat[2] - Y0) - (cheetah_final[2] - Y0)) <= TOLERANCE
+    ), "Drift y coordinate change doesn't match between PyAT and Cheetah"
+
+
+def test_quadrupole():
+    """Test quadrupole tracking"""
+    pyat_final, cheetah_final = run_element_tracking(
+        "quadrupole", ELEMENT_LENGTH, QUAD_K1, NUM_STEPS
+    )
+    assert_coordinates_close(pyat_final, cheetah_final, "quadrupole")
+
+
+def test_sextupole():
+    """Test sextupole tracking"""
+    pyat_final, cheetah_final = run_element_tracking(
+        "sextupole", ELEMENT_LENGTH, SEXT_K2, NUM_STEPS
+    )
+    assert_coordinates_close(pyat_final, cheetah_final, "sextupole")
+
+
+def test_dipole_long():
+    """Test long dipole tracking"""
+    pyat_final, cheetah_final = run_element_tracking(
+        "dipole", ELEMENT_LENGTH, DIPOLE_K0, NUM_STEPS
+    )
+    assert_coordinates_close(pyat_final, cheetah_final, "dipole_long")
+
+
+def test_dipole_short():
+    """Test short dipole tracking"""
+    pyat_final, cheetah_final = run_element_tracking(
+        "corrector", SHORT_LENGTH, DIPOLE_K0, NUM_STEPS
+    )
+    assert_coordinates_close(pyat_final, cheetah_final, "corrector")
+
+
+def test_dipole_long_radiation():
+    """Test long dipole tracking with radiation"""
+    pyat_final, cheetah_final = run_element_tracking(
+        "dipole", ELEMENT_LENGTH, DIPOLE_K0, NUM_STEPS, radiation=RADIATION_ON
+    )
+    assert_coordinates_close(pyat_final, cheetah_final, "dipole_long_radiation")
+
+
+def test_dipole_short_radiation():
+    """Test short dipole tracking with radiation"""
+    pyat_final, cheetah_final = run_element_tracking(
+        "corrector", SHORT_LENGTH, DIPOLE_K0, NUM_STEPS, radiation=RADIATION_ON
+    )
+    assert_coordinates_close(pyat_final, cheetah_final, "corrector_radiation")
+
+
+def test_quadrupole_fringe():
+    """Test quadrupole with fringe fields"""
+    pyat_final, cheetah_final = run_element_tracking(
+        "quadrupole", ELEMENT_LENGTH, QUAD_K1, NUM_STEPS, fringe=FRINGE_ON
+    )
+    assert_coordinates_close(pyat_final, cheetah_final, "quadrupole_fringe")
+
+
+def test_quadrupole_radiation():
+    """Test quadrupole with radiation"""
+    pyat_final, cheetah_final = run_element_tracking(
+        "quadrupole", ELEMENT_LENGTH, QUAD_K1, NUM_STEPS, radiation=RADIATION_ON
+    )
+    assert_coordinates_close(pyat_final, cheetah_final, "quadrupole_radiation")
+
+
+def test_sextupole_radiation():
+    """Test sextupole with radiation"""
+    pyat_final, cheetah_final = run_element_tracking(
+        "sextupole", ELEMENT_LENGTH, SEXT_K2, NUM_STEPS, radiation=RADIATION_ON
+    )
+    assert_coordinates_close(pyat_final, cheetah_final, "sextupole_radiation")
+
+
+def test_quadrupole_fringe_radiation():
+    """Test quadrupole with fringe fields and radiation"""
+    pyat_final, cheetah_final = run_element_tracking(
+        "quadrupole",
+        ELEMENT_LENGTH,
+        QUAD_K1,
+        NUM_STEPS,
+        fringe=FRINGE_ON,
+        radiation=RADIATION_ON,
+    )
+    assert_coordinates_close(pyat_final, cheetah_final, "quadrupole_fringe_radiation")
+
+
+if __name__ == "__main__":
+    # Provide backwards compatibility for running with python directly
+    import sys
+
+    import pytest as pytest_main
+
+    sys.exit(pytest_main.main(["-v", __file__]))
diff --git a/cheetah/__init__.py b/cheetah/__init__.py
index 5a37aa789..0de8f0d4e 100644
--- a/cheetah/__init__.py
+++ b/cheetah/__init__.py
@@ -9,6 +9,7 @@
     Element,
     HorizontalCorrector,
     Marker,
+    Multipole,
     Quadrupole,
     RBend,
     Screen,
diff --git a/cheetah/accelerator/__init__.py b/cheetah/accelerator/__init__.py
index bc783dcee..abedfe1ea 100644
--- a/cheetah/accelerator/__init__.py
+++ b/cheetah/accelerator/__init__.py
@@ -7,6 +7,7 @@
 from .element import Element  # noqa: F401
 from .horizontal_corrector import HorizontalCorrector  # noqa: F401
 from .marker import Marker  # noqa: F401
+from .multipole import Multipole  # noqa: F401
 from .quadrupole import Quadrupole  # noqa: F401
 from .rbend import RBend  # noqa: F401
 from .screen import Screen  # noqa: F401
diff --git a/cheetah/accelerator/multipole.py b/cheetah/accelerator/multipole.py
new file mode 100644
index 000000000..a19757813
--- /dev/null
+++ b/cheetah/accelerator/multipole.py
@@ -0,0 +1,883 @@
+from typing import Literal
+
+import matplotlib.pyplot as plt
+import torch
+from matplotlib.patches import Rectangle
+
+from cheetah.accelerator.element import Element
+from cheetah.particles import Beam, ParticleBeam
+from cheetah.utils import UniqueNameGenerator, verify_device_and_dtype
+
+generate_unique_name = UniqueNameGenerator(prefix="unnamed_element")
+
+
+class Multipole(Element):
+    """
+    Multipole magnet in a particle accelerator.
+
+    :param length: Length in meters.
+    :param polynom_a: Coefficients for skew multipole components (A_n) in 1/m^(n+1)
+    :param polynom_b: Coefficients for normal multipole components (B_n) in 1/m^(n+1)
+    :param max_order: Maximum order of the multipole field.
+    :param misalignment: Misalignment vector of the element in x- and y-directions.
+    :param tilt: Tilt angle of the element in x-y plane [rad].
+    :param num_steps: Number of integration steps.
+    :param fringe_quad_entrance: Whether to apply quadrupole fringe field at entrance
+        (0=no fringe, 1=Lee-Whiting, 2=Lee-Whiting+Elegant).
+    :param fringe_quad_exit: Whether to apply quadrupole fringe field at exit
+        (0=no fringe, 1=Lee-Whiting, 2=Lee-Whiting+Elegant).
+    :param fringe_int_m0: Fringe field integrals for entrance
+        (I0m/K1, I1m/K1, I2m/K1, I3m/K1, Lambda2m/K1).
+    :param fringe_int_p0: Fringe field integrals for exit
+        (I0p/K1, I1p/K1, I2p/K1, I3p/K1, Lambda2p/K1).
+    :param tracking_method: Method to use for tracking through the element.
+    :param name: Unique identifier of the element.
+    """
+
+    def __init__(
+        self,
+        length: torch.Tensor,
+        polynom_a: torch.Tensor | None = None,
+        polynom_b: torch.Tensor | None = None,
+        max_order: int = 1,
+        misalignment: torch.Tensor | None = None,
+        tilt: torch.Tensor | None = None,
+        num_steps: int = 1,
+        fringe_quad_entrance: int = 0,
+        fringe_quad_exit: int = 0,
+        fringe_int_m0: torch.Tensor | None = None,
+        fringe_int_p0: torch.Tensor | None = None,
+        tracking_method: Literal["symplectic4", "symplectic4_rad"] = "symplectic4",
+        name: str | None = None,
+        device: torch.device | None = None,
+        dtype: torch.dtype | None = None,
+    ) -> None:
+        device, dtype = verify_device_and_dtype(
+            [
+                length,
+                polynom_a,
+                polynom_b,
+                misalignment,
+                tilt,
+                fringe_int_m0,
+                fringe_int_p0,
+            ],
+            device,
+            dtype,
+        )
+        factory_kwargs = {"device": device, "dtype": dtype}
+        super().__init__(name=name, **factory_kwargs)
+
+        self.length = torch.as_tensor(length, **factory_kwargs)
+
+        # Set up multipole coefficients
+        if polynom_a is None:
+            # Initialize to zeros with size max_order+1
+            polynom_a = torch.zeros(max_order + 1, **factory_kwargs)
+        else:
+            polynom_a = torch.as_tensor(polynom_a, **factory_kwargs)
+            if polynom_a.size(0) <= max_order:
+                # Pad with zeros if needed
+                padding = torch.zeros(
+                    max_order + 1 - polynom_a.size(0), **factory_kwargs
+                )
+                polynom_a = torch.cat([polynom_a, padding])
+            else:
+                polynom_a = polynom_a[: max_order + 1]
+
+        if polynom_b is None:
+            # Initialize to zeros with size max_order+1
+            polynom_b = torch.zeros(max_order + 1, **factory_kwargs)
+        else:
+            polynom_b = torch.as_tensor(polynom_b, **factory_kwargs)
+            if polynom_b.size(0) <= max_order:
+                # Pad with zeros if needed
+                padding = torch.zeros(
+                    max_order + 1 - polynom_b.size(0), **factory_kwargs
+                )
+                polynom_b = torch.cat([polynom_b, padding])
+            else:
+                polynom_b = polynom_b[: max_order + 1]
+
+        self.register_buffer_or_parameter("polynom_a", polynom_a)
+        self.register_buffer_or_parameter("polynom_b", polynom_b)
+
+        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)
+        )
+
+        # Fringe field parameters
+        self.fringe_quad_entrance = fringe_quad_entrance
+        self.fringe_quad_exit = fringe_quad_exit
+
+        # Default fringe integrals if not provided
+        default_fringe_ints = torch.tensor([0.0, 0.5, 0.0, 0.0, 0.0], **factory_kwargs)
+
+        if fringe_int_m0 is not None:
+            self.register_buffer_or_parameter(
+                "fringe_int_m0", torch.as_tensor(fringe_int_m0, **factory_kwargs)
+            )
+        else:
+            self.register_buffer_or_parameter("fringe_int_m0", default_fringe_ints)
+
+        if fringe_int_p0 is not None:
+            self.register_buffer_or_parameter(
+                "fringe_int_p0", torch.as_tensor(fringe_int_p0, **factory_kwargs)
+            )
+        else:
+            self.register_buffer_or_parameter("fringe_int_p0", default_fringe_ints)
+
+        self.max_order = max_order
+        self.num_steps = num_steps
+        self.tracking_method = tracking_method
+
+    def track(self, incoming: Beam) -> Beam:
+        """
+        Track particles through the multipole element.
+
+        :param incoming: Beam entering the element.
+        :return: Beam exiting the element.
+        """
+        assert isinstance(
+            incoming, ParticleBeam
+        ), "Tracking is only supported for `ParticleBeam`."
+
+        if self.tracking_method == "symplectic4":
+            return self._track_symplectic4(incoming)
+        elif self.tracking_method == "symplectic4_rad":
+            return self._track_symplectic4_rad(incoming)
+        else:
+            raise ValueError(
+                f"Invalid tracking method {self.tracking_method}. "
+                + "Supported methods are 'symplectic4' and 'symplectic4_rad'."
+            )
+
+    def _track_symplectic4(self, incoming: ParticleBeam) -> ParticleBeam:
+        """
+        Track particles through the multipole element using the 4th order symplectic
+        integration.
+
+        :param incoming: ParticleBeam entering the element.
+        :return: ParticleBeam exiting the element.
+        """
+        # Constants for 4th order symplectic integration
+        DRIFT1 = 0.6756035959798286638
+        DRIFT2 = -0.1756035959798286639
+        KICK1 = 1.351207191959657328
+        KICK2 = -1.702414383919314656
+
+        # Get particle coordinates
+        x = incoming.x
+        px = incoming.px
+        y = incoming.y
+        py = incoming.py
+        tau = incoming.tau
+        delta = incoming.p
+
+        # Apply misalignment at entrance if needed
+        if torch.any(self.misalignment != 0) or torch.any(self.tilt != 0):
+            # For now, let's ignore tilt and just apply simple offset
+            x_offset = self.misalignment[..., 0]
+            y_offset = self.misalignment[..., 1]
+            x = x - x_offset
+            y = y - y_offset
+
+        # Check if quadrupole component exists (b2 != 0)
+        b2 = 0.0
+        if self.polynom_b.size(0) > 1:
+            b2 = self.polynom_b[1].item()
+
+        # Apply fringe field at entrance if enabled and we have a quadrupole component
+        if self.fringe_quad_entrance and b2 != 0:
+            if self.fringe_quad_entrance == 1:
+                # Apply simple Lee-Whiting fringe field
+                self._quad_fringe_pass_p(x, px, y, py, tau, delta, b2)
+            elif self.fringe_quad_entrance == 2:
+                # Apply Elegant-style fringe field
+                self._apply_linear_quad_fringe_entrance(x, px, y, py, tau, delta, b2)
+
+        # Prepare for tracking
+        SL = self.length / self.num_steps
+        L1 = SL * DRIFT1
+        L2 = SL * DRIFT2
+        K1 = SL * KICK1
+        K2 = SL * KICK2
+
+        # Symplectic integration loop
+        for _ in range(self.num_steps):
+            # First drift
+            norm = 1.0 / (1.0 + delta)
+            norm_l1 = L1 * norm
+            x += norm_l1 * px
+            y += norm_l1 * py
+            tau += norm_l1 * (px * px + py * py) / (2.0 * (1.0 + delta))
+
+            # First kick
+            self._apply_kick(x, y, px, py, K1)
+
+            # Second drift
+            norm_l2 = L2 * norm
+            x += norm_l2 * px
+            y += norm_l2 * py
+            tau += norm_l2 * (px * px + py * py) / (2.0 * (1.0 + delta))
+
+            # Second kick
+            self._apply_kick(x, y, px, py, K2)
+
+            # Third drift
+            x += norm_l2 * px
+            y += norm_l2 * py
+            tau += norm_l2 * (px * px + py * py) / (2.0 * (1.0 + delta))
+
+            # Third kick
+            self._apply_kick(x, y, px, py, K1)
+
+            # Fourth drift
+            x += norm_l1 * px
+            y += norm_l1 * py
+            tau += norm_l1 * (px * px + py * py) / (2.0 * (1.0 + delta))
+
+        # Apply fringe field at exit if enabled and we have a quadrupole component
+        if self.fringe_quad_exit and b2 != 0:
+            if self.fringe_quad_exit == 1:
+                # Apply simple Lee-Whiting fringe field
+                self._quad_fringe_pass_n(x, px, y, py, tau, delta, b2)
+            elif self.fringe_quad_exit == 2:
+                # Apply Elegant-style fringe field
+                self._apply_linear_quad_fringe_exit(x, px, y, py, tau, delta, b2)
+
+        # Apply misalignment at exit if needed
+        if torch.any(self.misalignment != 0) or torch.any(self.tilt != 0):
+            # Restore coordinates to laboratory frame
+            x = x + x_offset
+            y = y + y_offset
+
+        # Create output beam
+        outgoing_beam = ParticleBeam(
+            particles=torch.stack(
+                (x, px, y, py, tau, delta, torch.ones_like(x)), dim=-1
+            ),
+            energy=incoming.energy,
+            particle_charges=incoming.particle_charges,
+            survival_probabilities=incoming.survival_probabilities,
+            species=incoming.species,
+        )
+
+        return outgoing_beam
+
+    def _track_symplectic4_rad(self, incoming: ParticleBeam) -> ParticleBeam:
+        """
+        Track particles through the multipole element using the 4th order symplectic
+        integration with radiation effects included.
+
+        :param incoming: ParticleBeam entering the element.
+        :return: ParticleBeam exiting the element.
+        """
+        # Constants for 4th order symplectic integration
+        DRIFT1 = 0.6756035959798286638
+        DRIFT2 = -0.1756035959798286639
+        KICK1 = 1.351207191959657328
+        KICK2 = -1.702414383919314656
+
+        # Get particle coordinates
+        x = incoming.x
+        px = incoming.px
+        y = incoming.y
+        py = incoming.py
+        tau = incoming.tau
+        delta = incoming.p
+
+        # Apply misalignment at entrance if needed
+        if torch.any(self.misalignment != 0) or torch.any(self.tilt != 0):
+            # For now, let's ignore tilt and just apply simple offset
+            x_offset = self.misalignment[..., 0]
+            y_offset = self.misalignment[..., 1]
+            x = x - x_offset
+            y = y - y_offset
+
+        # Check if quadrupole component exists (b2 != 0)
+        b2 = 0.0
+        if self.polynom_b.size(0) > 1:
+            b2 = self.polynom_b[1].item()
+
+        # Apply fringe field at entrance if enabled and we have a quadrupole component
+        if self.fringe_quad_entrance and b2 != 0:
+            if self.fringe_quad_entrance == 1:
+                # Apply simple Lee-Whiting fringe field
+                self._quad_fringe_pass_p(x, px, y, py, tau, delta, b2)
+            elif self.fringe_quad_entrance == 2:
+                # Apply Elegant-style fringe field
+                self._apply_linear_quad_fringe_entrance(x, px, y, py, tau, delta, b2)
+
+        # Prepare for tracking
+        SL = self.length / self.num_steps
+        L1 = SL * DRIFT1
+        L2 = SL * DRIFT2
+        K1 = SL * KICK1
+        K2 = SL * KICK2
+
+        # Get the beam energy (needed for radiation effects)
+        energy = incoming.energy
+
+        # Symplectic integration loop
+        for _ in range(self.num_steps):
+            # First drift
+            self._drift6(x, px, y, py, tau, delta, L1)
+
+            # First kick with radiation
+            self._apply_kick_rad(x, y, px, py, delta, K1, energy)
+
+            # Second drift
+            self._drift6(x, px, y, py, tau, delta, L2)
+
+            # Second kick with radiation
+            self._apply_kick_rad(x, y, px, py, delta, K2, energy)
+
+            # Third drift
+            self._drift6(x, px, y, py, tau, delta, L2)
+
+            # Third kick with radiation
+            self._apply_kick_rad(x, y, px, py, delta, K1, energy)
+
+            # Fourth drift
+            self._drift6(x, px, y, py, tau, delta, L1)
+
+        # Apply fringe field at exit if enabled and we have a quadrupole component
+        if self.fringe_quad_exit and b2 != 0:
+            if self.fringe_quad_exit == 1:
+                # Apply simple Lee-Whiting fringe field
+                self._quad_fringe_pass_n(x, px, y, py, tau, delta, b2)
+            elif self.fringe_quad_exit == 2:
+                # Apply Elegant-style fringe field
+                self._apply_linear_quad_fringe_exit(x, px, y, py, tau, delta, b2)
+
+        # Apply misalignment at exit if needed
+        if torch.any(self.misalignment != 0) or torch.any(self.tilt != 0):
+            # Restore coordinates to laboratory frame
+            x = x + x_offset
+            y = y + y_offset
+
+        # Create output beam
+        outgoing_beam = ParticleBeam(
+            particles=torch.stack(
+                (x, px, y, py, tau, delta, torch.ones_like(x)), dim=-1
+            ),
+            energy=incoming.energy,
+            particle_charges=incoming.particle_charges,
+            survival_probabilities=incoming.survival_probabilities,
+            species=incoming.species,
+        )
+
+        return outgoing_beam
+
+    def _drift6(
+        self,
+        x: torch.Tensor,
+        px: torch.Tensor,
+        y: torch.Tensor,
+        py: torch.Tensor,
+        tau: torch.Tensor,
+        delta: torch.Tensor,
+        length: float,
+    ) -> None:
+        """
+        Apply a drift to the 6D phase space coordinates.
+
+        This method implements a drift space transformation in full 6D phase space.
+        The longitudinal coordinate (tau) update accounts for path length differences
+        due to transverse motion, which is essential for maintaining symplecticity
+        in the tracking. The implementation follows the approach used in PyAT.
+
+        :param x: Horizontal position
+        :param px: Horizontal momentum
+        :param y: Vertical position
+        :param py: Vertical momentum
+        :param tau: Longitudinal position
+        :param delta: Relative momentum deviation
+        :param length: Length of the drift
+        """
+        # Equivalent to ATdrift6/fastdrift in the C code
+        # Apply drift to positions
+        norm = 1.0 / (1.0 + delta)
+        norm_length = length * norm
+
+        x += norm_length * px
+        y += norm_length * py
+
+        # Update longitudinal coordinate - matching PyAT's implementation
+        # In PyAT, r[5] (tau) += NormL*(r[1]*r[1]+r[3]*r[3])/(2*(1+r[4])) where r[4] is
+        # delta
+        tau += norm_length * (px * px + py * py) / (2.0 * (1.0 + delta))
+
+    def _apply_kick_rad(
+        self,
+        x: torch.Tensor,
+        y: torch.Tensor,
+        px: torch.Tensor,
+        py: torch.Tensor,
+        delta: torch.Tensor,
+        kick_strength: float,
+        energy: torch.Tensor,
+    ) -> None:
+        """
+        Apply multipole kick to particles with radiation effects.
+
+        This method calculates the magnetic field components as in _apply_kick,
+        but also includes synchrotron radiation effects. The radiation calculation
+        follows the model described by M. Sands, which accounts for energy loss
+        due to photon emission in the magnetic field.
+
+        :param x: Horizontal position
+        :param y: Vertical position
+        :param px: Horizontal momentum
+        :param py: Vertical momentum
+        :param delta: Relative momentum deviation
+        :param kick_strength: Strength of the kick (includes step length factor)
+        :param energy: Beam energy in GeV
+        """
+        # BUGFIX: Properly implement PyAT's field calculation algorithm
+        # Always start with the value at max_order index, even if it's zero
+        # This is critical for consistent results regardless of max_order
+        ReSum = self.polynom_b[self.max_order].expand_as(x)
+        ImSum = self.polynom_a[self.max_order].expand_as(y)
+
+        # Apply the recursive algorithm exactly as in PyAT
+        for i in range(self.max_order - 1, -1, -1):
+            ReSumTemp = ReSum * x - ImSum * y + self.polynom_b[i]
+            ImSum = ImSum * x + ReSum * y + self.polynom_a[i]
+            ReSum = ReSumTemp
+
+        # Constants for radiation calculation
+        CGAMMA = 8.846e-5  # Radiation constant for electrons
+        TWOPI = 2 * torch.pi
+
+        # Store original delta for use in kick calculation (as in PyAT)
+        delta_orig = delta.clone()
+
+        # Calculate normalized velocities as in PyAT
+        p_norm = 1.0 / (1.0 + delta)
+        xpr = px * p_norm
+        ypr = py * p_norm
+
+        # Calculate B2P (perpendicular B-field squared) as in PyAT's StrB2perp
+        v_norm2 = 1.0 / (1.0 + xpr**2 + ypr**2)
+        bx = ImSum  # In PyAT, bx = ImSum
+        by = ReSum  # In PyAT, by = ReSum
+        B2P = (by**2 + bx**2 + (bx * ypr - by * xpr) ** 2) * v_norm2
+
+        # Calculate CRAD according to M.Sands (4.1) as in PyAT
+        CRAD = CGAMMA * energy**3 / (TWOPI * 1e27)  # [m]/[GeV^3]
+
+        # For straight elements, irho = 0
+        irho = 0.0
+
+        # Calculate and apply energy loss as in PyAT
+        delta.sub_(
+            CRAD
+            * (1.0 + delta) ** 2
+            * B2P
+            * (1.0 + (xpr**2 + ypr**2) / 2.0)
+            * kick_strength
+        )
+
+        # Recalculate momenta from angles after energy loss
+        p_norm = 1.0 / (1.0 + delta)
+        px.copy_(xpr / p_norm)
+        py.copy_(ypr / p_norm)
+
+        # Apply kicks using the original delta value as in PyAT
+        px.sub_(kick_strength * (ReSum - (delta_orig - x * irho) * irho))
+        py.add_(kick_strength * ImSum)
+
+    def _apply_kick(
+        self,
+        x: torch.Tensor,
+        y: torch.Tensor,
+        px: torch.Tensor,
+        py: torch.Tensor,
+        kick_strength: float,
+    ) -> None:
+        """
+        Apply multipole kick to particles.
+
+        This method calculates the magnetic field components from the multipole
+        expansion and applies the resulting momentum kicks to the particles. The
+        algorithm recursively computes the field components following the same approach
+        as in PyAT.
+
+        :param x: Horizontal position
+        :param px: Horizontal momentum
+        :param y: Vertical position
+        :param py: Vertical momentum
+        :param kick_strength: Strength of the kick (includes step length factor)
+        """
+        # BUGFIX: Properly implement PyAT's field calculation algorithm
+        # Always start with the value at max_order index, even if it's zero
+        # This is critical for consistent results regardless of max_order
+        ReSum = self.polynom_b[self.max_order].expand_as(x)
+        ImSum = self.polynom_a[self.max_order].expand_as(y)
+
+        # Apply the recursive algorithm exactly as in PyAT
+        for i in range(self.max_order - 1, -1, -1):
+            ReSumTemp = ReSum * x - ImSum * y + self.polynom_b[i]
+            ImSum = ImSum * x + ReSum * y + self.polynom_a[i]
+            ReSum = ReSumTemp
+
+        # Apply the kick
+        px -= kick_strength * ReSum
+        py += kick_strength * ImSum
+
+    @property
+    def is_skippable(self) -> bool:
+        """
+        Check if the element can be skipped during tracking.
+        For multipoles, it's skippable if all multipole coefficients are zero.
+        """
+        return (
+            torch.all(self.polynom_a == 0).item()
+            and torch.all(self.polynom_b == 0).item()
+        )
+
+    def split(self, resolution: torch.Tensor) -> list[Element]:
+        """
+        Split the multipole element into smaller pieces for higher-resolution tracking.
+
+        :param resolution: The desired resolution for splitting in meters.
+        :return: List of Multipole elements that together make up the original element.
+        """
+        num_splits = torch.ceil(torch.max(self.length) / resolution).int()
+
+        # For splits, only apply fringe at the very beginning and end
+        split_elements = []
+        for i in range(num_splits):
+            # Only apply fringe at first and last element
+            fringe_entrance = self.fringe_quad_entrance if i == 0 else 0
+            fringe_exit = self.fringe_quad_exit if i == num_splits - 1 else 0
+
+            split_elements.append(
+                Multipole(
+                    self.length / num_splits,
+                    polynom_a=self.polynom_a,
+                    polynom_b=self.polynom_b,
+                    max_order=self.max_order,
+                    misalignment=self.misalignment,
+                    tilt=self.tilt,
+                    num_steps=self.num_steps,
+                    fringe_quad_entrance=fringe_entrance,
+                    fringe_quad_exit=fringe_exit,
+                    fringe_int_m0=self.fringe_int_m0,
+                    fringe_int_p0=self.fringe_int_p0,
+                    tracking_method=self.tracking_method,
+                    dtype=self.length.dtype,
+                    device=self.length.device,
+                )
+            )
+
+        return split_elements
+
+    @property
+    def is_active(self) -> bool:
+        """Check if the element has any non-zero multipole coefficients."""
+        return (torch.any(self.polynom_a != 0) or torch.any(self.polynom_b != 0)).item()
+
+    def plot(self, ax: plt.Axes, s: float, vector_idx: tuple | None = None) -> None:
+        """Plot the multipole element."""
+        plot_s = s[vector_idx] if s.dim() > 0 else s
+        plot_length = self.length[vector_idx] if self.length.dim() > 0 else self.length
+
+        alpha = 1 if self.is_active else 0.2
+        # For general multipole, use a different color
+        patch = Rectangle(
+            (plot_s, 0), plot_length, 0.8, color="tab:purple", alpha=alpha, zorder=2
+        )
+        ax.add_patch(patch)
+
+    @property
+    def defining_features(self) -> list[str]:
+        return super().defining_features + [
+            "length",
+            "polynom_a",
+            "polynom_b",
+            "max_order",
+            "misalignment",
+            "tilt",
+            "fringe_quad_entrance",
+            "fringe_quad_exit",
+            "fringe_int_m0",
+            "fringe_int_p0",
+        ]
+
+    def __repr__(self) -> str:
+        return (
+            f"{self.__class__.__name__}(length={repr(self.length)}, "
+            + f"polynom_a={repr(self.polynom_a)}, "
+            + f"polynom_b={repr(self.polynom_b)}, "
+            + f"max_order={repr(self.max_order)}, "
+            + f"misalignment={repr(self.misalignment)}, "
+            + f"tilt={repr(self.tilt)}, "
+            + f"num_steps={repr(self.num_steps)}, "
+            + f"fringe_quad_entrance={repr(self.fringe_quad_entrance)}, "
+            + f"fringe_quad_exit={repr(self.fringe_quad_exit)}, "
+            + f"tracking_method={repr(self.tracking_method)}, "
+            + f"name={repr(self.name)})"
+        )
+
+    def _quad_fringe_pass_p(
+        self,
+        x: torch.Tensor,
+        px: torch.Tensor,
+        y: torch.Tensor,
+        py: torch.Tensor,
+        tau: torch.Tensor,
+        delta: torch.Tensor,
+        b2: float,
+    ) -> None:
+        """
+        Apply quadrupole fringe field effect at entrance using Lee-Whiting's formula.
+
+        This is a PyTorch implementation of QuadFringePassP from quadfringe.c
+
+        :param x: Horizontal position
+        :param px: Horizontal momentum
+        :param y: Vertical position
+        :param py: Vertical momentum
+        :param tau: Longitudinal position
+        :param delta: Relative momentum deviation
+        :param b2: Quadrupole strength coefficient (polynom_b[1])
+        """
+        # Calculate fringe field parameters
+        u = b2 / (12.0 * (1.0 + delta))
+        x2 = x * x
+        y2 = y * y
+        xy = x * y
+
+        # Calculate displacements
+        gx = u * (x2 + 3 * y2) * x
+        gy = u * (y2 + 3 * x2) * y
+
+        # Apply position changes
+        x.add_(gx)
+        y.sub_(gy)
+
+        # Calculate momentum changes
+        px_tmp = 3 * u * (2 * xy * py - (x2 + y2) * px)
+        py_tmp = 3 * u * (2 * xy * px - (x2 + y2) * py)
+
+        # Update longitudinal coordinate
+        tau.sub_((gy * py - gx * px) / (1.0 + delta))
+
+        # Apply momentum changes
+        px.add_(px_tmp)
+        py.sub_(py_tmp)
+
+    def _quad_fringe_pass_n(
+        self,
+        x: torch.Tensor,
+        px: torch.Tensor,
+        y: torch.Tensor,
+        py: torch.Tensor,
+        tau: torch.Tensor,
+        delta: torch.Tensor,
+        b2: float,
+    ) -> None:
+        """
+        Apply quadrupole fringe field effect at exit using Lee-Whiting's formula.
+
+        This is a PyTorch implementation of QuadFringePassN from quadfringe.c
+
+        :param x: Horizontal position
+        :param px: Horizontal momentum
+        :param y: Vertical position
+        :param py: Vertical momentum
+        :param tau: Longitudinal position
+        :param delta: Relative momentum deviation
+        :param b2: Quadrupole strength coefficient (polynom_b[1])
+        """
+        # Calculate fringe field parameters
+        u = b2 / (12.0 * (1.0 + delta))
+        x2 = x * x
+        y2 = y * y
+        xy = x * y
+
+        # Calculate displacements
+        gx = u * (x2 + 3 * y2) * x
+        gy = u * (y2 + 3 * x2) * y
+
+        # Apply position changes (opposite signs from entrance)
+        x.sub_(gx)
+        y.add_(gy)
+
+        # Calculate momentum changes
+        px_tmp = 3 * u * (2 * xy * py - (x2 + y2) * px)
+        py_tmp = 3 * u * (2 * xy * px - (x2 + y2) * py)
+
+        # Update longitudinal coordinate
+        tau.add_((gy * py - gx * px) / (1.0 + delta))
+
+        # Apply momentum changes (opposite signs from entrance)
+        px.sub_(px_tmp)
+        py.add_(py_tmp)
+
+    def _quad_partial_fringe_matrix(
+        self, K1: torch.Tensor, in_fringe: float, fringe_int: torch.Tensor, part: int
+    ) -> torch.Tensor:
+        """
+        Generate partial fringe matrix for quadrupole fringe effects (Elegant-style).
+
+        :param K1: Quadrupole strength divided by (1+delta)
+        :param in_fringe: Fringe direction (-1 for entrance, 1 for exit)
+        :param fringe_int: Fringe integrals
+        :param part: Part number (1 or 2)
+        :return: 6x6 transfer matrix
+        """
+        # Initialize 6x6 identity matrix
+        R = torch.eye(6, device=K1.device, dtype=K1.dtype).expand(K1.shape[0], 6, 6)
+
+        # Square of K1
+        K1sqr = K1 * K1
+
+        # Calculate J parameters according to part number
+        if part == 1:
+            J1x = in_fringe * (K1 * fringe_int[1] - 2 * K1sqr * fringe_int[3] / 3.0)
+            J2x = in_fringe * (K1 * fringe_int[2])
+            J3x = in_fringe * (K1sqr * (fringe_int[2] + fringe_int[4]))
+
+            # For y-plane, use negative K1
+            K1_y = -K1
+            J1y = in_fringe * (
+                K1_y * fringe_int[1] - 2 * K1_y * K1_y * fringe_int[3] / 3.0
+            )
+            J2y = -J2x
+            J3y = J3x
+        else:  # part == 2
+            J1x = in_fringe * (
+                K1 * fringe_int[1] + K1sqr * fringe_int[0] * fringe_int[2] / 2
+            )
+            J2x = in_fringe * (K1 * fringe_int[2])
+            J3x = in_fringe * (K1sqr * (fringe_int[4] - fringe_int[0] * fringe_int[1]))
+
+            # For y-plane, use negative K1
+            K1_y = -K1
+            J1y = in_fringe * (
+                K1_y * fringe_int[1] + K1_y * K1_y * fringe_int[0] * fringe_int[2]
+            )
+            J2y = -J2x
+            J3y = J3x
+
+        # Calculate matrix elements
+        exp_J1x = torch.exp(J1x)
+        R[:, 0, 0] = exp_J1x
+        R[:, 0, 1] = J2x / exp_J1x
+        R[:, 1, 0] = exp_J1x * J3x
+        R[:, 1, 1] = (1.0 + J2x * J3x) / exp_J1x
+
+        exp_J1y = torch.exp(J1y)
+        R[:, 2, 2] = exp_J1y
+        R[:, 2, 3] = J2y / exp_J1y
+        R[:, 3, 2] = exp_J1y * J3y
+        R[:, 3, 3] = (1.0 + J2y * J3y) / exp_J1y
+
+        return R
+
+    def _apply_linear_quad_fringe_entrance(
+        self,
+        x: torch.Tensor,
+        px: torch.Tensor,
+        y: torch.Tensor,
+        py: torch.Tensor,
+        tau: torch.Tensor,
+        delta: torch.Tensor,
+        b2: float,
+    ) -> None:
+        """
+        Apply quadrupole fringe field at entrance including linear Elegant-style
+        effects.
+
+        :param x: Horizontal position
+        :param px: Horizontal momentum
+        :param y: Vertical position
+        :param py: Vertical momentum
+        :param tau: Longitudinal position
+        :param delta: Relative momentum deviation
+        :param b2: Quadrupole strength coefficient (polynom_b[1])
+        """
+        # Elegant-style fringe field entrance (linear part)
+        in_fringe = -1.0
+        k1 = b2 / (1.0 + delta)
+
+        # First linear matrix
+        R = self._quad_partial_fringe_matrix(k1, in_fringe, self.fringe_int_p0, 1)
+
+        # Apply first matrix
+        x_new = R[:, 0, 0] * x + R[:, 0, 1] * px
+        px_new = R[:, 1, 0] * x + R[:, 1, 1] * px
+        y_new = R[:, 2, 2] * y + R[:, 2, 3] * py
+        py_new = R[:, 3, 2] * y + R[:, 3, 3] * py
+
+        x.copy_(x_new)
+        px.copy_(px_new)
+        y.copy_(y_new)
+        py.copy_(py_new)
+
+        # Apply nonlinear fringe field (AT code)
+        self._quad_fringe_pass_p(x, px, y, py, tau, delta, b2)
+
+        # Second linear matrix
+        R = self._quad_partial_fringe_matrix(k1, in_fringe, self.fringe_int_m0, 2)
+
+        # Apply second matrix
+        x_new = R[:, 0, 0] * x + R[:, 0, 1] * px
+        px_new = R[:, 1, 0] * x + R[:, 1, 1] * px
+        y_new = R[:, 2, 2] * y + R[:, 2, 3] * py
+        py_new = R[:, 3, 2] * y + R[:, 3, 3] * py
+
+        x.copy_(x_new)
+        px.copy_(px_new)
+        y.copy_(y_new)
+        py.copy_(py_new)
+
+    def _apply_linear_quad_fringe_exit(
+        self,
+        x: torch.Tensor,
+        px: torch.Tensor,
+        y: torch.Tensor,
+        py: torch.Tensor,
+        tau: torch.Tensor,
+        delta: torch.Tensor,
+        b2: float,
+    ) -> None:
+        """
+        Apply quadrupole fringe field at exit including linear Elegant-style effects.
+
+        :param x: Horizontal position
+        :param px: Horizontal momentum
+        :param y: Vertical position
+        :param py: Vertical momentum
+        :param tau: Longitudinal position
+        :param delta: Relative momentum deviation
+        :param b2: Quadrupole strength coefficient (polynom_b[1])
+        """
+        # Elegant-style fringe field exit (linear part)
+        in_fringe = 1.0
+        k1 = b2 / (1.0 + delta)
+
+        # First linear matrix
+        R = self._quad_partial_fringe_matrix(k1, in_fringe, self.fringe_int_m0, 1)
+
+        # Apply first matrix
+        x_new = R[:, 0, 0] * x + R[:, 0, 1] * px
+        px_new = R[:, 1, 0] * x + R[:, 1, 1] * px
+        y_new = R[:, 2, 2] * y + R[:, 2, 3] * py
+        py_new = R[:, 3, 2] * y + R[:, 3, 3] * py
+
+        x.copy_(x_new)
+        px.copy_(px_new)
+        y.copy_(y_new)
+        py.copy_(py_new)
diff --git a/tests/test_multipole.py b/tests/test_multipole.py
new file mode 100644
index 000000000..42604edff
--- /dev/null
+++ b/tests/test_multipole.py
@@ -0,0 +1,332 @@
+import pytest
+import torch
+
+import cheetah
+
+
+def test_multipole_as_drift():
+    """Test that a multipole with all coefficients zero approximates drift behaviour."""
+    multipole = cheetah.Multipole(length=torch.tensor(1.0))
+    drift = cheetah.Drift(length=torch.tensor(1.0))
+
+    incoming = cheetah.ParticleBeam.from_parameters(
+        num_particles=1_000, energy=torch.tensor(1e9)
+    )
+
+    outgoing_multipole = multipole.track(incoming)
+    outgoing_drift = drift.track(incoming)
+
+    assert torch.allclose(
+        outgoing_multipole.particles, outgoing_drift.particles, rtol=1e-6
+    )
+
+
+def test_multipole_as_quadrupole():
+    """
+    Compare a multipole configured as a quadrupole with a native quadrupole element.
+    """
+    k1 = torch.tensor(4.2)
+
+    multipole = cheetah.Multipole(
+        length=torch.tensor(1.0), polynom_b=torch.tensor([0.0, k1])  # B1 = k1
+    )
+    quadrupole = cheetah.Quadrupole(
+        length=torch.tensor(1.0),
+        k1=k1,
+        tracking_method="bmadx",
+        num_steps=10,  # Use multiple steps for better accuracy
+    )
+
+    incoming = cheetah.ParticleBeam.from_parameters(
+        num_particles=1_000, energy=torch.tensor(1e9), mu_x=torch.tensor(1e-3)
+    )
+
+    outgoing_multipole = multipole.track(incoming)
+    outgoing_quadrupole = quadrupole.track(incoming)
+
+    assert torch.allclose(
+        outgoing_multipole.particles, outgoing_quadrupole.particles, rtol=1e-6
+    )
+
+
+def test_multipole_as_horizontal_corrector():
+    """
+    Test that a multipole with B0 coefficient behaves like a horizontal corrector
+    element.
+    """
+    length = torch.tensor(1e-9)
+    angle = torch.tensor(2e-3)
+
+    corrector = cheetah.HorizontalCorrector(length=length, angle=angle)
+    multipole = cheetah.Multipole(
+        length=length,
+        polynom_b=torch.tensor([-angle / length, 0.0]),  # B0 = -angle/length
+        max_order=1,
+    )
+
+    beam = cheetah.ParticleBeam.from_parameters(
+        num_particles=1_000, energy=torch.tensor(1e9)
+    )
+
+    outgoing_corrector = corrector.track(beam)
+    outgoing_multipole = multipole.track(beam)
+
+    assert torch.allclose(
+        outgoing_corrector.particles, outgoing_multipole.particles, rtol=1e-6
+    )
+
+
+def test_multipole_with_misalignment():
+    """
+    Test that a multipole with misalignment behaves equivalently to a centered multipole
+    with an offset beam (with appropriate transformation).
+    """
+    length = torch.tensor(1.0)
+    k1 = torch.tensor(20.0)
+    offset_x = torch.tensor(0.001)
+    offset_y = torch.tensor(0.002)
+
+    # Create a quadrupole-like multipole with and without misalignment
+    multipole_misaligned = cheetah.Multipole(
+        length=length,
+        polynom_b=torch.tensor([0.0, k1]),
+        misalignment=torch.tensor([offset_x, offset_y]),
+    )
+    multipole_centered = cheetah.Multipole(
+        length=length,
+        polynom_b=torch.tensor([0.0, k1]),
+    )
+
+    # Create a beam centered at origin
+    centered_beam = cheetah.ParticleBeam.from_parameters(
+        num_particles=1_000, energy=torch.tensor(1e9)
+    )
+
+    # Create an offset beam with opposite offset relative to misalignment
+    offset_beam = cheetah.ParticleBeam.from_parameters(
+        num_particles=1_000, energy=torch.tensor(1e9), mu_x=-offset_x, mu_y=-offset_y
+    )
+
+    # Track through both scenarios
+    outbeam_misaligned_quad = multipole_misaligned.track(centered_beam)
+    outbeam_offset_beam = multipole_centered.track(offset_beam)
+
+    # Transform the offset beam results back to the lab frame
+    # For offset beam case, the beam was offset by (-offset_x, -offset_y)
+    # So we need to translate the results back by (offset_x, offset_y)
+    transformed_offset_beam_mu_x = outbeam_offset_beam.mu_x + offset_x
+    transformed_offset_beam_mu_y = outbeam_offset_beam.mu_y + offset_y
+
+    # Compare the results - they should be nearly identical after transformation
+    assert torch.allclose(
+        outbeam_misaligned_quad.mu_x,
+        transformed_offset_beam_mu_x,
+        atol=1e-5,  # Increased tolerance to account for numerical differences
+    )
+    assert torch.allclose(
+        outbeam_misaligned_quad.mu_y,
+        transformed_offset_beam_mu_y,
+        atol=1e-5,  # Increased tolerance to account for numerical differences
+    )
+    assert torch.allclose(
+        outbeam_misaligned_quad.mu_px,
+        outbeam_offset_beam.mu_px,
+        atol=1e-5,  # Increased tolerance to account for numerical differences
+    )
+    assert torch.allclose(
+        outbeam_misaligned_quad.mu_py,
+        outbeam_offset_beam.mu_py,
+        atol=1e-5,  # Increased tolerance to account for numerical differences
+    )
+
+
+def test_skew_multipole_behavior():
+    """
+    Test that skew multipole components (polynom_a) work as expected (like 90deg tilted
+    normal quadrupole).
+    """
+    # Normal quadrupole (using B1)
+    normal_quad = cheetah.Multipole(
+        length=torch.tensor(1.0),
+        polynom_b=torch.tensor([0.0, 5.0]),  # Strong for clear effect
+    )
+
+    # Skew quadrupole (using A1)
+    skew_quad = cheetah.Multipole(
+        length=torch.tensor(1.0),
+        polynom_a=torch.tensor([0.0, 5.0]),
+    )
+
+    # Create a beam with x offset only
+    incoming = cheetah.ParticleBeam.from_parameters(
+        num_particles=1_000,
+        energy=torch.tensor(1e9),
+        mu_x=torch.tensor(1e-3),
+        mu_y=torch.tensor(0.0),
+    )
+
+    # Track through both elements
+    out_normal = normal_quad.track(incoming)
+    out_skew = skew_quad.track(incoming)
+
+    # Normal and skew quadrupoles should have different coupling behaviors
+    # Skew quadrupole should have larger py change from x offset compared to normal quad
+    assert not torch.allclose(out_normal.mu_py, out_skew.mu_py, atol=1e-10)
+
+
+def test_multipole_max_order():
+    """
+    Test that the max_order parameter correctly limits the polynomial orders.
+    """
+    # Create a multipole with coefficients up to order 5
+    multipole_full = cheetah.Multipole(
+        length=torch.tensor(1.0),
+        polynom_b=torch.tensor([0.0, 1.0, 0.5, 0.2, 0.1, 0.05]),  # Orders 0-5
+        max_order=5,
+    )
+
+    # Create a multipole with the same coefficients but limited to order 2
+    multipole_limited = cheetah.Multipole(
+        length=torch.tensor(1.0),
+        polynom_b=torch.tensor(
+            [0.0, 1.0, 0.5, 0.2, 0.1, 0.05]
+        ),  # Only first 3 will be used
+        max_order=2,
+    )
+
+    # Verify that the coefficients were properly truncated
+    assert multipole_limited.polynom_b.size(0) == 3  # Should contain orders 0, 1, 2
+    assert torch.all(multipole_limited.polynom_b == torch.tensor([0.0, 1.0, 0.5]))
+
+    # Create another multipole with extra coefficients beyond max_order
+    polynom_b_extra = torch.zeros(10)
+    polynom_b_extra[1] = 1.0  # Set B1 = 1.0
+
+    multipole_extra = cheetah.Multipole(
+        length=torch.tensor(1.0),
+        polynom_b=polynom_b_extra,
+        max_order=2,
+    )
+
+    # Should be truncated to just 3 elements
+    assert multipole_extra.polynom_b.size(0) == 3
+
+
+@pytest.mark.parametrize("tracking_method", ["symplectic4", "symplectic4_rad"])
+def test_multipole_tracking_methods(tracking_method):
+    """
+    Test that different tracking methods don't crash.
+    """
+    # Create a multipole with both normal and skew components
+    multipole = cheetah.Multipole(
+        length=torch.tensor(1.0),
+        polynom_b=torch.tensor([0.0, 1.0, 1.0]),  # B1 and B2 components
+        polynom_a=torch.tensor([0.0, 0.5]),  # A1 component
+        tracking_method=tracking_method,
+        num_steps=10,  # Use multiple steps for better accuracy
+    )
+
+    # Create a beam with offset to see effects
+    incoming = cheetah.ParticleBeam.from_parameters(
+        num_particles=1_000, energy=torch.tensor(1e9), mu_x=torch.tensor(1e-3)
+    )
+
+    # Track - should run without errors
+    outgoing = multipole.track(incoming)
+
+    # Output should have expected dimensions
+    assert outgoing is not None
+    assert outgoing.particles.shape == incoming.particles.shape
+
+
+def test_multipole_split():
+    """
+    Test that splitting a multipole into smaller segments works correctly.
+    """
+    # Create a quadrupole-like multipole
+    original = cheetah.Multipole(
+        length=torch.tensor(1.0),
+        polynom_b=torch.tensor([0.0, 1.0]),
+        fringe_quad_entrance=1,
+        fringe_quad_exit=1,
+    )
+
+    # Split into 5 pieces
+    split_elements = original.split(torch.tensor(0.2))
+
+    assert len(split_elements) == 5
+    assert isinstance(split_elements[0], cheetah.Multipole)
+
+    # Check that fringe fields are only applied at the beginning and end
+    assert split_elements[0].fringe_quad_entrance == 1
+    assert split_elements[0].fringe_quad_exit == 0
+    assert split_elements[1].fringe_quad_entrance == 0
+    assert split_elements[1].fringe_quad_exit == 0
+    assert split_elements[4].fringe_quad_entrance == 0
+    assert split_elements[4].fringe_quad_exit == 1
+
+    # Check that the total length is preserved
+    total_length = sum(element.length.item() for element in split_elements)
+    assert torch.isclose(torch.tensor(total_length), original.length)
+
+    # Check that all other properties are preserved
+    for element in split_elements:
+        assert torch.all(element.polynom_b == original.polynom_b)
+        assert torch.all(element.polynom_a == original.polynom_a)
+        assert element.max_order == original.max_order
+        assert torch.all(element.misalignment == original.misalignment)
+        assert torch.all(element.tilt == original.tilt)
+
+
+def test_multipole_parameter_validation():
+    """
+    Test that the multipole element properly validates and initializes parameters.
+    """
+    # Test with default parameters
+    multipole = cheetah.Multipole(length=torch.tensor(1.0))
+    assert multipole.polynom_a.size(0) == 2  # Should contain orders 0, 1
+    assert multipole.polynom_b.size(0) == 2
+    assert multipole.max_order == 1
+
+    # Test with custom max_order
+    multipole = cheetah.Multipole(length=torch.tensor(1.0), max_order=3)
+    assert multipole.polynom_a.size(0) == 4  # Should contain orders 0, 1, 2, 3
+    assert multipole.polynom_b.size(0) == 4
+
+    # Test with partial polynomial coefficients
+    multipole = cheetah.Multipole(
+        length=torch.tensor(1.0), polynom_b=torch.tensor([0.0, 1.0]), max_order=3
+    )
+    assert multipole.polynom_b.size(0) == 4  # Should be padded to max_order+1
+    assert multipole.polynom_b[0] == 0.0
+    assert multipole.polynom_b[1] == 1.0
+    assert multipole.polynom_b[2] == 0.0
+    assert multipole.polynom_b[3] == 0.0
+
+
+def test_multipole_vectorization():
+    """Test that a multipole with vectorized parameters works correctly."""
+    # Create a multipole with vectorized k1 values (quadrupole strengths)
+    k1_values = torch.tensor([1.0, 2.0, 3.0])
+    multipole = cheetah.Multipole(
+        length=torch.tensor(1.0),
+        polynom_b=torch.tensor([0.0, 1.0]).unsqueeze(0).repeat(3, 1)
+        * k1_values.unsqueeze(1),
+    )
+
+    # Create a simple beam
+    incoming = cheetah.ParticleBeam.from_parameters(
+        num_particles=1_000,
+        energy=torch.tensor(1e9),
+        mu_x=torch.tensor(1e-3),  # Small offset to see focusing effects
+    )
+
+    # Track through the multipole
+    outgoing = multipole.track(incoming)
+
+    # Verify the output shape matches the vectorization
+    assert outgoing.mu_px.shape == (3,)
+
+    # Verify that different k1 values produce different focusing effects
+    # Stronger k1 should result in stronger focusing
+    assert abs(outgoing.mu_px[2]) > abs(outgoing.mu_px[1]) > abs(outgoing.mu_px[0])