When you measure the FWHM of a diffraction peak and plug it into the Scherrer equation, you get a number. That number is not the crystallite size — it is the crystallite size assuming zero microstrain. If your sample has lattice distortions (and most real materials do), the Scherrer result underestimates the true crystallite size because part of the broadening you measured comes from strain, not from finite crystal dimensions. The Williamson-Hall (W-H) method fixes this by using the angular dependence of broadening to separate the two contributions.
This is the third and final post in my XRD processing series. We have covered peak finding and peak fitting — now we use the FWHM values from those fits to extract both crystallite size \(D\) and microstrain \(\varepsilon\) simultaneously.
1. The Limitation of Scherrer Alone
The Scherrer equation gives the crystallite size from a single peak:
where \(K\) is the shape factor (typically 0.9 for spherical crystallites), \(\lambda\) is the X-ray wavelength, \(\beta_D\) is the size-broadened FWHM in radians, and \(\theta\) is the Bragg angle. The problem is that the observed broadening \(\beta\) has two contributions:
Size Broadening
Inversely proportional to crystallite size. The angular dependence follows \(\beta_D \propto 1/\cos\theta\). Dominant for nanocrystalline materials with \(D < 100\) nm.
Strain Broadening
Proportional to the microstrain \(\varepsilon\) and follows \(\beta_\varepsilon \propto \tan\theta\). Arises from non-uniform lattice distortions, dislocations, point defects, grain boundaries, and compositional gradients.
The Scherrer equation lumps both into a single \(\beta\) and attributes all broadening to size. For well-annealed bulk ceramics with negligible strain this is fine, but for nanoparticles — especially sol-gel or mechanically milled samples — the strain contribution can be substantial. Applying Scherrer naively will give you a crystallite size that is systematically too small.
2. The Williamson-Hall Equation
Williamson and Hall (1953) recognized that size and strain broadening have different angular dependencies. If we assume both contributions add linearly (valid for Lorentzian profiles), the total broadening is:
Multiplying both sides by \(\cos\theta\):
This is a linear equation in the form \(y = b + mx\), where \(y = \beta\cos\theta\), \(x = 4\sin\theta\), the slope \(m = \varepsilon\), and the y-intercept \(b = K\lambda/D\). By plotting all peaks on this graph and performing a linear regression, you extract both \(D\) and \(\varepsilon\) simultaneously.
Size broadening scales as \(1/\cos\theta\) while strain broadening scales as \(\tan\theta\). After multiplying by \(\cos\theta\), the size term becomes a constant (independent of angle) while the strain term becomes proportional to \(\sin\theta\). The plot separates the two because they respond differently to the Bragg angle — and that is the core insight of the Williamson-Hall method.
3. Step-by-Step Procedure
Collect FWHM values from peak fitting
For each diffraction peak, you need the fitted \(2\theta_0\) (peak center) and \(\beta_\text{obs}\) (FWHM from the pseudo-Voigt fit). Convert the FWHM from degrees to radians: \(\beta_\text{rad} = \beta_\text{deg} \times \pi/180\).
Subtract instrumental broadening
Measure a standard reference (LaB6, Si, CeO2) under identical conditions to obtain \(\beta_\text{inst}\) at each \(2\theta\). For Lorentzian profiles, use \(\beta = \beta_\text{obs} - \beta_\text{inst}\). For Gaussian profiles, use \(\beta^2 = \beta_\text{obs}^2 - \beta_\text{inst}^2\). If your peak profiles are pseudo-Voigt, the Gaussian deconvolution is more commonly used.
Calculate the plotting variables
For each peak: compute \(\theta = 2\theta_0 / 2\) in radians, then \(x = 4\sin\theta\) and \(y = \beta\cos\theta\).
Perform linear regression
Fit \(y = b + mx\) using least squares. The y-intercept gives \(D = K\lambda/b\) and the slope gives \(\varepsilon = m\).
Assess the quality of the fit
Check \(R^2\) and the scatter of points around the line. Poor linearity may indicate anisotropic broadening (different crystallite sizes or strain along different crystallographic directions), which the isotropic W-H model cannot capture.
4. Python Implementation
Python — Williamson-Hall plot and analysisimport matplotlib.pyplot as plt
from scipy import stats
# ── Input: peak centers (2θ in degrees) and FWHM (degrees)
# These come from your pseudo-Voigt fitting results
two_theta_peaks = np.array([22.85, 32.52, 40.10, 46.60,
52.48, 57.98, 68.10])
fwhm_deg = np.array([0.182, 0.195, 0.210, 0.228,
0.243, 0.261, 0.295])
# ── Instrumental broadening from standard (LaB6)
# Interpolate to your peak positions if needed
fwhm_inst_deg = np.array([0.060, 0.062, 0.065, 0.068,
0.072, 0.076, 0.085])
# ── Convert to radians
beta_obs = fwhm_deg * np.pi / 180
beta_inst = fwhm_inst_deg * np.pi / 180
# ── Instrumental correction (Gaussian deconvolution)
beta = np.sqrt(beta_obs**2 - beta_inst**2)
# ── Calculate W-H variables
theta = np.radians(two_theta_peaks / 2) # θ in radians
x = 4 * np.sin(theta) # 4sinθ
y = beta * np.cos(theta) # βcosθ
# ── Linear regression
slope, intercept, r_value, p_value, std_err = stats.linregress(x, y)
# ── Extract physical quantities
K = 0.9 # Scherrer shape factor
lam = 1.5406e-10 # Cu Kα in meters
D = K * lam / intercept # crystallite size in meters
D_nm = D * 1e9 # convert to nm
epsilon = slope # microstrain (dimensionless)
print(f'Crystallite size D = {D_nm:.1f} nm')
print(f'Microstrain ε = {epsilon:.4f}')
print(f'R² = {r_value**2:.4f}')
# Data points
ax.plot(x, y, 'ko', markersize=8, label='Experimental')
# Linear fit
x_fit = np.linspace(0, x.max() * 1.1, 100)
y_fit = slope * x_fit + intercept
ax.plot(x_fit, y_fit, 'r--', lw=1.5,
label=f'Linear fit (R² = {r_value**2:.3f})')
# Annotate
ax.annotate(f'D = {D_nm:.1f} nm\nε = {epsilon:.4f}',
xy=(0.05, 0.95), xycoords='axes fraction',
fontsize=11, va='top',
bbox=dict(boxstyle='round', fc='wheat', alpha=0.5))
ax.set_xlabel(r'$4\sin\theta$', fontsize=13)
ax.set_ylabel(r'$\beta\cos\theta$ (rad)', fontsize=13)
ax.set_title('Williamson-Hall Plot')
ax.legend(fontsize=10)
ax.set_xlim(left=0)
ax.set_ylim(bottom=0)
plt.tight_layout()
plt.savefig('williamson_hall.png', dpi=200)
plt.show()
5. Scherrer vs. Williamson-Hall: When Does It Matter?
| Criterion | Scherrer | Williamson-Hall |
|---|---|---|
| Peaks used | Single peak | All peaks simultaneously |
| Strain included? | No — assumed zero | Yes — separated from size |
| Anisotropy? | Per-peak values possible | Isotropic assumption (UDM) |
| When adequate | Well-annealed, strain-free | Nanocrystalline, strained samples |
| Reliability | Overestimates broadening → underestimates size | More accurate if assumptions hold |
A quick diagnostic: if you apply Scherrer to each peak individually and the resulting "sizes" vary systematically with \(2\theta\) (increasing at higher angles), then strain is contributing significantly and you need the W-H method. If all peaks give roughly the same size, strain is negligible and Scherrer is sufficient.
6. Modified Williamson-Hall Methods
The standard W-H equation (Section 2) is called the Uniform Deformation Model (UDM) because it assumes the strain is uniform in all crystallographic directions. In real materials — especially anisotropic structures like orthorhombic perovskites — this assumption often fails, producing poor \(R^2\) values and scattered W-H plots. Several modifications address this:
Uniform Stress Deformation Model (USDM)
If the lattice strain follows Hooke's law, the stress \(\sigma\) and strain \(\varepsilon\) are related by the Young's modulus \(E_{hkl}\) along each crystallographic direction: \(\sigma = E_{hkl}\varepsilon_{hkl}\). The modified equation becomes:
Here the x-axis becomes \(4\sin\theta / E_{hkl}\) and the slope gives the stress \(\sigma\) directly. This requires knowing the elastic constants of your material. For perovskite manganites, these can be estimated from first-principles calculations or taken from literature for the parent compound (e.g., LaMnO3).
Uniform Deformation Energy Density Model (UDEDM)
This model uses the strain energy density \(u\) rather than stress or strain:
The advantage of UDEDM is that it does not assume a particular stress state. The x-axis is \(4\sin\theta / \sqrt{E_{hkl}}\) and the slope squared divided by 2 gives the energy density \(u\).
Start with UDM. If \(R^2 < 0.9\) and the points scatter without a clear trend, try USDM or UDEDM — the anisotropy correction often dramatically improves the linearity. For cubic or nearly cubic structures, UDM is usually sufficient. For orthorhombic Nd0.7Sr0.3MnO3, we have found that the scatter in UDM plots is noticeable, and USDM gives more consistent results across different synthesis conditions.
7. Interpreting the Results
A few rules of thumb for interpreting W-H output:
- Positive slope: Tensile microstrain. The lattice is stretched on average. Common in nanoparticles due to surface tension effects.
- Negative slope: Compressive microstrain. Less common but possible in thin films under substrate-induced stress or in core-shell nanostructures.
- Near-zero slope: Negligible strain. Scherrer alone would have been fine.
- Negative intercept: Unphysical — indicates a problem with your data (wrong instrumental correction, insufficient number of peaks, or dominant anisotropy). Do not report a negative intercept as a valid result.
- Large scatter (low \(R^2\)): Either your FWHM values have large uncertainties (weak peaks, poor fits) or the isotropic strain assumption is violated. Try the modified models or report results with appropriate caveats.
8. Comparison with Rietveld Refinement
Rietveld refinement extracts crystallite size and strain from full-pattern fitting rather than individual peak widths. In FullProf or GSAS-II, size and strain parameters enter through the Thompson-Cox-Hastings profile function, and the refinement simultaneously handles overlapping peaks, preferred orientation, and asymmetry. The W-H method is simpler, more transparent, and requires no structural model — but Rietveld results are generally more reliable because they use all the information in the pattern, not just peak widths.
In practice, I use W-H as a quick sanity check before investing time in a full Rietveld refinement. If the W-H size and strain agree with the Rietveld values (within 10–20%), I am confident that the peak broadening analysis is consistent. If they disagree significantly, it usually points to anisotropic broadening that the isotropic W-H model cannot capture but Rietveld handles through anisotropic strain parameters.
9. Common Pitfalls
A linear fit with 3 data points is not meaningful. You need at least 5–6 peaks spanning a wide \(2\theta\) range to produce a reliable W-H plot. If your phase has fewer visible peaks, the W-H method may not be applicable — fall back on Scherrer with an explicit caveat about strain.
If your sample contains an impurity phase, do not include its peaks in the W-H plot for the main phase. Each phase has its own size and strain, and mixing them will produce nonsensical results. Identify and exclude impurity peaks before plotting.
The W-H method is even more sensitive to instrumental broadening than Scherrer, because the strain information lives in the variation of FWHM with angle. If you do not subtract the instrumental contribution, the trend of \(\beta\cos\theta\) vs. \(4\sin\theta\) will be biased by the angle-dependent instrumental function, and both \(D\) and \(\varepsilon\) will be wrong.
If your W-H plot shows systematic deviations (e.g., some families of reflections above the line and others below), this is evidence of anisotropic strain or anisotropic crystallite shape — not noise. Do not force a line through the data and report the slope and intercept as if the isotropic model is valid. Either switch to the modified models (USDM/UDEDM) or report the anisotropy explicitly.
The Williamson-Hall equation as written above uses the FWHM (\(\beta\)). Some formulations use the integral breadth \(\beta_\text{int}\) (= area/height) instead, which gives slightly different numerical results. Make sure you are consistent: if your fitting code outputs integral breadth, the Scherrer constant \(K\) should be adjusted accordingly (typically \(K = 1.0\) for integral breadth vs. \(K = 0.9\) for FWHM).
10. Complete Pipeline: From FWHM to Size and Strain
Python — Full Williamson-Hall analysis pipelineimport matplotlib.pyplot as plt
from scipy import stats
def williamson_hall(two_theta_deg, fwhm_deg,
fwhm_inst_deg=None,
wavelength=1.5406e-10,
K=0.9,
deconv='gaussian'):
"""
Williamson-Hall analysis (UDM).
Parameters
----------
two_theta_deg : array, peak positions in degrees
fwhm_deg : array, observed FWHM in degrees
fwhm_inst_deg : array or None, instrumental FWHM
wavelength : float, X-ray wavelength in meters
K : float, Scherrer shape factor
deconv : str, 'gaussian' or 'lorentzian'
Returns
-------
dict with D_nm, epsilon, R2, and plot data
"""
# Convert to radians
beta_obs = np.radians(fwhm_deg)
theta = np.radians(two_theta_deg / 2)
# Instrumental correction
if fwhm_inst_deg is not None:
beta_inst = np.radians(fwhm_inst_deg)
if deconv == 'gaussian':
beta = np.sqrt(beta_obs**2 - beta_inst**2)
else: # lorentzian
beta = beta_obs - beta_inst
else:
beta = beta_obs
# W-H variables
x = 4 * np.sin(theta)
y = beta * np.cos(theta)
# Linear regression
res = stats.linregress(x, y)
# Extract quantities
D = K * wavelength / res.intercept
D_nm = D * 1e9
epsilon = res.slope
return {
'D_nm': D_nm,
'epsilon': epsilon,
'R2': res.rvalue**2,
'intercept': res.intercept,
'slope': res.slope,
'slope_err': res.stderr,
'intercept_err': res.intercept_stderr,
'x': x, 'y': y
}
# ── Example usage ──
two_theta = np.array([22.85, 32.52, 40.10, 46.60,
52.48, 57.98, 68.10])
fwhm = np.array([0.182, 0.195, 0.210, 0.228,
0.243, 0.261, 0.295])
fwhm_inst = np.array([0.060, 0.062, 0.065, 0.068,
0.072, 0.076, 0.085])
wh = williamson_hall(two_theta, fwhm, fwhm_inst)
print(f'Crystallite size: {wh["D_nm"]:.1f} nm')
print(f'Microstrain: {wh["epsilon"]:.5f}')
print(f'R²: {wh["R2"]:.4f}')
# ── Also run Scherrer for comparison ──
theta_rad = np.radians(two_theta / 2)
beta_corr = np.sqrt(np.radians(fwhm)**2 - np.radians(fwhm_inst)**2)
D_scherrer = 0.9 * 1.5406e-10 / (beta_corr * np.cos(theta_rad)) * 1e9
print(f'\nScherrer sizes per peak:')
for tt, d in zip(two_theta, D_scherrer):
print(f' 2θ = {tt:.1f}° → D = {d:.1f} nm')
print(f' Mean Scherrer: {D_scherrer.mean():.1f} ± {D_scherrer.std():.1f} nm')
11. Summary: The XRD Processing Pipeline
Across these three posts, we have built a complete workflow from raw diffractometer output to quantitative microstructural parameters:
- Peak Finding — load data, subtract background, smooth, detect peaks, convert to d-spacings.
- Peak Fitting — fit Gaussian/Lorentzian/pseudo-Voigt profiles to extract precise positions, FWHM, and integrated intensities.
- Williamson-Hall — plot \(\beta\cos\theta\) vs. \(4\sin\theta\) to separate crystallite size from microstrain.
Each step feeds into the next, and errors propagate forward. A missed peak in step 1 means a missing data point in the W-H plot. A poorly fitted FWHM in step 2 means a scattered W-H plot with low \(R^2\). Getting the early steps right is not just good practice — it is the foundation that makes the final analysis meaningful.
