diff --git a/autogalaxy/profiles/light/standard/__init__.py b/autogalaxy/profiles/light/standard/__init__.py
index 7de832bc..710e1575 100644
--- a/autogalaxy/profiles/light/standard/__init__.py
+++ b/autogalaxy/profiles/light/standard/__init__.py
@@ -13,6 +13,7 @@
     ElsonFreeFall,
     ElsonFreeFallSph,
 )
+from .multipole import SersicMultipole, GaussianMultipole
 from .shapelets.polar import ShapeletPolarSph, ShapeletPolar
 from .shapelets.cartesian import ShapeletCartesianSph, ShapeletCartesian
 from .shapelets.exponential import ShapeletExponentialSph, ShapeletExponential
diff --git a/autogalaxy/profiles/light/standard/multipole.py b/autogalaxy/profiles/light/standard/multipole.py
new file mode 100644
index 00000000..d66fc9ef
--- /dev/null
+++ b/autogalaxy/profiles/light/standard/multipole.py
@@ -0,0 +1,246 @@
+"""
+Multipole light profiles.
+
+These profiles extend the standard Sersic and Gaussian light profiles by perturbing the
+eccentric radius with a Fourier series at orders m=3 and m=4 before evaluating the
+underlying radial profile:
+
+    r_perturbed = r * (1 + c_3 cos(3 theta) + s_3 sin(3 theta)
+                          + c_4 cos(4 theta) + s_4 sin(4 theta))
+
+where ``theta`` is the polar angle in the profile's elliptical reference frame. With all
+multipole components set to zero (the default), the profile reduces exactly to its base
+class (``Sersic`` or ``Gaussian``).
+
+The (cos, sin) components are passed as a 2-tuple ``multipole_<m>_comps = (c, s)``,
+matching the parameter convention used by ``EllipseMultipole`` and ``PowerLawMultipole``.
+
+Note: this is the "elliptical multipole" flavour — the perturbation rides on the
+already-elliptical ``grid_radii``, so the deformation follows the profile's ellipticity.
+This is distinct from the polar-frame perturbation in ``EllipseMultipole``.
+"""
+
+from numpy import seterr
+from typing import Optional, Tuple
+
+import numpy as np
+
+import autoarray as aa
+
+from autogalaxy.profiles.light.decorators import check_operated_only
+from autogalaxy.profiles.light.standard.gaussian import Gaussian
+from autogalaxy.profiles.light.standard.sersic import Sersic
+
+
+class _LightProfileMultipoleMixin:
+    """
+    Shared `perturbed_radii_from` helper for light profiles with m=3 and m=4 Fourier
+    perturbations applied to the eccentric radius.
+
+    Concrete subclasses must set ``multipole_3_comps`` and ``multipole_4_comps`` on
+    ``self`` and inherit from a profile that provides ``eccentric_radii_grid_from``.
+    """
+
+    multipole_3_comps: Tuple[float, float]
+    multipole_4_comps: Tuple[float, float]
+
+    def perturbed_radii_from(
+        self,
+        grid: aa.type.Grid2DLike,
+        xp=np,
+        **kwargs,
+    ) -> np.ndarray:
+        """
+        Returns the eccentric radii of ``grid`` perturbed by the m=3 and m=4 Fourier
+        multipoles, as a raw backend array (numpy or jax.numpy).
+
+        The result is floored at 1e-8 to keep the Sersic ``r ** (1/n)`` evaluation
+        finite when a large multipole component would otherwise drive the radius
+        through zero.
+        """
+        grid_radii = self.eccentric_radii_grid_from(grid=grid, xp=xp, **kwargs)
+        grid_radii_array = (
+            grid_radii.array if hasattr(grid_radii, "array") else grid_radii
+        )
+        y = grid.array[:, 0]
+        x = grid.array[:, 1]
+        theta = xp.arctan2(y, x)
+        c3, s3 = self.multipole_3_comps
+        c4, s4 = self.multipole_4_comps
+        perturbation = (
+            1.0
+            + c3 * xp.cos(3.0 * theta)
+            + s3 * xp.sin(3.0 * theta)
+            + c4 * xp.cos(4.0 * theta)
+            + s4 * xp.sin(4.0 * theta)
+        )
+        return xp.maximum(xp.multiply(grid_radii_array, perturbation), 1e-8)
+
+
+class SersicMultipole(_LightProfileMultipoleMixin, Sersic):
+    def __init__(
+        self,
+        centre: Tuple[float, float] = (0.0, 0.0),
+        ell_comps: Tuple[float, float] = (0.0, 0.0),
+        intensity: float = 0.1,
+        effective_radius: float = 0.6,
+        sersic_index: float = 4.0,
+        multipole_3_comps: Tuple[float, float] = (0.0, 0.0),
+        multipole_4_comps: Tuple[float, float] = (0.0, 0.0),
+    ):
+        """
+        The elliptical Sersic light profile with m=3 and m=4 Fourier multipole
+        perturbations on the eccentric radius.
+
+        Parameters
+        ----------
+        centre
+            The (y,x) arc-second coordinates of the profile centre.
+        ell_comps
+            The first and second ellipticity components of the elliptical coordinate
+            system. The multipole perturbation is applied to the eccentric radius and
+            therefore follows this ellipticity.
+        intensity
+            Overall intensity normalisation of the light profile.
+        effective_radius
+            The circular radius containing half the light of this profile (in the
+            unperturbed limit ``multipole_*_comps = (0, 0)``).
+        sersic_index
+            Controls the concentration of the profile.
+        multipole_3_comps
+            The ``(cos, sin)`` components of the m=3 Fourier perturbation. Defaults to
+            ``(0.0, 0.0)`` which reduces the profile to ``Sersic``.
+        multipole_4_comps
+            The ``(cos, sin)`` components of the m=4 Fourier perturbation. Defaults to
+            ``(0.0, 0.0)`` which reduces the profile to ``Sersic``.
+        """
+        super().__init__(
+            centre=centre,
+            ell_comps=ell_comps,
+            intensity=intensity,
+            effective_radius=effective_radius,
+            sersic_index=sersic_index,
+        )
+        self.multipole_3_comps = multipole_3_comps
+        self.multipole_4_comps = multipole_4_comps
+
+    def image_2d_via_radii_from(
+        self, grid_radii: np.ndarray, xp=np, **kwargs
+    ) -> np.ndarray:
+        """
+        Returns the 2D Sersic image evaluated at the input radial values.
+
+        Unlike ``Sersic.image_2d_via_radii_from``, this override accepts a raw backend
+        array (the output of ``perturbed_radii_from``) rather than an autoarray-wrapped
+        grid, since the perturbation step strips the wrapper.
+        """
+        seterr(all="ignore")
+        return xp.multiply(
+            self._intensity,
+            xp.exp(
+                xp.multiply(
+                    -self.sersic_constant,
+                    xp.add(
+                        xp.power(
+                            xp.divide(grid_radii, self.effective_radius),
+                            1.0 / self.sersic_index,
+                        ),
+                        -1,
+                    ),
+                )
+            ),
+        )
+
+    @aa.over_sample
+    @aa.decorators.to_array
+    @check_operated_only
+    @aa.decorators.transform
+    def image_2d_from(
+        self,
+        grid: aa.type.Grid2DLike,
+        xp=np,
+        operated_only: Optional[bool] = None,
+        **kwargs,
+    ) -> aa.Array2D:
+        """
+        Returns the 2D image of the multipole-perturbed Sersic profile.
+        """
+        perturbed_radii = self.perturbed_radii_from(grid=grid, xp=xp, **kwargs)
+        return self.image_2d_via_radii_from(grid_radii=perturbed_radii, xp=xp, **kwargs)
+
+
+class GaussianMultipole(_LightProfileMultipoleMixin, Gaussian):
+    def __init__(
+        self,
+        centre: Tuple[float, float] = (0.0, 0.0),
+        ell_comps: Tuple[float, float] = (0.0, 0.0),
+        intensity: float = 0.1,
+        sigma: float = 1.0,
+        multipole_3_comps: Tuple[float, float] = (0.0, 0.0),
+        multipole_4_comps: Tuple[float, float] = (0.0, 0.0),
+    ):
+        """
+        The elliptical Gaussian light profile with m=3 and m=4 Fourier multipole
+        perturbations on the eccentric radius.
+
+        Parameters
+        ----------
+        centre
+            The (y,x) arc-second coordinates of the profile centre.
+        ell_comps
+            The first and second ellipticity components of the elliptical coordinate
+            system. The multipole perturbation is applied to the eccentric radius and
+            therefore follows this ellipticity.
+        intensity
+            Overall intensity normalisation of the light profile.
+        sigma
+            The sigma value of the Gaussian.
+        multipole_3_comps
+            The ``(cos, sin)`` components of the m=3 Fourier perturbation. Defaults to
+            ``(0.0, 0.0)`` which reduces the profile to ``Gaussian``.
+        multipole_4_comps
+            The ``(cos, sin)`` components of the m=4 Fourier perturbation. Defaults to
+            ``(0.0, 0.0)`` which reduces the profile to ``Gaussian``.
+        """
+        super().__init__(
+            centre=centre, ell_comps=ell_comps, intensity=intensity, sigma=sigma
+        )
+        self.multipole_3_comps = multipole_3_comps
+        self.multipole_4_comps = multipole_4_comps
+
+    def image_2d_via_radii_from(
+        self, grid_radii: np.ndarray, xp=np, **kwargs
+    ) -> np.ndarray:
+        """
+        Returns the 2D Gaussian image evaluated at the input radial values.
+
+        Unlike ``Gaussian.image_2d_via_radii_from``, this override accepts a raw backend
+        array (the output of ``perturbed_radii_from``) rather than an autoarray-wrapped
+        grid, since the perturbation step strips the wrapper.
+        """
+        return xp.multiply(
+            self._intensity,
+            xp.exp(
+                -0.5
+                * xp.square(
+                    xp.divide(grid_radii, self.sigma / xp.sqrt(self.axis_ratio(xp)))
+                )
+            ),
+        )
+
+    @aa.over_sample
+    @aa.decorators.to_array
+    @check_operated_only
+    @aa.decorators.transform
+    def image_2d_from(
+        self,
+        grid: aa.type.Grid2DLike,
+        xp=np,
+        operated_only: Optional[bool] = None,
+        **kwargs,
+    ) -> aa.Array2D:
+        """
+        Returns the 2D image of the multipole-perturbed Gaussian profile.
+        """
+        perturbed_radii = self.perturbed_radii_from(grid=grid, xp=xp, **kwargs)
+        return self.image_2d_via_radii_from(grid_radii=perturbed_radii, xp=xp, **kwargs)
diff --git a/docs/api/light.rst b/docs/api/light.rst
index fcaf8935..712f1359 100644
--- a/docs/api/light.rst
+++ b/docs/api/light.rst
@@ -16,8 +16,10 @@ Standard parametric light profiles whose ``intensity`` is a free parameter of th
 
    Gaussian
    GaussianSph
+   GaussianMultipole
    Sersic
    SersicSph
+   SersicMultipole
    Exponential
    ExponentialSph
    DevVaucouleurs
diff --git a/test_autogalaxy/profiles/light/standard/test_multipole.py b/test_autogalaxy/profiles/light/standard/test_multipole.py
new file mode 100644
index 00000000..3f283418
--- /dev/null
+++ b/test_autogalaxy/profiles/light/standard/test_multipole.py
@@ -0,0 +1,221 @@
+from __future__ import division, print_function
+
+import numpy as np
+import pytest
+
+import autogalaxy as ag
+
+
+grid = ag.Grid2DIrregular([[1.0, 0.0], [0.5, 0.5], [2.0, 1.0], [-1.0, 0.3]])
+
+
+def test__sersic_multipole__zero_perturbation__matches_sersic_exactly():
+    kwargs = dict(
+        centre=(0.0, 0.0),
+        ell_comps=(0.1, 0.2),
+        intensity=2.5,
+        effective_radius=1.3,
+        sersic_index=3.0,
+    )
+    base = ag.lp.Sersic(**kwargs)
+    mp = ag.lp.SersicMultipole(
+        **kwargs,
+        multipole_3_comps=(0.0, 0.0),
+        multipole_4_comps=(0.0, 0.0),
+    )
+
+    image_base = np.asarray(base.image_2d_from(grid=grid))
+    image_mp = np.asarray(mp.image_2d_from(grid=grid))
+
+    assert image_base == pytest.approx(image_mp, rel=1e-12, abs=1e-12)
+
+
+def test__sersic_multipole__defaults__match_default_sersic():
+    image_base = np.asarray(ag.lp.Sersic().image_2d_from(grid=grid))
+    image_mp = np.asarray(ag.lp.SersicMultipole().image_2d_from(grid=grid))
+
+    assert image_base == pytest.approx(image_mp, rel=1e-12, abs=1e-12)
+
+
+def test__sersic_multipole__m4_only__has_fourfold_symmetry_on_circular_base():
+    mp = ag.lp.SersicMultipole(
+        centre=(0.0, 0.0),
+        ell_comps=(0.0, 0.0),
+        intensity=1.0,
+        effective_radius=0.6,
+        sersic_index=4.0,
+        multipole_4_comps=(0.05, 0.02),
+    )
+    pts = np.array([[1.0, 0.5], [0.3, 1.2], [-0.5, 0.8]])
+    rotated = np.column_stack([-pts[:, 1], pts[:, 0]])
+
+    image_a = np.asarray(mp.image_2d_from(grid=ag.Grid2DIrregular(pts)))
+    image_b = np.asarray(mp.image_2d_from(grid=ag.Grid2DIrregular(rotated)))
+
+    assert image_a == pytest.approx(image_b, rel=1e-10, abs=1e-12)
+
+
+def test__sersic_multipole__m3_only__has_threefold_symmetry_on_circular_base():
+    mp = ag.lp.SersicMultipole(
+        centre=(0.0, 0.0),
+        ell_comps=(0.0, 0.0),
+        intensity=1.0,
+        effective_radius=0.6,
+        sersic_index=4.0,
+        multipole_3_comps=(0.03, 0.04),
+    )
+    c120 = np.cos(2.0 * np.pi / 3.0)
+    s120 = np.sin(2.0 * np.pi / 3.0)
+    pts = np.array([[1.0, 0.5], [0.3, 1.2], [-0.5, 0.8]])
+    rotated = np.column_stack(
+        [
+            c120 * pts[:, 0] - s120 * pts[:, 1],
+            s120 * pts[:, 0] + c120 * pts[:, 1],
+        ]
+    )
+
+    image_a = np.asarray(mp.image_2d_from(grid=ag.Grid2DIrregular(pts)))
+    image_b = np.asarray(mp.image_2d_from(grid=ag.Grid2DIrregular(rotated)))
+
+    assert image_a == pytest.approx(image_b, rel=1e-10, abs=1e-12)
+
+
+def test__sersic_multipole__centre_translation__shifts_image():
+    mp_origin = ag.lp.SersicMultipole(
+        centre=(0.0, 0.0),
+        ell_comps=(0.0, 0.0),
+        intensity=1.0,
+        effective_radius=0.6,
+        sersic_index=4.0,
+        multipole_4_comps=(0.05, 0.0),
+    )
+    mp_translated = ag.lp.SersicMultipole(
+        centre=(0.3, -0.2),
+        ell_comps=(0.0, 0.0),
+        intensity=1.0,
+        effective_radius=0.6,
+        sersic_index=4.0,
+        multipole_4_comps=(0.05, 0.0),
+    )
+    pts = np.array([[1.0, 0.5], [-0.4, 0.7]])
+    pts_shifted = pts + np.array([0.3, -0.2])
+
+    image_origin = np.asarray(mp_origin.image_2d_from(grid=ag.Grid2DIrregular(pts)))
+    image_translated = np.asarray(
+        mp_translated.image_2d_from(grid=ag.Grid2DIrregular(pts_shifted))
+    )
+
+    assert image_origin == pytest.approx(image_translated, rel=1e-10, abs=1e-12)
+
+
+def test__sersic_multipole__nonzero_multipole__differs_from_sersic():
+    sersic = ag.lp.Sersic(
+        centre=(0.0, 0.0),
+        ell_comps=(0.0, 0.0),
+        intensity=1.0,
+        effective_radius=0.6,
+        sersic_index=4.0,
+    )
+    mp = ag.lp.SersicMultipole(
+        centre=(0.0, 0.0),
+        ell_comps=(0.0, 0.0),
+        intensity=1.0,
+        effective_radius=0.6,
+        sersic_index=4.0,
+        multipole_4_comps=(0.1, 0.0),
+    )
+    diff = np.asarray(mp.image_2d_from(grid=grid)) - np.asarray(
+        sersic.image_2d_from(grid=grid)
+    )
+
+    assert np.max(np.abs(diff)) > 1e-5
+
+
+def test__sersic_multipole__finite_and_nonnegative_under_large_perturbation():
+    mp = ag.lp.SersicMultipole(
+        centre=(0.0, 0.0),
+        ell_comps=(0.1, 0.1),
+        intensity=1.0,
+        effective_radius=0.6,
+        sersic_index=4.0,
+        multipole_3_comps=(0.2, -0.15),
+        multipole_4_comps=(0.3, 0.1),
+    )
+    uniform_grid = ag.Grid2D.uniform(shape_native=(11, 11), pixel_scales=0.2)
+    image = np.asarray(mp.image_2d_from(grid=uniform_grid))
+
+    assert np.all(np.isfinite(image))
+    assert np.all(image >= 0.0)
+
+
+def test__gaussian_multipole__zero_perturbation__matches_gaussian_exactly():
+    kwargs = dict(centre=(0.0, 0.0), ell_comps=(0.1, 0.2), intensity=1.5, sigma=0.7)
+    base = ag.lp.Gaussian(**kwargs)
+    mp = ag.lp.GaussianMultipole(
+        **kwargs,
+        multipole_3_comps=(0.0, 0.0),
+        multipole_4_comps=(0.0, 0.0),
+    )
+
+    image_base = np.asarray(base.image_2d_from(grid=grid))
+    image_mp = np.asarray(mp.image_2d_from(grid=grid))
+
+    assert image_base == pytest.approx(image_mp, rel=1e-12, abs=1e-12)
+
+
+def test__gaussian_multipole__defaults__match_default_gaussian():
+    image_base = np.asarray(ag.lp.Gaussian().image_2d_from(grid=grid))
+    image_mp = np.asarray(ag.lp.GaussianMultipole().image_2d_from(grid=grid))
+
+    assert image_base == pytest.approx(image_mp, rel=1e-12, abs=1e-12)
+
+
+def test__gaussian_multipole__m4_only__has_fourfold_symmetry_on_circular_base():
+    mp = ag.lp.GaussianMultipole(
+        centre=(0.0, 0.0),
+        ell_comps=(0.0, 0.0),
+        intensity=1.0,
+        sigma=0.5,
+        multipole_4_comps=(0.05, 0.02),
+    )
+    pts = np.array([[1.0, 0.5], [0.3, 1.2], [-0.5, 0.8]])
+    rotated = np.column_stack([-pts[:, 1], pts[:, 0]])
+
+    image_a = np.asarray(mp.image_2d_from(grid=ag.Grid2DIrregular(pts)))
+    image_b = np.asarray(mp.image_2d_from(grid=ag.Grid2DIrregular(rotated)))
+
+    assert image_a == pytest.approx(image_b, rel=1e-10, abs=1e-12)
+
+
+def test__gaussian_multipole__nonzero_multipole__differs_from_gaussian():
+    gauss = ag.lp.Gaussian(
+        centre=(0.0, 0.0), ell_comps=(0.0, 0.0), intensity=1.0, sigma=0.5
+    )
+    mp = ag.lp.GaussianMultipole(
+        centre=(0.0, 0.0),
+        ell_comps=(0.0, 0.0),
+        intensity=1.0,
+        sigma=0.5,
+        multipole_4_comps=(0.1, 0.0),
+    )
+    diff = np.asarray(mp.image_2d_from(grid=grid)) - np.asarray(
+        gauss.image_2d_from(grid=grid)
+    )
+
+    assert np.max(np.abs(diff)) > 1e-5
+
+
+def test__gaussian_multipole__finite_and_nonnegative_under_large_perturbation():
+    mp = ag.lp.GaussianMultipole(
+        centre=(0.0, 0.0),
+        ell_comps=(0.1, 0.1),
+        intensity=1.0,
+        sigma=0.5,
+        multipole_3_comps=(0.2, -0.15),
+        multipole_4_comps=(0.3, 0.1),
+    )
+    uniform_grid = ag.Grid2D.uniform(shape_native=(11, 11), pixel_scales=0.2)
+    image = np.asarray(mp.image_2d_from(grid=uniform_grid))
+
+    assert np.all(np.isfinite(image))
+    assert np.all(image >= 0.0)