Lorentzian Peak Fitting

In NMR spectroscopy, the peaks that appear at characteristic chemical shifts are not Gaussian — they have heavier tails because the underlying physics (Lorentzian broadening from homogeneous line width) produces a Cauchy distribution. Fitting a Gaussian to a Lorentzian peak underestimates the tails and gives wrong width and area measurements. The Lorentzian model replaces the Gaussian's exponential decay with an algebraic `1/(1 + x²)` decay, which falls off much more slowly. In practice, real peaks are often a mix (*Voigt profile*) — but fitting each model separately and comparing the residuals tells you which is more appropriate for your data.

### The Lorentzian function

The Lorentzian has amplitude `A`, centre `x₀`, and half-width at half-maximum `γ`. Its tails decay as `1/x²` rather than the Gaussian's `exp(-x²)`.

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-10, 10, 500)

def lorentzian(x, amplitude, center, gamma):
    """gamma: half-width at half-maximum (HWHM)"""
    return amplitude * (gamma**2 / ((x - center)**2 + gamma**2))

def gaussian(x, amplitude, center, sigma):
    return amplitude * np.exp(-(x - center)**2 / (2 * sigma**2))

# Same FWHM: Lorentzian FWHM = 2*gamma, Gaussian FWHM = 2*sqrt(2*ln2)*sigma
gamma = 1.0
sigma = gamma / np.sqrt(2 * np.log(2))  # matching FWHM

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(x, lorentzian(x, 1, 0, gamma), color="tomato",    linewidth=2, label="Lorentzian (γ=1)")
ax.plot(x, gaussian(x,   1, 0, sigma), color="steelblue", linewidth=2, label=f"Gaussian (σ={sigma:.2f})")
ax.set_xlabel("x"); ax.set_ylabel("Intensity")
ax.set_title("Lorentzian vs Gaussian — same FWHM")
ax.legend()
ax.set_ylim(-0.05, 1.1)
plt.tight_layout()
plt.show()

print(f"Lorentzian FWHM: {2 * gamma:.3f}")
print(f"Gaussian FWHM:   {2 * np.sqrt(2 * np.log(2)) * sigma:.3f}")

Lorentzian FWHM: 2.000
Gaussian FWHM:   2.000

- The Lorentzian's tails are visibly heavier — the wings of the peak extend much further from the centre than the Gaussian of the same FWHM.
- For NMR, the Lorentzian width `γ` relates directly to the transverse relaxation time T₂: `γ = 1 / (π T₂)`.
- The total area of a Lorentzian is `π × A × γ`; for a Gaussian it is `A × σ × √(2π)`. Same FWHM doesn't mean same area.

### Fitting a Lorentzian to noisy data

`curve_fit` fits the Lorentzian exactly as it fits a Gaussian — only the model function changes.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

rng = np.random.default_rng(42)
x = np.linspace(-8, 8, 300)

def lorentzian(x, amplitude, center, gamma):
    return amplitude * (gamma**2 / ((x - center)**2 + gamma**2))

true_amp, true_center, true_gamma = 3.0, 0.3, 1.2
y_true  = lorentzian(x, true_amp, true_center, true_gamma)
y_noisy = y_true + rng.normal(0, 0.15, len(x))

p0 = [2.0, 0.0, 1.0]
popt, pcov = curve_fit(lorentzian, x, y_noisy, p0=p0,
                       bounds=([0, -np.inf, 0.01], [np.inf, np.inf, np.inf]))
perr = np.sqrt(np.diag(pcov))

print(f"{'Parameter':>12}  {'True':>8}  {'Fitted':>8}  {'±1σ':>8}")
print("-" * 44)
for name, true, fit, err in zip(
        ["Amplitude", "Centre", "Gamma"],
        [true_amp, true_center, true_gamma],
        popt, perr):
    print(f"{name:>12}  {true:>8.3f}  {fit:>8.3f}  {err:>8.4f}")

print(f"\nFitted FWHM: {2 * popt[2]:.3f}  (true: {2 * true_gamma:.3f})")

fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(x, y_noisy, color="lightgray",  linewidth=0.8, label="Data")
ax.plot(x, y_true,  color="black",      linewidth=1.0, linestyle="--", label="True")
ax.plot(x, lorentzian(x, *popt), color="tomato", linewidth=2, label="Lorentzian fit")
ax.set_xlabel("x"); ax.set_ylabel("Intensity")
ax.set_title("Lorentzian peak fitting")
ax.legend()
plt.tight_layout()
plt.show()

   Parameter      True    Fitted       ±1σ
--------------------------------------------
   Amplitude     3.000     2.989    0.0334
      Centre     0.300     0.310    0.0134
       Gamma     1.200     1.199    0.0190

Fitted FWHM: 2.398  (true: 2.400)

- `bounds=([0, -inf, 0.01], [inf, inf, inf])` enforces a positive amplitude and positive `γ` (HWHM must be > 0), preventing the optimiser from finding physically meaningless solutions.
- The fitted `γ` gives the FWHM directly as `2γ` — no conversion factor needed (unlike the Gaussian's `2.355σ`).
- `pcov` and `perr` work identically to the Gaussian case: the diagonal gives parameter variance, off-diagonal terms give correlations between parameters.

### Comparing Gaussian vs Lorentzian fits

Fitting both models to the same data and comparing residuals tells you which profile better describes the peaks.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

rng = np.random.default_rng(42)
x = np.linspace(-8, 8, 300)

def lorentzian(x, a, c, g):
    return a * (g**2 / ((x - c)**2 + g**2))

def gaussian(x, a, c, s):
    return a * np.exp(-(x - c)**2 / (2 * s**2))

# Simulate Lorentzian data
y_true  = lorentzian(x, 3.0, 0.3, 1.2)
y_noisy = y_true + rng.normal(0, 0.15, len(x))

popt_lor, _ = curve_fit(lorentzian, x, y_noisy, p0=[2, 0, 1],
                         bounds=([0,-np.inf,0.01],[np.inf,np.inf,np.inf]))
popt_gau, _ = curve_fit(gaussian,   x, y_noisy, p0=[2, 0, 1],
                         bounds=([0,-np.inf,0.01],[np.inf,np.inf,np.inf]))

resid_lor = y_noisy - lorentzian(x, *popt_lor)
resid_gau = y_noisy - gaussian(x,   *popt_gau)

rmse_lor = np.sqrt(np.mean(resid_lor**2))
rmse_gau = np.sqrt(np.mean(resid_gau**2))

fig, axes = plt.subplots(2, 1, figsize=(10, 6), sharex=True)
axes[0].plot(x, y_noisy, color="lightgray", linewidth=0.8, label="Data")
axes[0].plot(x, lorentzian(x, *popt_lor), color="tomato",    linewidth=2,
             label=f"Lorentzian (RMSE={rmse_lor:.4f})")
axes[0].plot(x, gaussian(x,   *popt_gau), color="steelblue", linewidth=2, linestyle="--",
             label=f"Gaussian (RMSE={rmse_gau:.4f})")
axes[0].legend(fontsize=9); axes[0].set_ylabel("Intensity")
axes[0].set_title("Fit comparison — Lorentzian data")

axes[1].plot(x, resid_lor, color="tomato",    linewidth=1.0, label="Lorentzian residuals")
axes[1].plot(x, resid_gau, color="steelblue", linewidth=1.0, linestyle="--", label="Gaussian residuals")
axes[1].axhline(0, color="black", linewidth=0.8)
axes[1].legend(fontsize=9); axes[1].set_ylabel("Residual"); axes[1].set_xlabel("x")
axes[1].set_title("Residuals — Gaussian shows systematic tails error")

plt.tight_layout()
plt.show()
print(f"Lorentzian RMSE: {rmse_lor:.4f}")
print(f"Gaussian RMSE:   {rmse_gau:.4f}")

Lorentzian RMSE: 0.1393
Gaussian RMSE:   0.2076

- Lower RMSE confirms that the Lorentzian model fits the data better when the data was generated from a Lorentzian process.
- The residuals panel is the key diagnostic: Gaussian residuals on Lorentzian data show a characteristic pattern — positive residuals in the tails, negative near the peak — because the Gaussian underestimates the tails.
- If both models give similar residuals and RMSE, either is acceptable and you should choose the simpler or more theoretically justified one.

### Conclusion

Use a Lorentzian model when your peaks arise from resonance phenomena (NMR, Raman, optical transitions) where homogeneous broadening dominates. Always plot residuals — they immediately reveal whether the chosen model is appropriate. For peaks that are neither purely Gaussian nor Lorentzian, a pseudo-Voigt (linear combination of both) is the practical middle ground used in powder diffraction and spectroscopy.

For the Gaussian alternative, see [Gaussian peak fitting](/tutorials/gaussian-peak-fitting). For detecting peak positions without fitting a parametric model, see [finding peaks with scipy](/tutorials/find-peaks-with-scipy).