In this notebook, we will expose the basic theory of scattering in the distorted wave Born approximation for scalar interactions. This means we will not consider magnetic neutron scattering. A seperate lecture will be devoted to polarized neutron scattering.
After giving a general formalism for scattering, using the Green's function method, we first show how this can be applied to transmission geometry experiments. This will lead us to the Born approximation.
The major part will consider small angle scattering for multilayer samples, where a different approximation will be necessary. We will also show how the concepts of form factor and interference function will arise in this context.
Finally, we conclude with a note on how the detected intensity in a detector pixel is related to the differential cross section.
The stationary Schrödinger equation: $$\left[\frac{\hbar^2}{2m}\nabla^2 + E \right]\Psi(\mathbf r) = V(\mathbf r)\Psi(\mathbf r)$$
This can be rewritten to explicitly show the presence of the Helmholtz operator:
where $k$ is the wavenumber and $-v(\mathbf r)/4\pi$ is the scattering lenght density.
Suppose we need to solve the following inhomogeneous PDE: $$L\,y(\mathbf r) = f(\mathbf r)$$
where $L$ is a linear operator. We also assume here that appropriate boundary conditions are given.
We can first try to solve a version of this PDE without $f(\mathbf r)$ and use its solution to find $y(\mathbf r)$ in the presence of $f(\mathbf r)$.
The method of Green's functions does exactly that. The Green's function for the operator $L$ is defined as the solution to:
One can then construct a solution $y(\mathbf r)$ to the original problem: $$y(\mathbf r) = \int d^3s \, G(\mathbf r,\mathbf s)f(\mathbf s)$$
This will now satisfy the original PDE, since: $$L\,y(\mathbf r) = L\,\int d^3s\, G(\mathbf r,\mathbf s)f(\mathbf s) = \int\, d^3s\, \delta(\mathbf r -\mathbf s)f(\mathbf s) = f(\mathbf r)$$
The only thing remaining is to ensure that the solution $y(\mathbf r)$ obeys the right boundary conditions. There is a more general formulation of the Green's function method that automatically creates a solution with the correct boundary conditions. However, often this can be handled in an ad hoc manner, which is what we will also do in the following sections.
When the potential $V(\mathbf r)$ is non-zero only in a small localized region (where the sample sits), we need to find the Green's function of the Helmholtz operator.
The Green's function for the Helmholtz operator is explicitly known: $$G_0^+(\mathbf r,\mathbf s) = -\frac{e^{ik \left\| \mathbf r -\mathbf s \right\|}}{4\pi \left\|\mathbf r -\mathbf s\right\|}$$
This the advanced or causal Green's function. The retarded Green's function has a similar expression, but with a negative exponent.
We can apply here a far field approximation by putting the origin of the coordinate system in the sample (where $v(\mathbf s)\neq 0$) and assuming that we measure the neutron at a much larger distance $\mathbf r$: $$G_0^+(\mathbf r,\mathbf s) \approx -\frac{e^{ikr}\, e^{-i\mathbf{k_f} \cdot \mathbf s}}{4\pi r}$$ with $\mathbf{k_f}=k \mathbf r /r$ the outgoing wavevector.
Using the Green's function. the solution to the Schrödinger equation can be written as:
with $\phi_0(\mathbf r)$ a solution of the homogeneous equation (when $v(\mathbf r)=0$) that will determine the boundary condition satisfied by $\Psi(\mathbf r)$.
Using the far field approximation for the Green's function then produces the Lippmann-Schwinger equation:
At this point, it is interesting to note that the Green's function can be written as
with $\phi_0(\mathbf s; -\mathbf{k_f}) = e^{-i\mathbf{k_f} \cdot \mathbf s}$ a solution for the Schrödinger equation with $v(\mathbf r)=0$.
When the total wavefunction can be written as: $$\Psi(\mathbf r) = \phi_0(\mathbf r) + f(\theta,\varphi)\frac{e^{ikr}}{r}$$
with $\phi_0(\mathbf r)$ the solution without scattering, one defines the scattering amplitude as the factor $f(\theta,\varphi)$, which depends only the angular coordinates.
Clearly, our solution already has this form and the scattering amplitude is
From the expression, it is also clear that the amplitude only depends on the angular coordinates of $\mathbf r$ (through $\mathbf{k_f}$).
If we calculate the scattered particle flux (normalized to the incoming flux), we get the differential scattering cross section: $$\frac{d\sigma}{d\Omega} = \left|f(\theta,\varphi)\right|^2$$
The plane waves $\phi_0(\mathbf r) = e^{i\mathbf{k_i}\cdot\mathbf r}$, with $\|\mathbf{k_i}\|=k$, form a basis for the homogeneous solutions, as can be easily seen by plugging them into the original Schrödinger equation with $v(\mathbf r)=0$.
If we assume the scattering term is weak, compared to the incoming wave, we obtain a first order approximation of the total wavefunction:
with $\mathbf q = \mathbf{k_i} - \mathbf{k_f}$ the wavevector transfer of the scattering process.
This gives the scattering amplitude $$f(\mathbf q) = -\frac{1}{4\pi}\int d^3s\, e^{i\mathbf q \cdot \mathbf s} v(\mathbf s)$$
which is just the three dimensional Fourier transform of the scattering length density. Accordingly, the differential scattering cross section is:
In GISAS experiments, the samples consist of a multilayer structure, particles and possibly rough interfaces between the layers. Because the effect of the multilayer itself on the scattering can be quite strong (eg. total reflection), we must find a different way of calculating the cross section.
We can split the scattering potential in two parts: one that is constant in each layer and a part that contains the perturbations:
where we defined the z-axis to be perpendicular to the sample plane.
In the distorted wave Born approximation, the multilayer structure (without roughness or particles) is directly incorporated into both the Green's functions and the homogeneous solutions to the wave equation.
More details will be provided in a separate lecture on reflectivity. Here we only provide the result. The homogeneous solutions can be written as a sum of two plane waves in each layer $j$: $$\phi(\mathbf r) = \sum_j A_j^+ e^{i\mathbf{k_j^+} \cdot \mathbf r} + A_j^- e^{i\mathbf{k_j^-} \cdot \mathbf r}$$ where the amplitudes $A_j^\pm$ and wavevectors $\mathbf{k_j^\pm}$ are to be determined by solving the boundary conditions. In the expression above, we assume that these amplitudes are only defined inside their corresponding layer (equivalently, one can set $A_j^\pm = 0$ outside layer $j$).
To denote the dependence of these solutions to the boundary condition of the incoming plane wave, we put the incoming wavevector $\mathbf{k}$ explicitly:
In the region where the beam enters (i.e. above the sample surface), it obeys the following boundary condition:
where $R$ is now the familiar reflectivity amplitude.
For reasons that have to do with the self-adjointness of the Hamiltonian on the Hilbert space of wavefunctions, the far field approximation of the Green's function for the multilayer can also be written as $$G^+(\mathbf r,\mathbf s) \approx -\frac{e^{ikr}}{4\pi r}\phi(\mathbf s; -\mathbf{k_f})$$
with $\phi(\mathbf s; -\mathbf{k_f})$ a solution for the Schrödinger equation with $\delta v(\mathbf r)=0$ and incoming wavevector $-\mathbf{k_f}$.
This leads to a similar expression of the scattering amplitude as in the case of the Born approximation:
This is the scattering amplitude in the distorted wave Born approximation. Compared to the Born approximation, we now need to consider four different terms for each layer, each consisting of a Fourier transform with a certain combination of wavevectors.
Because of translation invariance in the in-plane directions, all wavevectors (in each layer $j$) $\mathbf{k_j^\pm}$ of a homogeneous solution $\phi(\mathbf r;\mathbf k)$ will have the same projections in the sample plane: $$\mathbf{k_{jx}^\pm} = \mathbf{k_x}$$ $$\mathbf{k_{jy}^\pm} = \mathbf{k_y}$$ where $\mathbf{k}$ is again the wavevector of the incoming plane wave.
For a single particle, inside a layer $j$, the scattering amplitude becomes:
where the $\mathbf{q_m}$ are the four different wavevectors transfers between the two incoming and two outgoing plane waves in layer $j$ and the $C_m$ are products of their complex amplitudes.
For a particle with a constant scattering length density, the scattering amplitude is a weighted sum of Fourier transforms of the shape function, which is defined as being zero outside and one inside the shape. We will denote this as a DWBA form factor:
with $V_S$ the particle's volume and $S(\mathbf 0)$ denoting the particle at a default position (zero displacment).
Because all the wavevector tranfers have the same in plane components, a displacement along the in plane directions $\mathbf{R_\parallel}$ can be written as
For $N$ particles with form factors $\mathcal F_i$ and sitting at positions $\mathbf R_i$, the scattering amplitude is then:
where we dropped a few indices and dependences for ease of notation.
The differential scattering cross section is: $$\frac{d\sigma}{d\Omega}(\mathbf{k_i}, \mathbf{k_f}) = \sum_{i,j} \mathcal F_i \mathcal F_j^* e^{i\mathbf q\cdot (\mathbf{R_i}-\mathbf{R_j})}$$
If we take the expectation value of this cross section over the types and positions of the particles, we have:
where $P_{i,j}(\alpha, \beta, \mathbf{R_1}, \mathbf{R_2})$ is the probability density function for the particles $i$ and $j$ to be of types $\alpha$ and $\beta$ and to be at positions $\mathbf{R_1}$ and $\mathbf{R_2}$.
When the positions of the particles are decoupled with their types, the probability density factorizes:
The expected cross section then becomes:
where $\langle . \rangle_\alpha$ denotes an expectation value over the type $\alpha$.
This can be rewritten as: $$\frac{1}{N}\left\langle\frac{d\sigma}{d\Omega}(\mathbf{k_i}, \mathbf{k_f})\right\rangle = \left\langle \left| \mathcal F_\alpha \right|^2 \right\rangle_\alpha - \left|\left\langle \mathcal F_\alpha \right\rangle_\alpha \right|^2 + \left|\left\langle \mathcal F_\alpha \right\rangle_\alpha \right|^2 \left( 1 + \frac{1}{N}\sum_{i\neq j} \int dP_{i,j}(\mathbf{R_1}, \mathbf{R_2})\, e^{i\mathbf q\cdot (\mathbf{R_1}-\mathbf{R_2})}\right)$$
The factor in the last term is called the interference function $S(\mathbf q)$ and encodes the information about the relative positions of the particles:
Because the interference function only contains the in plane components of the wavevector transfer, it depends only on the difference $\mathbf q = \mathbf{k_i} - \mathbf{k_f}$.
The previous expression of the differential cross section is called the decoupling approximation since it assumes a decoupling between the types and positions of the particles.
We can also assume that there is only coupling between particles of the same type:
In this case, the differential cross section becomes: $$\frac{1}{N}\left\langle\frac{d\sigma}{d\Omega}(\mathbf{k_i}, \mathbf{k_f})\right\rangle = \left\langle \left|\mathcal F_\alpha \right|^2S_\alpha(\mathbf q)\right\rangle_\alpha$$
This is called the local monodisperse approximation. It can also be interpreted as consisting of different domains, each with their own particle type and interference function, that have no coherent scattering between different domains.
For rough interfaces, we can again use the expression we found for the scattering amplitude:
However, now the scattering perturbation $\delta v(\mathbf s)$ will have parts above and below the interface. This forces us to consider not just the usual four terms in the distorted wave Born approximation, but eight terms (four for each side of the interface).
We will not go into the details here, but to first order in the lateral correlation coefficient, the expectation value of the differential cross section can be found exactly for a bivariate normal distribution of the roughness height.
In an experiment, each detector pixel records the number of scattered particles, which is the integral over the pixel of the differential cross section times a normalization factor: $$I_p = I_0 \int_{\Delta\Omega} \frac{d\sigma}{d\Omega}$$ where $\Delta\Omega$ is the total solid angle of the pixel.
If the cross section does not vary over the short range of solid angle inside a single pixel, we can assume it to be constant:
where we calculate the cross section at some preferential point inside the pixel, eg. the midpoint. In BornAgain, this is referred to as the analytical computation method, since it doesn't require a non-deterministic calculation.
If the cross section fluctuates considerably inside a pixel, we can not use this approach. In BornAgain, we implemented a Monte Carlo integration over the pixel's solid angle to account for this effect.