Deconvolution of distribution of diffraction reflexes

Question

I'm a chemist stuck in a mathematical problem. Please bear with me as I'm trying to express myself in Math language.

Let me explain in short terms the experimental method I'm using: X-ray diffraction. For this method, atoms in a solid can be viewed as ordered balls in space. X-rays, having similar wavelength as the distance between atoms, are diffracted by those balls and can be detected as diffraction fringes. The distance between the diffracting atoms (or planes of atoms, aka lattice distance) can then be calculated using Bragg's law We use the method to measure the symmetry of the atom ordering, their distances and placement in space.

As chemists, we like to deduce so-called lattice parameters (the lengths of the smallest repeating unit of atoms in all three cartesian coordinates) from three or more characteristic atomic plane distances. In my case the repeating unit has the form of a cuboid. The measured diffraction fringes are equivalent to the points of the fourier transform of the spacing of the atoms in real space and describe the repetition of those atoms. They are linked to the cuboid edge lengths via:

1/d² = h²/a² + k²/b² + l²/c²

h, k and l are discrete values, characteristic to the specific diffraction fringe, d is the measured lattice distance.

Unfortunately, the cuboid is not stiff, but a, b and c have some distribution. As such I have measured an asymmetric distribution of the lattice distances d. I would like to find out how the cuboid edge lengths contribute to the measured distribution.

I fitted a normal distribution to my experimental data for 1/d², for which it is easy to deconvolute to get distributions of 1/a² and so on. This doesn't reflect my data correctly, though, since I observe an asymmetric distribution of 1/d².

Do you have an idea how to extract information of the distribution of a, b and c from those three distributions of lattice spacings? In a first approximation, these three edge lengths are independent variables.

Answer 1

What helps is to re-express everything in variables that add linearly. Set $$ X=\frac{1}{a^{2}},\qquad Y=\frac{1}{b^{2}},\qquad Z=\frac{1}{c^{2}}, $$ and for a given reflection $(h,k,l)$ write $$ W=\frac{1}{d^{2}}=h^{2}X+k^{2}Y+l^{2}Z. $$ If $a,b,c$ vary independently from crystallite to crystallite (or over whatever ensemble you are sampling), then $X,Y,Z$ are independent positive random variables and each measured peak gives you the distribution of a different linear combination $W$ with known nonnegative coefficients. Your asymmetric $W$ is then entirely natural: even if $X,Y,Z$ were symmetric, the sum with unequal coefficients will generally be skewed.

There are two complementary ways to go from the observed families $\{W_{r}\}$ to information about $X,Y,Z$. The first uses moments (or, better, cumulants). Independence implies additivity of cumulants for linear combinations: for any order $m\ge1$, $$ \kappa_m(W)=h^{2m}\kappa_m(X)+k^{2m}\kappa_m(Y)+l^{2m}\kappa_m(Z). $$ Thus, if you estimate $\kappa_m(W)$ from the line shape of several distinct reflections $(h,k,l)$, you obtain (for each fixed $m$) a linear system in the three unknowns $\kappa_m(X),\kappa_m(Y),\kappa_m(Z)$. Three reflections with linearly independent coefficient triples $(h^{2m},k^{2m},l^{2m})$ are enough to solve for the $m$-th cumulants; using many reflections gives an overdetermined system you can solve in least squares to reduce noise. Doing this for $m=1,2,3,\dots$ gives you the mean, variance, skewness, etc. of $X,Y,Z$. With those cumulants in hand, you can either report them directly or fit a plausible positive family to each of $X,Y,Z$ (gamma, lognormal, inverse-Gaussian, …) by matching a few low-order cumulants. This immediately accounts for the asymmetry you see in $W$: e.g. if $a,b,c$ are approximately lognormal, then $X,Y,Z$ are also lognormal (since $\log X=-2\log a$), and a scaled sum of independent lognormals is well known to be skewed to the right.

The second route works in the Fourier (characteristic-function) domain and is essentially a deconvolution with known scalings. Let $\phi_U(t)=\mathbb{E}[e^{itU}]$ be the characteristic function. For each measured reflection $r$ with indices $(h_r,k_r,l_r)$, $$ \phi_{W_r}(t)=\phi_X(h_r^{2}t)\,\phi_Y(k_r^{2}t)\,\phi_Z(l_r^{2}t). $$ Taking logs turns products into sums: $\log\phi_{W_r}(t)=\log\phi_X(h_r^{2}t)+\log\phi_Y(k_r^{2}t)+\log\phi_Z(l_r^{2}t)$. If you estimate $\phi_{W_r}$ (or $\log\phi_{W_r}$) empirically from each line shape, these are functional equations for the three unknown functions $\log\phi_X,\log\phi_Y,\log\phi_Z$ evaluated at scaled arguments. With a suite of reflections spanning different $(h^2,k^2,l^2)$, you can fit $\log\phi_X,\log\phi_Y,\log\phi_Z$ on a grid of $t$-values by least squares with smoothness/positivity regularization; then invert the Fourier transform to obtain nonparametric estimates of the densities of $X,Y,Z$. This is a “blind” but identified deconvolution because the scalings are known and differ across reflections; the log transform makes the estimation additive and stabilizes products.

Both approaches extend naturally to error models and finite resolution. In practice a robust workflow is: estimate the first two or three cumulants $\kappa_m(W_r)$ for a moderate number of well-separated reflections; solve the linear systems for $\kappa_m(X),\kappa_m(Y),\kappa_m(Z)$ using independence; fit a positive family to each of $X,Y,Z$ (gamma and lognormal are physically reasonable for positive geometric quantities); and finally predict the full peak shapes $W_r=h_r^{2}X+k_r^{2}Y+l_r^{2}Z$ from the fitted $X,Y,Z$ as a check. If the parametric fit is inadequate, switch to the characteristic-function fitting to recover $X,Y,Z$ nonparametrically, enforcing smoothness and nonnegativity. In either case you are no longer forcing normality on $1/d^{2}$; the observed asymmetry is explained by the positivity and skewness of the underlying $1/a^{2},1/b^{2},1/c^{2}$, and the multiple reflections provide exactly the independent linear combinations needed to recover their distributions.