Theoretical Background · Dynamical XRD 2θ–ω

READING MAP

What this guide will help you build

This is a detailed theory-to-code pathway for high-resolution coupled 2θ–ω X-ray diffraction scans of layered crystals. Unlike a textbook overview that ends at Bragg’s law, we follow complex amplitudes all the way from material data to a detector-ready intensity.

Atomic physics

CIF geometry, structure factors, anomalous dispersion, susceptibility and the diffraction vector.

Wave propagation

Takagi–Taupin coupling, normalized Riccati equation, branch choice and transfer matrices.

III

Layer stacks

Substrate boundary conditions, reverse layer recursion, strain and coherent attenuation.

Measured scans

Angular convolution, numerical resolution, reproducible tests and physically meaningful fitting.

Assumptions used throughout. Plane-wave illumination, elastic scattering, a near-Bragg two-wave approximation and material properties that vary along depth z but not laterally. We distinguish the model's normalized coordinate T from physical depth z; no formula should be moved between conventions without explicitly transforming it.

Core notation and units

X-ray wavelength; Å for crystal data.

θ, 2θ

Sample angle and detector angle; convert degrees → radians before trig functions.

H = (hkl)

Reciprocal lattice vector; H is not simply a Cartesian direction in non-cubic crystals.

χ₀, χ_±H

Complex Fourier components of electric susceptibility.

D₀, D_H

Complex transmitted and diffracted wave envelopes.

η, T, X

Normalized detuning, optical thickness and amplitude ratio within a declared convention.

γ₀, γ_H

Projections of propagation directions on the surface normal.

Coherent amplitude reduction; distinguish microscopic Debye–Waller from phenomenological damage models.

01 · THE KINEMATIC STARTING POINT

From scattering vectors to Bragg diffraction

An X-ray plane wave produces an elastically scattered wave whose phase depends on the scattering position. The momentum-transfer vector is the difference between the scattered and incident wave vectors; for elastic scattering both have magnitude k = 2π/λ.

$Equation: \mathbf{q}=\mathbf{k}_{H}-\mathbf{k}_{0},\qquad |\mathbf{q}|=\frac{4\pi\sin\theta}{\lambda}$

(1)

For atoms at positions r_j within unit cells translated by R_n, the coherent scattering amplitude factors into a lattice term and an atomic basis term:

$Equation: A(\mathbf{q})=\sum_n e^{i\mathbf{q}\cdot\mathbf{R}_n}\sum_j f_j(\mathbf{q})e^{i\mathbf{q}\cdot\mathbf{r}_j}$

(2)

Constructive interference occurs when the scattering vector meets a reciprocal lattice point, q = H = ha* + kb* + lc*. In a symmetric scan, |q| = 4π sin θ / λ and |H| = 2π/d_hkl; equating them gives Bragg’s law.

$Equation: 2d_{hkl}\sin\theta_B=m\lambda$

(3)

Geometry of Bragg diffraction: Two coherent beams interact with a crystalline film stack above a semi-infinite substrate. Geometrical schematic; not to scale. — **Geometry of Bragg diffraction.** Two coherent beams interact with a crystalline film stack above a semi-infinite substrate. Geometrical schematic; not to scale.

Why a finite film produces interference fringes

A uniform, sufficiently thin film of thickness t accumulates an extra phase Δq_zz away from the exact reciprocal-lattice condition. Integrating over depth gives the familiar thickness-amplitude formula, and the squared modulus produces the sinc-shaped envelope.

$Equation: A_f(\Delta q_z)\propto\frac{e^{i\Delta q_z t}-1}{i\Delta q_z}$

(4)

$Equation: I_f(\Delta q_z)\propto t^2\left[\frac{\sin(\Delta q_z t/2)}{\Delta q_z t/2}\right]^2$

(5)

Kinematic thickness fringes: Numerically evaluated sinc-squared envelopes for relative film thicknesses 1, 2 and 3. Calculated curves, not experimental scans. — **Kinematic thickness fringes.** Numerically evaluated sinc-squared envelopes for relative film thicknesses 1, 2 and 3. Calculated curves, not experimental scans.

The plot is useful for intuition: thicker films give narrower main peaks and more closely spaced fringes. But it ignores multiple scattering, extinction, refraction and strong substrate reflections. Near an intense Bragg reflection, adding separate intensities for layers is fundamentally wrong; the correct object to propagate is the complex amplitude.

Implementation insight. As a rough angular-grid estimate, differentiating q_z = 4π sin θ/λ gives:

$Equation: \Delta(2\theta)_{\rm fringe}\simeq\frac{\lambda}{t\cos\theta_B}\quad(\mathrm{radians})$

(6)

This is a kinematic estimate in radians for the 2θ axis, not a substitute for numerical convergence tests.

02 · MATERIALS AND CRYSTALLOGRAPHY

Structure factors become a complex susceptibility

Every reflection depends on the atom positions, occupancies, wavelength and atomic form factors. A useful crystallographic implementation separates these inputs from thickness and other sample-stack parameters. The scattering factor includes non-resonant and anomalous terms:

$Equation: f_j(s,\lambda)=f_j^0(s)+f_j\prime(\lambda)+if_j\prime\prime(\lambda),\quad s=\sin\theta/\lambda$

(7)

The crystal structure factor for reciprocal vector H = (hkl) is a coherent sum over atoms in the unit cell. Occupancy o_j, Cartesian-equivalent fractional coordinates and atomic displacement factors all enter before the dynamical solver:

$Equation: F_H=\sum_j o_j f_j\exp[2\pi i(hx_j+ky_j+lz_j)]\,e^{-W_j(H)}$

(8)

Here W_j is the crystallographic displacement exponent. For a compound or alloy, occupancies must follow the modelled composition. For a monoclinic or hexagonal crystal, use the full unit-cell metric to compute interplanar distances and the scattering direction; do not use cubic-only shortcuts.

The connection to Maxwell’s wave equation is made through the Fourier components of the complex electric susceptibility:

$Equation: \chi_H=-\frac{r_e\lambda^2}{\pi V_c}F_H,\quad \chi_0=-\frac{r_e\lambda^2}{\pi V_c}F_{000}$

(9)

Here r_e is the classical electron radius, V_c the unit-cell volume and F₀₀₀ forward scattering. Under the convention shown, χ₀ produces refraction and absorption; χ_H and χ_−H provide coherent Bragg coupling.

Important with anomalous absorption: do not blindly substitute F_−H = F_H*. The relation generally fails when the atomic scattering factors are complex. Likewise, a change of Fourier sign or temporal wave convention changes the sign of the displayed susceptibility formula.

Cache the right quantities

For many trial thicknesses or strain profiles of the same material, the undistorted atomic positions and scattering tables can be reused. However, d_hkl, Bragg geometry, wavelength-dependent factors and any angle-dependent scattering approximation must be invalidated when their physical inputs change. The MROX Studio 41 manual distinguishes full and fast structure-factor evaluation; the latter must still be keyed by the relevant per-layer Bragg condition.

03 · DYNAMICAL THEORY

Deriving the Takagi–Taupin equations

When a wave repeatedly scatters between the transmitted and diffracted channels, the fields are coupled. Begin with the scalar Helmholtz equation for one selected polarization (a Maxwell/vector treatment is required when the scalar approximation fails):

$Equation: \nabla^2\Psi+k^2[1+\chi(\mathbf{r})]\Psi=0,\qquad k=2\pi/\lambda$

(10)

Close to a single reciprocal-lattice vector H, retain the average and two resonant harmonics of the susceptibility. Distortion is represented by a displacement field u(r) which modulates the lattice phase:

$Equation: \chi(\mathbf{r})\simeq\chi_0+\chi_H e^{i\mathbf{H}\cdot(\mathbf{r}-\mathbf{u})}+\chi_{-H}e^{-i\mathbf{H}\cdot(\mathbf{r}-\mathbf{u})}$

(11)

Decompose the internal field into two nearly resonant carriers with slowly varying envelopes, keeping the transmitted and Bragg-diffracted beams explicitly:

$Equation: \Psi(\mathbf{r})=D_0(\mathbf{r})e^{i\mathbf{k}_0\cdot\mathbf{r}}+D_H(\mathbf{r})e^{i(\mathbf{k}_0+\mathbf{H})\cdot\mathbf{r}}$

(12)

Derivation in four steps

Insert the two-wave expression into the Helmholtz equation and collect coefficients of the two resonant exponential carriers.
Discard terms oscillating at non-resonant reciprocal vectors (two-wave approximation).
Apply the slowly varying envelope approximation: neglect ∇²D relative to 2ik·∇D.
Project the remaining directional derivatives along transmitted and diffracted ray coordinates s₀ and s_H.

In a generic sign convention, the first-order coupled system has the structure:

$Equation: \partial_{s_0}D_0=ia_0D_0+ib_0D_H e^{-i\mathbf{H}\cdot\mathbf{u}}$

(13)

$Equation: \partial_{s_H}D_H=ia_HD_H+ib_HD_0 e^{+i\mathbf{H}\cdot\mathbf{u}}$

(14)

The coefficients a₀, a_H contain refraction and deviation from exact Bragg diffraction, while b₀, b_H contain the Fourier susceptibilities, polarization factor and geometry. The exact assignment of χ_H versus χ_−H and displacement-phase signs follows the chosen Fourier convention: derive both equations consistently rather than mixing sources.

For a laterally uniform stratified crystal, all material properties depend on depth only. The coupled PDEs then reduce to ordinary differential equations. This 1D reduction enables exact constant-slab propagation and is the relevant model for the present 2θ–ω page.

When not to use this page’s 1D solver. Dislocations with strong lateral gradients, diffraction imaging, bending with in-plane spatial dependence, broad angular divergence or multiple simultaneously excited reciprocal vectors may require full 2D/3D TT or many-beam dynamical diffraction.

04 · A SINGLE COMPLEX UNKNOWN

Why the amplitude ratio satisfies Riccati

A consistent change of coordinate, phase and amplitude normalization transforms the two-wave system into a simple canonical pair. This is an explicitly declared mathematical convention, not a claim that the legacy MROX variable Xt has exactly these signs.

$Equation: -i\frac{dU}{dT}=-V,\qquad -i\frac{dV}{dT}=U+2\eta V$

(15)

Introduce the normalized ratio X(T) = V(T)/U(T). Differentiate this quotient:

$Equation: X=\frac{V}{U},\qquad \frac{dX}{dT}=\frac{UV\prime-VU\prime}{U^2}$

(16)

Using U′ = −iV and V′ = iU + 2iηV, we obtain −iX′ = (U² + 2ηUV + V²)/U². Hence:

$Equation: -i\frac{dX}{dT}=X^2+2\eta X+1$

(17)

Notice that absolute phase and intensity scale have disappeared, but coherent interference is retained because X remains complex. This reduction is the mathematical basis for efficient amplitude-ratio recursions in 1D Bragg diffraction.

Linearisation provides an independent check

From the canonical relation U′ = −iV, we have X = iU′/U. Substitution into the Riccati equation cancels the quadratic derivative terms and produces a second-order linear equation:

$Equation: X=i\,U\prime/U\quad\Rightarrow\quad U\prime\prime-2i\eta U\prime-U=0$

(18)

Solving this linear equation and forming V/U must agree with a direct Riccati solver and a 2×2 matrix exponential. That three-way equivalence is a valuable unit test.

05 · ANALYTIC SOLUTION

Fixed points, branches and numerical stability

For constant η, the two stationary solutions are the roots of the Riccati quadratic:

$Equation: X_{\pm}=-\eta\pm\sqrt{\eta^2-1},\qquad \xi=\sqrt{\eta^2-1}$

(19)

For real η without absorption, the interval |η| < 1 is the ideal Darwin region: the square root is imaginary and strong coupling occurs. Outside it, the roots are real. In an absorbing crystal, η and the square root are complex; the physically acceptable root is selected by analytic continuation and the non-growing wave condition, not by an arbitrary sign test alone.

$Equation: \frac{X(T)-X_+}{X(T)-X_-}=\frac{X(0)-X_+}{X(0)-X_-}\,e^{2i\xi T}$

(20)

This exact expression propagates a starting ratio through a homogeneous slab. It is compact but not always the most stable numerical representation: close to a root, both numerator and denominator may become small, and exponential factors can overflow for thick or absorbing layers. Matrix/ratio/scattering formulations should be cross-checked.

Ideal Darwin plateau: Analytic perfect-substrate, lossless, two-beam reflectivity as a function of normalized detuning. The plateau is an ideal limit. — **Ideal Darwin plateau.** Analytic perfect-substrate, lossless, two-beam reflectivity as a function of normalized detuning. The plateau is an ideal limit.

Exact finite-slab response: Calculated perfect slab two-wave reflectivity for four optical thicknesses. No absorption, instrumental convolution or material-specific coefficients. — **Exact finite-slab response.** Calculated perfect slab two-wave reflectivity for four optical thicknesses. No absorption, instrumental convolution or material-specific coefficients.

INTERACTIVE MATHEMATICAL EXAMPLE

Build a rocking curve by changing the optical thickness

Ideal · no absorption

Optical thickness T1.5

Computed from the exact dimensionless two-wave slab transfer matrix with X_in = 0, not from a fitted experimental diffractogram.

06 · EXACT HOMOGENEOUS SLABS

Transfer matrices and the Möbius map

Write the canonical two-wave system as Y′ = A Y where Y = (U, V)^T and A is constant inside a homogeneous slab:

$Equation: \frac{d\mathbf{Y}}{dT}=\mathbf{A}\mathbf{Y},\quad A_{11}=0,\ A_{12}=-i,\ A_{21}=i,\ A_{22}=2i\eta$

(21)

Then Y_out = exp(AΔT)Y_in. The matrix explicitly includes coherent coupling, so no interference phases need to be approximated by hand. Instead of propagating arbitrarily scaled fields, form their ratio. Any 2×2 linear transform induces the fractional linear (Möbius) map:

$Equation: X_{\rm out}=\frac{M_{21}+M_{22}X_{\rm in}}{M_{11}+M_{12}X_{\rm in}}$

(22)

The key computational advantage is that X_out depends only on the starting ratio and the slab parameters, so the whole stack can be traversed without storing absolute fields at every interface. However, a singular or near-zero denominator must be handled carefully, especially for absorbing and very thick slabs.

Verified dimensionless reference · PythonDownload complete example ↓

import numpy as np

def propagate(x_in, eta, T):
    """Exact normalized two-wave transfer; not production MROX physics."""
    eta = complex(eta)
    s = np.sqrt(1.0 - eta*eta)
    h = T if abs(s) < 1e-10 else np.sinh(s*T)/s
    c = np.cosh(s*T)
    # Overall exp(i*eta*T) cancels from the complex ratio.
    m11 = c - 1j*eta*h
    m12 = -1j*h
    m21 = 1j*h
    m22 = c + 1j*eta*h
    return (m21 + m22*x_in)/(m11 + m12*x_in)

# Two-beam, lossless finite slab; dimensionless detuning.
print(abs(propagate(0j, eta=0., T=1.5))**2)

Why does the code use hyperbolic functions? The trace-shifted matrix A − iηI squares to (1−η²)I. Therefore the matrix exponential is a combination of cosh(sT) and sinh(sT)/s, where s = √(1−η²). The s → 0 limit should be handled analytically as sinh(sT)/s → T.

07 · THE LOWER BOUNDARY CONDITION

The semi-infinite substrate and the Darwin plateau

The bottom of a Bragg stack is usually treated as a semi-infinite crystal. One does not prescribe an unrelated diffracted wave arriving from infinite depth. Instead, choose the eigenmode that meets the physically appropriate boundedness condition at depth, including absorption and the direction of energy flow.

$Equation: X_{\infty}=-\eta\pm\sqrt{\eta^2-1}\qquad(\mathrm{normalised\ convention})$

(23)

That expression shows the canonical roots. In the MROX Studio 41 implementation described in the supplied manual, the substrate ratio is written in a different normalization as Xt = η ∓ √(η²−1). The change of sign is not a bug by itself; it is a change of convention. A developer must transform η, T and X together, then verify the predicted observable.

In the ideal nonabsorbing, symmetric and properly flux-normalized case, the substrate intensity shows the Darwin plateau. Within the two-wave framework a common surface intensity expression includes the structure-factor ratio:

$Equation: I_H\propto\left|\frac{F_H}{F_{-H}}\right|\,|X_{\rm surface}|^2$

(24)

Its flat high-reflectivity region is a defining distinction from the thin-crystal kinematical approximation. In real XRD data the plateau is modified by absorption, asymmetric geometry, polarization, instrumental response and crystal perfection.

Do not sum intensities across slabs. When an epitaxial film sits on a substrate, all interface reflections modify the complex ratio at the surface. Squaring each slab's amplitude separately and adding the results destroys the phase-sensitive Pendellösung and thickness-fringe patterns.

08 · THE BARTELS IDEA

Multilayer propagation from substrate to surface

Bartels, Hornstra and Lobeek showed how dynamical diffraction from complex multilayers can be computed by a recursive amplitude-ratio treatment derived from Takagi–Taupin equations. The underlying physical boundary condition is known at the semi-infinite substrate, so the recursion proceeds upward through the films.

Complete two-wave 2θ–ω calculation pipeline: Physical flow from crystallography and scan geometry through complex amplitude recurrence to instrument response and comparison with data. — **Complete two-wave 2θ–ω calculation pipeline.** Physical flow from crystallography and scan geometry through complex amplitude recurrence to instrument response and comparison with data.

For each homogeneous layer j and each angular point, prepare its detuning η_j, dynamical thickness T_j, complex coupling factors and the ratio arriving from the layer below. The supplied MROX Studio 41 manual documents this specific Bartels-family algebra:

$Equation: S_1=(X_{\rm in}-\eta+\xi)e^{-iT\xi},\quad S_2=(X_{\rm in}-\eta-\xi)e^{iT\xi}$

(25)

$Equation: X_{\rm out}=\eta+\xi\frac{S_1+S_2}{S_1-S_2}$

(26)

The dynamical thickness in that documented implementation can be written as

$Equation: T_j=\frac{\pi C_j\gamma_{M,j}\,t_j\sqrt{F_{H,j}F_{-H,j}}}{\lambda\sqrt{|\gamma_{0,j}\gamma_{H,j}|}}$

(27)

where γ_M = r_e λ²/(π V_c); t_j and λ must use identical length units (for example, 1 nm = 10 Å). Polarization enters through C; γ₀ and γ_H contain the beam/surface geometry. A characteristic dynamical coupling length is proportional to:

$Equation: \ell_{\rm dyn}\sim\frac{\lambda\sqrt{|\gamma_0\gamma_H|}}{\pi|C|\sqrt{|\chi_H\chi_{-H}|}}$

(28)

The proportionality and exact “extinction length” definition vary between authors and geometry conventions. The useful interpretation is universal: the larger the coupling and the thicker the crystal relative to its dynamical length scale, the less reliable the single-scattering approximation.

Developer warning — no unverified formula swapping. The Bartels-family expressions above are transcribed as the manual documents them; the canonical Riccati derivation on this page uses a different X/η/T normalization. They have the same two-wave/Möbius mathematical structure, but literal numerical equivalence requires an explicit convention mapping and reference tests. Do not simply replace the recurrence in working code.

A safe stack algorithm

Prepare: read the crystallography, wavelength, chosen hkl, geometry and the 2θ grid.
Calculate: F₀, F_H, F_−H, χ, γ and layer-specific η/T at each angular point.
Boundary: select the semi-infinite substrate eigenmode on the physically correct analytic branch.
Recurse: propagate the complex ratio from the deepest finite slab to the first film.
Observe: convert the surface complex ratio to diffraction intensity with the convention-specific flux/structure normalization.
Instrument: convolve intensity, add correctly modelled spectral components and background only after coherent propagation.

09 · DEFORMED CRYSTALS

How strain moves an XRD reflection

A spatially varying displacement field causes a phase shift H·u(z) in the diffracted wave and changes the local condition for Bragg diffraction. At first order, the relative change in interplanar spacing is the strain tensor projected along the actual reciprocal-space normal:

$Equation: \frac{\Delta d_{hkl}}{d_{hkl}}\simeq\hat{\mathbf{n}}^{T}\boldsymbol{\varepsilon}\hat{\mathbf{n}}$

(29)

Differentiate Bragg’s law at fixed wavelength to obtain the angular shift. When plotting a 2θ scan, the factor of two comes from the detector axis:

$Equation: \Delta(2\theta_B)\simeq-2\tan\theta_B\,\frac{\Delta d_{hkl}}{d_{hkl}}\quad(\mathrm{radians})$

(30)

Bragg shift under normal strain: First-order differentiation of Bragg’s law: positive strain decreases 2θ. Curves are theoretical estimates. — **Bragg shift under normal strain.** First-order differentiation of Bragg’s law: positive strain decreases 2θ. Curves are theoretical estimates.

A positive normal strain increases the plane spacing, shifting the reflection to smaller 2θ; a contraction shifts it upward. The plot shows the first-order rule for several example strain values. At larger strain or near complex material geometries, update the complete lattice metric instead of relying on the derivative.

Elasticity and pseudomorphic layers

The in-plane strain of an epitaxial film can be constrained by the substrate lattice; the out-of-plane strain follows from the elastic tensor and the appropriate free-surface stress conditions. In tensor/Voigt notation:

$Equation: \sigma_i=\sum_j C_{ij}\varepsilon_j,\qquad \sigma_{nn}=0\quad(\mathrm{free\ surface})$

(31)

For cubic (001) layers, textbook biaxial relations are often sufficient, but they do not generalise blindly to monoclinic β-Ga₂O₃, non-cubic cuts or arbitrary off-axis reflections. A robust solver builds the direct and reciprocal lattice bases from the CIF, rotates the elastic stiffness tensor into sample coordinates, then solves the plane-stress conditions.

Units trap: a 0.5% deformation is 0.005 as a fraction. The supplied Studio 41 manual notes that its layer-state interface stores strain in percent, whereas internal mathematical expressions require fractions. Mixing them shifts peaks by a factor of 100.

10 · COHERENCE AND DISORDER

Debye–Waller: atomic meaning versus fit parameter

Thermal motion and static disorder make the atomic scattering phases fluctuate. For a zero-mean Gaussian displacement vector with covariance tensor U, the coherent average scattering amplitude at reciprocal vector H is:

$Equation: \left\langle e^{i\mathbf{H}\cdot\mathbf{u}}\right\rangle=e^{-\frac{1}{2}\mathbf{H}^T\mathbf{U}\mathbf{H}}$

(32)

For an isotropic displacement along H, this reduces to exp(−½|H|²σ_u²). Intensity is reduced by the squared modulus, not by the same factor again unless that is how the chosen parameter is explicitly defined.

Microscopic Debye–Waller damping: Coherent scattering amplitude versus projected root-mean-square atomic displacement at different |H| values, computed from the Gaussian displacement model. — **Microscopic Debye–Waller damping.** Coherent scattering amplitude versus projected root-mean-square atomic displacement at different |H| values, computed from the Gaussian displacement model.

This dependence on |H|² is the key microscopic prediction: higher-order Bragg reflections may be more sensitive to displacement disorder. It is not equivalent to an arbitrary empirical attenuation that happens to be called “DW”. A fitted Debye–Waller-like damage parameter can also absorb mosaicity, inhomogeneity or instrumental effects if the model is insufficiently constrained.

A subtle forward-scattering issue

In standard crystallographic displacement theory the Debye–Waller factor damps nonzero reciprocal harmonics F_H but does not reduce F₀₀₀ by the same mechanism, because H=0 makes the phase fluctuation vanish. The supplied manual flags that the examined MROX Studio 41 implementation multiplies F_H, F_−H and F₀₀₀ by a layer DW factor. It records this as a physics question for testing, not as a validated microscopic identity.

For independent implementations: decide whether DW means microscopic mean-square displacement, an empirical coherent-volume fraction, or a phenomenological implantation damage profile. Give the model a distinct name and document exactly which amplitudes it multiplies.

11 · CONTINUOUS DEPTH PROFILES

Representing ion implantation with thin slabs

Ion implantation changes lattice parameter, static disorder and sometimes the coherent diffracting fraction as a function of depth. A piecewise-constant discretisation approximates continuous profiles by many small slabs, each with its own strain, atomic parameters and coherent attenuation.

The supplied manual describes a sum of asymmetric Gaussian strain components, with a potentially different width on each side of the maximum:

$Equation: \varepsilon(z)=\sum_{k=1}^{K}A_k\exp\left[-\frac{(z-\mu_k)^2}{2\sigma_k(z)^2}\right]$

(33)

$Equation: \sigma_k(z)=\sigma_{L,k}\ [z\leq\mu_k];\qquad \sigma_{R,k}\ [z>\mu_k]$

(34)

SCIENTIFIC FIGURE · ILLUSTRATIVE IMPLANTATION STRAIN DEPTH PROFILECalculated / theoretical

Illustrative implantation strain depth profile. Two asymmetric Gaussian functions and their sum, evaluated as a mathematical example. Not an SRIM measurement or an MROX fit.

The amplitudes A_k in that model are expressed in percent strain; the centres and widths are distances. Convert displayed FWHM values to Gaussian standard deviations using FWHM = 2√(2 ln 2) σ. When discretising, take values consistently at slab midpoints or use slab-average integration if higher accuracy is required.

To maintain physical comparability under refinement, changing the slab step Δz must leave the total affected thickness and continuous profile unchanged. Shorten the final slab if necessary rather than silently rounding the thickness upward. Check the complete stack after applying a profile: the material, strain, thickness and DW of unrelated layers must remain unchanged.

One project-specific operational choice is not a physical law. The manual documents a “selected layer only” application mode which can extend that layer’s effective thickness if the selected implantation profile extends deeper. For real samples, ion ranges do not create additional matter. A new program should decide explicitly how to handle a profile crossing an actual interface, rather than copying this UI behaviour without thinking about sample physics.

12 · FROM IDEAL TO MEASURED

Polarisation, spectral lines and instrumental broadening

The two-wave result describes an ideal coherent diffraction intensity. The experimental detector sees a response affected by source spectrum, divergence, angular acceptance and background. Typically the instrument acts on intensity after coherent wave propagation:

$Equation: I_{\rm measured}(x)=B(x)+A\int I_{\rm dyn}(u)R(x-u)\,du,\quad x=2\theta$

(35)

A widely used parameterisation of an angular response is the pseudo-Voigt function, mixing Gaussian and Lorentzian line shapes. Define both components with the same FWHM w and each normalised to unit integrated area:

$Equation: PV(x)=(1-\beta)G(x;w)+\beta L(x;w),\quad 0\leq\beta\leq1$

(36)

$Equation: G(x;w)=\frac{2\sqrt{\ln2}}{w\sqrt{\pi}}\exp\left(-4\ln2\frac{x^2}{w^2}\right)$

(37)

$Equation: L(x;w)=\frac{2}{\pi w}\frac{1}{1+4x^2/w^2}$

(38)

Gaussian, Lorentzian and pseudo-Voigt: Illustrative peak-normalised line shapes, showing how their tails differ. The displayed example uses consistent shape widths. — **Gaussian, Lorentzian and pseudo-Voigt.** Illustrative peak-normalised line shapes, showing how their tails differ. The displayed example uses consistent shape widths.

Why the instrument changes a diffraction curve: Synthetic ideal peak family before and after a Gaussian intensity convolution; this is a signal-processing example. — **Why the instrument changes a diffraction curve.** Synthetic ideal peak family before and after a Gaussian intensity convolution; this is a signal-processing example.

Within a single two-wave scalar approximation, a common polarisation coupling factor is C_σ = 1 or C_π = cos(2θ_B). Some intensity-only codes use its absolute value inside their coupling normalisation; if this is done, keep phase conventions consistent throughout propagation.

What changes for a Kα doublet?

The Cu Kα₁ and Kα₂ photons have different wavelengths. Therefore, their peak positions must be calculated from their respective wavelengths via Bragg’s law, and the intensities added incoherently according to the source spectrum and optical selection. Reusing a single simulated profile and applying an arbitrary angular translation is only an approximation whose accuracy depends on scan range and material.

Beware double counting. If complex χ₀ already contains absorption, applying an additional Beer–Lambert attenuation to the final reflected intensity without a separate physical reason can count absorption twice.

13 · ANGLE-DEPENDENT FILTERS

Resolution that changes across the scan

The supplied manual separately describes a legacy-inspired empirical model in which Gaussian smoothing increases with angular distance from the substrate Bragg position. It is a useful fit model but is not a first-principles model of every diffractometer:

$Equation: \sigma(x)=\sigma_{\min}+(\sigma_{\max}-\sigma_{\min})\left(\frac{|x-x_B|}{D_{\max}}\right)^p$

(39)

Here x = 2θ, D_max = max_x|x−x_B|, and the width interpolates between σ_min and σ_max with exponent p > 0. Because each output point has its own kernel width, ordinary stationary convolution is not exactly the same operation.

$Equation: I_{\rm res}(x_i)=\frac{\sum_jI(x_j)\exp[-(x_j-x_i)^2/(2\sigma_i^2)]}{\sum_j\exp[-(x_j-x_i)^2/(2\sigma_i^2)]}$

(40)

Empirical variable-width Gaussian resolution: Example width law that grows with distance from the substrate peak, illustrating the scan-window dependence of its parameters. — **Empirical variable-width Gaussian resolution.** Example width law that grows with distance from the substrate peak, illustrating the scan-window dependence of its parameters.

The displayed discrete version is a locally normalised Gaussian weighted average; it preserves a constant signal. For uneven angular grids, sampling density must be accounted for with integration/quadrature weights if one wants a proper continuous integral approximation. For a scan with very small σ compared with the angular step, oversampling is needed before claiming accurate sub-step broadening.

Window dependence: D_max is computed from the current scan limits. Cropping the same experimental scan therefore changes the effective width profile if the parameters are held fixed. Never compare fitted resolution parameters across different angular windows without reporting that definition.

14 · A PRACTICAL DEVELOPER ROADMAP

Implement a robust 2θ–ω solver

For someone developing a new program, it helps to separate a minimal mathematically verifiable solver from a fully physical crystal simulator. The first proves that the amplitude algebra is consistent; the second adds material, geometry, boundary conditions and instrument physics. Implement these layers in this order.

Declare the conventions before the code

Write down the time dependence (e^−iωt or e^+iωt), reciprocal-space Fourier sign, direction of depth, orientation of surface normal, transmitted/diffracted amplitude normalisation and branch criteria. Fix a single meaning for η, T and X.

Implement and test the dimensionless two-wave kernel

Implement Y′ = AY and the Möbius map first. Compare matrix propagation with the analytic Riccati formula, including T→0, η=±1 and long slabs. The included riccati_reference.py performs some of these tests.

Introduce real crystallography

Read CIF unit-cell parameters, symmetry and fractional coordinates; form the reciprocal metric, d_hkl, V_c and all relevant complex F. Verify the results against known reference reflections before entering the TT code.

Map material quantities to propagation

Calculate susceptibilities and the geometry-dependent coefficients at each angular point. Transform into the exact normalized variables of the chosen recurrence. Establish and record how the flux factor converts surface X into measured reflectivity.

Choose the substrate branch physically

Use radiation/decay conditions and analytic continuation under absorption. Test continuity as detuning passes through the Darwin edges and compare against independently published perfect-crystal results.

Add layers, strain and damage

Recurse bottom-to-top, preserving complex phases. Update the full lattice metric under strain; discretise profiles by depth while retaining physical thickness. Separate microscopic Debye–Waller from empirical coherent attenuation.

Convolve only after coherent propagation

Apply source spectrum, polarisation, calibrated resolution, background and any explicit angular-zero correction. Preserve the original experimental angular grid; do not move measured points to improve agreement.

Keep fitting outside the forward model

Use a pure function simulate(parameters, angles) that does not mutate UI state or files. Cache by dependencies, not by convenient but incomplete labels. Let fitting algorithms, uncertainties and optional Bayesian inference call the same forward engine.

Reference software and regeneration scripts

01Riccati / transfer matrixFully disclosed dimensionless Python test implementation↓02All 12 scientific figuresPython, NumPy, SciPy and Matplotlib; generates SVGs↓03Typeset equationsRebuilds 44 mathematical SVG formulae without CDN↓

All visual assets are supplied locally. The mathematical demonstration code is open for study and modification; it is not a reconstructed MROX production solver and cannot predict a specific material before susceptibility, geometry and boundary conditions are implemented.

15 · VALIDATION AND NUMERICAL AUDIT

Tests that separate correct physics from plausible plots

A simulated rocking curve can look beautiful while being physically wrong. The following tests should be automatic and run when any sign convention, absorption formula or propagation kernel changes.

Test	Expected property	What it catches
Zero slab thickness	X_out → X_in as T → 0.	Incorrect transfer order, signs and singular limits.
Analytic Riccati vs matrix	Same complex amplitude in the same declared convention.	Formula transpositions and eigenvalue errors.
Ideal finite slab	Finite, nonnegative intensity; flux conservation in lossless normalised cases.	Wrong intensity factor or unstable propagation.
Substrate branch	Physically decaying/causal solution with smooth absorption limit.	Wrong square-root sheet and unphysical gain.
Identical adjacent layers	Splitting a homogeneous slab must not change reflectivity.	Incorrect interface and matrix composition.
Zero diffraction coupling	Bragg-reflected amplitude tends to zero appropriately.	Spurious kinematic/dynamical mixing.
Angular grid refinement	Curves and peak positions converge under halving angular spacing.	Aliased fringe periods and false peak maxima.
Depth-grid refinement	Profiled-stack reflectivity converges as Δz is halved.	Under-resolved implantation gradients.
Instrument OFF	Raw coherent intensity unchanged.	Hidden filtering and double convolution.
Constant-input filter	Local normalised Gaussian returns the same constant.	Bad boundary normalisation.
Absorption and anomalies	No unphysical amplification; complex F_−H handled independently.	Conjugation assumptions and branch errors.
Strain sign and units	Expanded spacing shifts the ideal peak towards lower 2θ.	Percent/fraction errors and flipped geometry.

A useful slab-grid convergence metric is:

$Equation: E_{\Delta z}=\frac{\max_i|I_{\Delta z}(x_i)-I_{\Delta z/2}(x_i)|}{\max(\epsilon,\max_iI_{\Delta z/2}(x_i))}$

(41)

Here the two curves are evaluated on the same angular points and with the same underlying continuous strain/DW profile. Relative to full-scale maximum intensity, E_Δz is a reproducible discretisation estimate, though weak features may need a log-intensity or local criterion.

Numerical hazards worth planning for

Square-root branch

Complex square roots change sign across branch cuts. Track the physical eigenmode continuously, with an explicit absorption sign convention.

Large exponentials

Transfer matrices can overflow or lose significant digits for thick layers. Consider normalised Möbius maps and scattering matrices.

Near-zero denominators

The fractional transform can become ill-conditioned even when the physical reflectivity remains finite. Add tests and safe alternative forms.

Finite angular sampling

Sub-step fringes are not accurately recovered by a smoothing filter. Oversample and integrate appropriately.

Fitting is not physical validation

As an example, the supplied manual records a deterministic objective combining robust logarithmic error, overall intensity correlation and first-difference correlation:

$Equation: F=0.60L_{\rm Huber}+0.25(1-\rho_I)+0.15(1-\rho_{\Delta I})$

(42)

$Equation: H_{\delta}(e)=e^2/2\ [|e|\leq\delta];\quad \delta(|e|-\delta/2)\ [|e|>\delta]$

(43)

This is a project-specific optimisation score, not a normalised probabilistic likelihood. If a new solver reports Bayesian credible intervals, it must instead specify a statistical residual model and prior. For independent Gaussian log-intensity errors with common scale σ, the log-likelihood has the generic form:

$Equation: \log p(D|\boldsymbol{\vartheta},\sigma)=-N\log\sigma-\frac{1}{2\sigma^2}\sum_i r_i^2+\mathrm{constant}$

(44)

Instrumental width, strain, roughness/damage and thickness can be highly correlated when fitting only one symmetric reflection. Confidence in extracted physical parameters is improved by independent measurements, more reflections, explicit uncertainty analysis and robustness across initial guesses.

16 · WHAT THIS REFERENCE CLAIMS

Boundaries of the theory and source

What is mathematical and reproducible

Bragg/kinematic relations, susceptibility mapping under the stated sign convention, the dimensionless Riccati identity, constant-slab matrix propagation and the figures computed from those equations. The supplied reference Python script checks the analytic/matrix identity numerically.

What requires physical validation

Actual MROX Studio 41 material-specific numerical outputs, branch mapping between its Xt and canonical X, the phenomenological DW treatment of F₀₀₀, resolution calibration and any inferred sample-specific strain or disorder.

The supplied manual explicitly describes a static code audit rather than experimental or MATLAB runtime validation. The present page follows that distinction. It explains scientific methods openly and teaches reproducible mathematics; it does not publish unverifiable experimental results or claim identity with proprietary algorithms solely on the basis of similar-looking expressions.

SELECTED SCIENTIFIC SOURCES

References and further reading

Foundational references for independent readers who want to inspect derivations and build their own software. DOI links lead to the publishers’ original records.

W. J. Bartels, J. Hornstra & D. J. W. Lobeek (1986). “X-ray diffraction of multilayers and superlattices.” Acta Crystallographica A 42, 539–545. doi:10.1107/S0108767386098768 ↗. Central multilayer recurrence.
J.-P. Guigay & M. Sanchez del Rio (2024). “Formulation of perfect-crystal diffraction from Takagi–Taupin equations: numerical implementation in the crystalpy library.” Journal of Synchrotron Radiation 31, 1469–1480. doi:10.1107/S160057752400924X ↗. Independent matrix/scattering formulation.
M. Carlsen & H. Simons (2022). “A finite difference scheme for integrating the Takagi–Taupin equations on an arbitrary orthogonal grid.” Acta Crystallographica A 78, 395–401. doi:10.1107/S2053273322004934 ↗. Beyond 1D slab geometry.
IUCr Commission on Crystallographic Nomenclature. Atomic displacement parameter nomenclature ↗. Crystallographic definitions of thermal/disorder parameters.
MROX Studio 41 — 2θ–ω scientific manual (3 November 2026; author-provided source). Principal source for the MROX-specific implementation notes, the documented Bartels-family recurrence and its declared validation limitations.

CONTINUE EXPLORING

From equations to research applications

Explore how XRD/XRR modelling is used in real scientific publications, or study the reproducible example code included with this page.

MROX Publications ↗Reference Python code ↓

Dynamical diffraction,from first principles