diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 9e72ff2..7810758 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -14,8 +14,10 @@ The rules for spectrumkit's CHANGELOG file:
 Unreleased
 ----------
 
+- Add analytical test for dielectric spectrum of a Debye dipole (#24)
+- Fix susceptibility being a factor of 2 too large (#24)
 - Fix logarithmic binning with ``binafter=0`` (#23).
-- Require ``maicos>=0.13``, drop the ``mdacli<=0.1.35`` pin. (#20)
+- Require ``maicos>=0.13``, drop the ``mdacli<=0.1.35`` pin. (#23)
 
 v0.1.0 (2026/05/20)
 -------------------
diff --git a/examples/validation/README.rst b/examples/validation/README.rst
new file mode 100644
index 0000000..3692b03
--- /dev/null
+++ b/examples/validation/README.rst
@@ -0,0 +1,5 @@
+Validation
+----------
+
+Comparisons of the spectrum routes against systems with an analytically known spectrum.
+Each example follows one test of the test suite.
diff --git a/examples/validation/debye_relaxation.py b/examples/validation/debye_relaxation.py
new file mode 100644
index 0000000..c276234
--- /dev/null
+++ b/examples/validation/debye_relaxation.py
@@ -0,0 +1,194 @@
+r""".. _howto-debye-relaxation:
+
+Debye relaxation of an Ornstein-Uhlenbeck dipole
+------------------------------------------------
+
+The test suite checks the spectrum routes against a dipole whose susceptibility is known
+exactly. This example walks through that test, ``TestDebyeRelaxation`` in
+``tests/tests.py``, and shows what it compares.
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+import scipy.constants
+from scipy.signal import lfilter
+
+from spectrumkit import calculate_spectrum_from_current, calculate_spectrum_from_dipole
+
+BLUE, ORANGE, INK, SHADE = "#2a78d6", "#eb6834", "#0b0b0b", "#e1e0d9"
+
+# sphinx_gallery_thumbnail_number = 2
+
+# %%
+# The model
+# ~~~~~~~~~
+#
+# Each Cartesian component of the dipole :math:`\boldsymbol{M}` follows the Langevin
+# equation
+#
+# .. math::
+#
+#     \dot M_\alpha(t) = -\frac{M_\alpha(t)}{\tau}
+#     + \sqrt{\frac{2 \sigma^2}{\tau}} \, \xi_\alpha(t),
+#     \qquad \langle \xi_\alpha(t) \xi_\beta(t') \rangle
+#     = \delta_{\alpha \beta} \, \delta(t - t').
+#
+# The equation is linear and the noise is uncorrelated with the past, so multiplying by
+# :math:`M_\alpha(0)` and averaging leaves :math:`\dot C = -C / \tau` for the
+# autocorrelation:
+#
+# .. math::
+#
+#     C(t) = \langle \boldsymbol{M}(0) \cdot \boldsymbol{M}(t) \rangle
+#     = 3 \sigma^2 e^{-|t| / \tau}.
+#
+# The classical fluctuation-dissipation theorem turns this into the Debye
+# susceptibility of a system with volume :math:`V` at temperature :math:`T`,
+#
+# .. math::
+#
+#     \chi(\omega) = \frac{\chi_0}{1 - i \omega \tau}, \qquad
+#     \chi''(\omega) = \chi_0 \, \frac{\omega \tau}{1 + \omega^2 \tau^2}, \qquad
+#     \chi_0 = \frac{\langle \boldsymbol{M}^2 \rangle}{3 V \varepsilon_0 k_B T},
+#
+# with the loss peak at :math:`\omega \tau = 1`.
+#
+# Generating the dipole
+# ~~~~~~~~~~~~~~~~~~~~~
+#
+# Integrating the Langevin equation over one time step :math:`\Delta t` gives the
+# recursion
+#
+# .. math::
+#
+#     M_{n+1} = a M_n + \sigma \sqrt{1 - a^2} \, \eta_n, \qquad
+#     a = e^{-\Delta t / \tau},
+#
+# with standard normal :math:`\eta_n`. It is exact for any time step, and drawing the
+# first frame from the stationary distribution avoids an equilibration phase. The test
+# uses 2 ns of a dipole with a relaxation time of 1 ps, sampled every 10 fs. (The test
+# suite generates it with the helper ``ornstein_uhlenbeck`` in
+# ``tests/data/synthetic.py``.)
+
+dt = 0.01  # ps
+tau = 1.0  # ps
+sigma = 1.0  # e·Å
+volume = 1000.0  # Å^3
+temperature = 300.0  # K
+n_frames = 200_000
+
+rng = np.random.default_rng(0)
+a = np.exp(-dt / tau)
+noise = rng.normal(0, sigma * np.sqrt(1 - a**2), (n_frames, 3))
+noise[0] = rng.normal(0, sigma, 3)
+dipole = lfilter([1], [1, -a], noise, axis=0)
+
+# %%
+# The trajectory looks like noise, but its autocorrelation decays as the model predicts.
+
+lags = np.arange(0, int(4 * tau / dt), 5)
+acf = np.array(
+    [np.mean(np.sum(dipole[: n_frames - k] * dipole[k:], axis=1)) for k in lags]
+)
+
+time = dt * np.arange(n_frames)
+shown = time < 20
+
+fig, (ax_traj, ax_acf) = plt.subplots(1, 2, figsize=(8, 3), layout="constrained")
+ax_traj.plot(time[shown], dipole[shown, 0], color=BLUE, lw=1)
+ax_traj.set(xlabel="$t$ / ps", ylabel=r"$M_x$ / e$\cdot$Å", title="Trajectory")
+
+ax_acf.plot(lags * dt / tau, acf / acf[0], "o", color=BLUE, ms=4, label="trajectory")
+ax_acf.plot(
+    lags * dt / tau, np.exp(-lags * dt / tau), color=INK, label=r"$e^{-t/\tau}$"
+)
+ax_acf.set(xlabel=r"$t / \tau$", ylabel="$C(t) / C(0)$", title="Autocorrelation")
+ax_acf.legend(frameon=False)
+plt.show()
+
+# %%
+# Spectra from the dipole and from the current
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#
+# :func:`~spectrumkit.dielectricspectrum.calculate_spectrum_from_dipole` works on the
+# dipole, :func:`~spectrumkit.dielectricspectrum.calculate_spectrum_from_current` on the
+# current :math:`\boldsymbol{J} = \dot{\boldsymbol{M}}`, here a finite difference. Both
+# average the periodogram over 100 segments of 20 ps.
+
+dipole_route = calculate_spectrum_from_dipole(dipole, dt, volume, temperature, segs=100)
+current = np.diff(dipole, axis=0) / dt
+current_route = calculate_spectrum_from_current(
+    current, dt, volume, temperature, segs=100
+)
+
+# %%
+# The reference is the Debye susceptibility. :math:`\chi_0` uses the variance of the
+# generated dipole, so the comparison checks the shape and normalization of the
+# spectrum, not the sampling error of :math:`\langle \boldsymbol{M}^2 \rangle`.
+
+pref = scipy.constants.e**2 * scipy.constants.angstrom**2
+pref /= 3 * volume * scipy.constants.angstrom**3
+pref /= scipy.constants.k * temperature * scipy.constants.epsilon_0
+chi_0 = pref * np.sum(np.var(dipole, axis=0))
+
+
+def debye(nu):
+    """Debye susceptibility at the frequency ``nu`` in THz."""
+    return chi_0 / (1 - 2j * np.pi * nu * tau)
+
+
+# %%
+# The test compares the imaginary part in the band :math:`0.5 < \omega \tau < 2` around
+# the loss peak (shaded) and passes when the mean ratio to the Debye curve is within
+# 5 % of one (dotted lines).
+
+band = np.array([0.5, 2]) / (2 * np.pi * tau)
+nu = np.logspace(-2, np.log10(50), 300)
+
+fig, (ax, ax_ratio) = plt.subplots(
+    2, 1, sharex=True, figsize=(7, 5.5), height_ratios=(2, 1), layout="constrained"
+)
+ax.loglog(nu, debye(nu).imag, color=INK, ls="--", label="Debye")
+for res, label, color in [
+    (dipole_route, "from dipole", BLUE),
+    (current_route, "from current", ORANGE),
+]:
+    nu_res, chi_res = res["nu_binned"][1:], res["susc_binned"].imag[1:]
+    ax.loglog(nu_res, chi_res, color=color, lw=1.5, label=label)
+    ax_ratio.semilogx(nu_res, chi_res / debye(nu_res).imag, color=color, lw=1.5)
+
+for axis in (ax, ax_ratio):
+    axis.axvspan(*band, color=SHADE, zorder=0)
+ax_ratio.axhline(1, color=INK, ls="--", lw=1)
+for limit in (0.95, 1.05):
+    ax_ratio.axhline(limit, color=INK, ls=":", lw=1)
+
+ax.set(ylabel=r"$\chi''$", title="Dielectric loss")
+ax.legend(frameon=False)
+ax_ratio.set(xlabel=r"$\nu$ / THz", ylabel=r"$\chi'' / \chi''_\mathrm{Debye}$")
+ax_ratio.set_ylim(0.6, 1.4)
+plt.show()
+
+# %%
+# Inside the band, the mean ratios are
+
+for name, res in [("dipole", dipole_route), ("current", current_route)]:
+    in_band = (res["nu"] > band[0]) & (res["nu"] < band[1])
+    ratio = np.mean(res["susc"].imag[in_band] / debye(res["nu"][in_band]).imag)
+    print(f"from {name}: {ratio:.3f}")
+
+# %%
+# Outside the band, the curves deviate in ways the test leaves out:
+#
+# - Above about 10 THz the dipole route rises. This is aliasing, not an error of the
+#   route: sampled every 10 fs, the series has a higher spectrum near the Nyquist
+#   frequency of 50 THz than the continuous process.
+# - Between 1 and 10 THz the dipole route lies about 5 % high, and the current route
+#   overshoots at its lowest frequencies. Both are spectral leakage: the periodogram
+#   treats each 20 ps segment as periodic, and the jump between its ends adds a
+#   spurious contribution.
+#
+# The test runs the same check for
+# :func:`~spectrumkit.dielectricspectrum.WienerKhinchin` with
+# ``spectrum_type="dipole"`` and ``"flux"``. It covers only the imaginary part; the real
+# part comes from a numerical Kramers-Kronig transform and is not tested here.
diff --git a/src/spectrumkit/dielectricspectrum.py b/src/spectrumkit/dielectricspectrum.py
index f50b98c..2c2fc1e 100644
--- a/src/spectrumkit/dielectricspectrum.py
+++ b/src/spectrumkit/dielectricspectrum.py
@@ -156,9 +156,10 @@ def calculate_spectrum_from_dipole(
     dsusc.real = np.sqrt(dsusc.real)
     dsusc.imag = np.sqrt(dsusc.imag)
 
-    # Normalization factor
-    susc *= pref / (seglen * segs * dt)
-    dsusc *= pref / (seglen * segs * dt)
+    # Normalization: 1/T turns |FT(P)|² into the power spectrum S(ω), and the
+    # classical fluctuation-dissipation theorem gives χ''(ω) ∝ ω S(ω) / 2.
+    susc *= pref / (2 * seglen * segs * dt)
+    dsusc *= pref / (2 * seglen * segs * dt)
 
     # Convert to THz
     nu = nu / (2 * np.pi)
@@ -328,9 +329,10 @@ def calculate_spectrum_from_current(
     dsusc.real = np.sqrt(dsusc.real)
     dsusc.imag = np.sqrt(dsusc.imag)
 
-    # Normalization factor
-    susc *= pref / (seglen * segs * dt)
-    dsusc *= pref / (seglen * segs * dt)
+    # Normalization: 1/T turns |FT(P)|² into the power spectrum S(ω), and the
+    # classical fluctuation-dissipation theorem gives χ''(ω) ∝ ω S(ω) / 2.
+    susc *= pref / (2 * seglen * segs * dt)
+    dsusc *= pref / (2 * seglen * segs * dt)
 
     # Convert to THz
     nu = nu / (2 * np.pi)
@@ -826,7 +828,8 @@ def WienerKhinchin(time, timeseries, volume, temperature, spectrum_type="dipole"
     pref /= 3 * volume * scipy.constants.angstrom**3
     pref /= scipy.constants.epsilon_0 * scipy.constants.k * temperature
 
-    dielectricspectrum_imag = powerspectrum * pref
+    # Classical fluctuation-dissipation theorem: χ''(ω) = ω S_PP(ω) / 2
+    dielectricspectrum_imag = powerspectrum * pref / 2
 
     if spectrum_type == "dipole":
         dielectricspectrum_imag *= nu
diff --git a/tests/data/synthetic.py b/tests/data/synthetic.py
new file mode 100644
index 0000000..cc2d599
--- /dev/null
+++ b/tests/data/synthetic.py
@@ -0,0 +1,69 @@
+#
+# Copyright (c) 2026 Authors and contributors
+# (see the AUTHORS.rst file for the full list of names)
+#
+# Released under the GNU Public Licence, v3 or any higher version
+# SPDX-License-Identifier: GPL-3.0-or-later
+"""Synthetic time series with analytically known spectra."""
+
+import numpy as np
+from scipy.signal import lfilter
+
+
+def ornstein_uhlenbeck(
+    n_frames: int,
+    n_dim: int = 3,
+    dt: float = 1.0,
+    tau: float = 1.0,
+    sigma: float = 1.0,
+    seed: int = 0,
+) -> np.ndarray:
+    r"""Ornstein-Uhlenbeck process, e.g. a dipole with Debye relaxation.
+
+    Each of the ``n_dim`` independent components has zero mean, variance
+    ``sigma**2`` and an exponential autocorrelation with correlation time ``tau``.
+    Its spectra are known in closed form (see Notes), which makes it a reference for
+    testing spectrum estimators. The process is sampled exactly as an AR(1) recursion
+    and starts from the stationary distribution.
+
+    Parameters
+    ----------
+    n_frames : int
+        Number of time steps.
+    n_dim : int, default=3
+        Number of independent components.
+    dt : float, default=1.0
+        Time step.
+    tau : float, default=1.0
+        Correlation time, in the same unit as ``dt``.
+    sigma : float, default=1.0
+        Standard deviation of each component.
+    seed : int, default=0
+        Seed of the random number generator.
+
+    Returns
+    -------
+    numpy.ndarray
+        Time series with shape (n_frames, n_dim).
+
+    Notes
+    -----
+    Each component solves the Langevin equation
+
+    .. math::
+
+        \dot x(t) = -\frac{x(t)}{\tau} + \sqrt{\frac{2 \sigma^2}{\tau}} \, \xi(t),
+        \qquad \langle \xi(t) \xi(t') \rangle = \delta(t - t').
+
+    The equation is linear and the noise is uncorrelated with the past, so multiplying
+    by :math:`x(0)` and averaging leaves :math:`\dot C = -C / \tau` for the
+    autocorrelation, i.e. :math:`C(t) = \sigma^2 e^{-|t| / \tau}`. Its Fourier
+    transform is a Lorentzian, and the fluctuation-dissipation theorem turns it into
+    the Debye susceptibility :math:`\chi_0 / (1 - i \omega \tau)`.
+
+    """
+    rng = np.random.default_rng(seed)
+    a = np.exp(-dt / tau)
+    noise = rng.normal(0, sigma * np.sqrt(1 - a**2), (n_frames, n_dim))
+    noise[0] = rng.normal(0, sigma, n_dim)  # start from the stationary distribution
+    return lfilter([1], [1, -a], noise, axis=0)
diff --git a/tests/tests.py b/tests/tests.py
index 4ac43aa..584844a 100644
--- a/tests/tests.py
+++ b/tests/tests.py
@@ -11,12 +11,19 @@
 import MDAnalysis as mda
 import numpy as np
 import pytest
+import scipy.constants
 from numpy.testing import assert_allclose
 
-from spectrumkit import DielectricSpectrum, calculate_spectrum_from_dipole
+from spectrumkit import (
+    DielectricSpectrum,
+    WienerKhinchin,
+    calculate_spectrum_from_current,
+    calculate_spectrum_from_dipole,
+)
 
 sys.path.append(str(Path(__file__).parents[1]))
 from data import WATER_TPR_NPT, WATER_TRR_NPT
+from data.synthetic import ornstein_uhlenbeck
 
 
 class TestDielectricSpectrum:
@@ -91,8 +98,8 @@ def test_output(self, ag, monkeypatch, tmp_path):
 
         V = 1559814.4
         nu = [0.0, 0.5, 1.0]
-        susc = [2.7 + 0.0j, 4.3 + 5.660148j, 2.7 + 5.6j]
-        dsusc = [0.5 + 0.0j, 0.7 + 0.765955j, 0.5 + 0.9j]
+        susc = [1.35 + 0.0j, 2.15 + 2.830074j, 1.35 + 2.8j]
+        dsusc = [0.25 + 0.0j, 0.35 + 0.3829775j, 0.25 + 0.45j]
 
         assert_allclose(ds.V, V, rtol=1e-1)
         assert_allclose(ds.results.nu, nu, rtol=1)
@@ -177,3 +184,97 @@ def test_consistency_with_dielectric_spectrum(self, monkeypatch, tmp_path):
         assert_allclose(ds.results.nu, results["nu"], rtol=1e-10)
         assert_allclose(ds.results.susc, results["susc"], rtol=1e-10)
         assert_allclose(ds.results.dsusc, results["dsusc"], rtol=1e-10)
+
+
+class TestDebyeRelaxation:
+    """Compare the spectrum routes to an Ornstein-Uhlenbeck dipole.
+
+    A dipole with exponential autocorrelation <P(0)P(t)> = <P²> exp(-t/τ) has the
+    Debye susceptibility χ(ω) = χ_0 / (1 - iωτ) with χ_0 = <P²> / (3 V k_B T ε_0).
+    """
+
+    dt = 0.01  # ps
+    tau = 1.0  # ps
+    volume = 1000.0  # Å^3
+    temperature = 300.0  # K
+
+    @pytest.fixture(scope="class")
+    @classmethod
+    def dipole(cls):
+        """Ornstein-Uhlenbeck dipole time series (n_frames, 3) in e·Å."""
+        return ornstein_uhlenbeck(200_000, dt=cls.dt, tau=cls.tau)
+
+    def debye(self, dipole, nu):
+        """Analytic Debye susceptibility for the sampled dipole variance."""
+        pref = scipy.constants.e**2 * scipy.constants.angstrom**2
+        pref /= 3 * self.volume * scipy.constants.angstrom**3
+        pref /= scipy.constants.k * self.temperature * scipy.constants.epsilon_0
+        chi_0 = pref * np.sum(np.var(dipole, axis=0))
+        return chi_0 / (1 - 2j * np.pi * nu * self.tau)
+
+    def band(self, nu):
+        """Mask for 0.5 < ωτ < 2, around the loss peak."""
+        omega_tau = 2 * np.pi * nu * self.tau
+        return (omega_tau > 0.5) & (omega_tau < 2)
+
+    @pytest.mark.parametrize(("n_dim", "sigma"), [(1, 1.0), (3, 2.5)])
+    def test_ornstein_uhlenbeck(self, n_dim, sigma):
+        """The reference process has variance σ² and decays as exp(-t/τ).
+
+        Checked with plain numpy, independent of spectrumkit.
+        """
+        x = ornstein_uhlenbeck(
+            200_000, n_dim=n_dim, dt=self.dt, tau=self.tau, sigma=sigma
+        )
+        assert x.shape == (200_000, n_dim)
+
+        x = x - x.mean(axis=0)
+        lags = np.array([0, 50, 100, 200])
+        corr = np.array(
+            [np.mean(np.sum(x[: len(x) - k] * x[k:], axis=1)) for k in lags]
+        )
+        assert_allclose(corr[0] / n_dim, sigma**2, rtol=0.1)
+        assert_allclose(corr / corr[0], np.exp(-lags * self.dt / self.tau), atol=0.05)
+
+    def test_periodogram_dipole(self, dipole):
+        """Periodogram of the dipole reproduces Im(χ)."""
+        res = calculate_spectrum_from_dipole(
+            dipole, self.dt, self.volume, self.temperature, segs=100, nobin=True
+        )
+        m = self.band(res["nu"])
+        expected = self.debye(dipole, res["nu"][m])
+        assert_allclose(np.mean(res["susc"].imag[m] / expected.imag), 1, rtol=0.05)
+
+    def test_periodogram_current(self, dipole):
+        """Periodogram of the current reproduces Im(χ)."""
+        current = np.diff(dipole, axis=0) / self.dt
+        res = calculate_spectrum_from_current(
+            current, self.dt, self.volume, self.temperature, segs=100, nobin=True
+        )
+        m = self.band(res["nu"])
+        expected = self.debye(dipole, res["nu"][m])
+        assert_allclose(np.mean(res["susc"].imag[m] / expected.imag), 1, rtol=0.05)
+
+    @pytest.mark.parametrize("spectrum_type", ["dipole", "flux"])
+    def test_wiener_khinchin(self, dipole, spectrum_type):
+        """WienerKhinchin, averaged over blocks, reproduces Im(χ)."""
+        series = dipole
+        if spectrum_type == "flux":
+            series = np.diff(dipole, axis=0) / self.dt
+        block = 2000
+        time = self.dt * np.arange(block)
+        results = [
+            WienerKhinchin(
+                time,
+                series[i : i + block],
+                self.volume,
+                self.temperature,
+                spectrum_type=spectrum_type,
+            )
+            for i in range(0, len(series) - block + 1, block)
+        ]
+        nu = results[0][0] / (2 * np.pi)
+        susc = np.mean([s for _, s in results], axis=0)
+        m = self.band(nu)
+        expected = self.debye(dipole, nu[m])
+        assert_allclose(np.mean(susc.imag[m] / expected.imag), 1, rtol=0.05)