%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%   BornAgain Developers Reference
%%
%%   homepage:   http://www.bornagainproject.org
%%
%%   copyright:  Forschungszentrum Jülich GmbH 2015-
%%
%%   license:    Creative Commons CC-BY-SA
%%
%%   authors:    Scientific Computing Group at MLZ Garching
%%   editor:     Joachim Wuttke
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{Flat multilayer systems}\label{sec:Multilayers}%
\index{Multilayer|(}%
\index{Layered structure|see{Multilayer}}

This chapter specializes the DWBA for a multilayer system with
$\overline{v}(\r)=\overline{v}(z)$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Wave propagation and scattering in layered samples}\label{Swave21}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%===============================================================================
\subsection{Wave propagation in 2+1 dimensions}\label{Sgrazingwave}
%===============================================================================

We now specialize the results from~\cref{SSca} to wave propagation
in a sample that is, on average, translationally invariant in 2~dimensions.
Following standard convention,
we choose the surface of the sample in the $xy$~plane,
\index{Sample plane}%
and its normal along~$z$.
\index{Sample normal}%
\nomenclature[2x020]{$x$}{Horizontal coordinate, in the sample plane}%
\nomenclature[2y020]{$y$}{Horizontal coordinate, in the sample plane}%
\nomenclature[2z020]{$z$}{Vertical coordinate, along the sample normal}%
\index{Horizontal plane}%
\index{Vertical direction}%
In visualizations, we will always represent the $xy$~plane as \E{horizontal},
and the $z$~axis as upward \E{vertical},
altough there are ``horizontal'' reflectometers
where the sample is upright to allow for a horizontal scattering plane.
\index{Reflectometer!vertical vs horizontal}%

Scattering from such systems will be studied in distorted-wave Born approximation.
To determine the neutron scattering cross section~\cref{Exsection},
we need to determine the incident and final wavefunctions
$\psi_\si$ and~$\psi_\sf$.
Vertical variations of the refractive index $n(z)$
\index{Refractive index!vertical variation}%
cause refraction and reflection.
\index{Glancing angle}%
\index{Refraction}%
\index{Reflection}%
For waves propagating at small glancing angles,
the reflectance can take any value between 0 and~1,
even though $1-n$ is only of the order $10^{-5}$ or smaller.
Such zeroth-order effects cannot be accounted for
by perturbative scattering theory.
Instead, we need to deal with refraction and reflection
from the onset, in the wave propagation equation.
Accordingly, the SLD decomposition~\cref{Edecompose} takes the form
\begin{equation}\label{Edecompose_z}
  v(\r) = \mv(z) + \delta v(\r),
\end{equation}
\index{Wave propagation!in multilayer|(}
and the unperturbed distorted wave equation~\cref{EDPsi0} becomes
\begin{equation}\label{EWaveZ}
  \left\{\Nabla^2+k(z)^2\right\}\psi(\r) = 0.
\end{equation}
Below and above the sample,
$k(z)=\text{const}$:
in these regions, $\psi(\r)$~is a superposition of plane waves.
The exciting wavefunction is
\begin{equation}\label{Epsiminus}
  \psi_\se(\r) = \e^{i\k_\plll\r_\plll+ik_{\perp\se}z},
\end{equation}
\nomenclature[0$\plll$]{$\plll$}{Parallel to the $xy$ sample plane}%
\nomenclature[0$\perp$]{$\perp$}{Normal to the $xy$ sample plane}%
\nomenclature[2k021\perp]{$k_\perp$}{Component of $\k$ along the sample normal}%
\nomenclature[2k041\plll]{$\k_\plll$}{Projection of $\k$ onto the sample plane}%
The subscripts $\plll$ and~$\perp$ refer to the sample $xy$ plane.
The wavevector components $\k_\plll$ and $k_{\perp}$ must fulfill
\begin{equation}
  k(z)^2=\k_\plll^2+k_{\perp}^2.
\end{equation}
\index{Wavenumber!vertical}%
\index{Vertical wavenumber}%
Continuity across the sample implies
\begin{equation}\label{Ekconst}
  \k_\plll = \text{const}.
\end{equation}
\index{Wavevector!horizontal}%
\index{Horizontal wavevector}%
From here on, we abbreviate
\begin{equation}\label{Dkappa}
  \kappa \coloneqq k_\perp.
\end{equation}
When the incident wave hits the sample,
it is wholly or partly reflected.
Therefore, the full the solution of~\cref{EWaveZ} in the half space
of the radiation source is
\begin{equation}\label{Eref1}
  \psi(\r) = \e^{i\k_\plll\r_\plll+i\kappa_\se z} +
      R\, \e^{i\k_\plll\r_\plll-i\kappa_\se z}
\end{equation}
with a complex reflection coefficient~$R$.
\index{Reflection!coefficient}%
The reflected flux is given by the re\-flect\-an\-ce $|R|^2$.
\index{Reflectance}%
\index{Flux!reflected}%
In the opposite halfspace, the solution of~\cref{EWaveZ} is simply
\begin{equation}\label{Etra1}
  \psi(\r) = T\, \e^{i\k_\plll\r_\plll+i\kappa_\se z}
\end{equation}
with a complex transmission coefficient~$T$.
The transmitted flux is given by the transmittance $|T|^2$.
\index{Transmittance}%
\index{Flux!transmitted}%
As before, subscript $\se$ stands for the exciting wave in vacuum outside the sample.

Within the sample, the wave equation~\cref{EWaveZ}
is solved by the factorization ansatz
\begin{equation}\label{Ekpar}
\psi(\r) = \e^{i \k_\plll\r_\plll} \phi(z).
\end{equation}
\nomenclature[1φ020 0z020]{$\phi(z)$}{$z$-dependent factor of $\psi(\r)$}%
The vertical wavefunction~$\phi(z)$
is governed by the one-dimensional wave equation
\begin{equation}\label{Ewavez}
\left\{\partial_z^2 + k(z)^2 - k_\plll^2 \right\} \phi(z) = 0.
\end{equation}
As solution of a differential equation of second degree,
$\phi(z)$~can be written as superposition
of a downward travelling wave $\phi^-(z)$
and an upward travelling wave $\phi^+(z)$.
Accordingly, the three-dimensional wavefunction can be written as
\begin{equation}\label{Epsisumpm}
  \psi(\r) = \psi^-(\r)+\psi^+(\r).
\end{equation}
\nomenclature[1ψ041 0\pm 2r040]{$\psi^\pm(\r)$}{Upward ($+$) or downward ($-$) propagating component of $\psi(\r)$}%
\nomenclature[0\pm]{$\pm$}{Upward ($+$) or downward ($-$) propagating}%

%===============================================================================
\subsection{The four DWBA terms}\label{Sdwba4terms}
%===============================================================================

All the above holds not only for the incident wavefunction~$\psi_\si$,
but also for the wavefunction~$\psi_\sf$
that is tracked back from a detector pixel towards the sample.
Therefore the scattering matrix element
involves two incident and two final partial wavefunctions.
The resulting sum
\index{Wave propagation!in multilayer|)}
\begin{equation}\label{Edwba4}
  \braket{\psi_\si|\delta v|\psi_\sf}
  = \braket{\psi^-_\si|\delta v|\psi^-_\sf}
  + \braket{\psi^-_\si|\delta v|\psi^+_\sf}
  + \braket{\psi^+_\si|\delta v|\psi^-_\sf}
  + \braket{\psi^+_\si|\delta v|\psi^+_\sf}
\end{equation}
is depicted in \Cref{Fdwba4terms}.
It can be written in an obvious shorthand notation
\begin{equation}\label{Edwba}
  \braket{\psi_\si|\delta v|\psi_\sf}
  = \sum_{\pm_\si} \sum_{\pm_\sf}
    \braket{\psi^\pm_\si|\delta v|\psi^\pm_\sf}.
\end{equation}
This equation contains the essence of
the DWBA for GISAS,
and is the base for all scattering models implemented in \BornAgain.
Since $\braket{\psi_\si|\delta v|\psi_\sf}$
appears as a squared modulus
in the differential cross section~\cref{Exsection},
the four terms of \cref{Edwba} can interfere with each other,
which adds to the complexity of GISAS patterns.

%--------------------------------------------------------------------------------
\begin{figure}[tb]
\begin{center}
\includegraphics[width=1\textwidth]{fig/drawing/dwba_4terms.ps}
\end{center}
\caption{The four terms in the DWBA scattering matrix element~\cref{Edwba}.
Note that this is a highly simplified visualization.
In particular, it does not show multiple reflections
of incoming or scattered radiation,
though they are properly accounted for by DWBA theory and by all simulation software.}
\label{Fdwba4terms}
\end{figure}
%--------------------------------------------------------------------------------

BornAgain supports multilayer samples
with refractive index discontinuities at layer interfaces.
Conventions for layer numbers and interface coordinates are introduced in~\Cref{Fdefz}.
\index{Coordinates!sample}%
\index{Interface!coordinate}%
\index{Multilayer!numbering}%
\index{Multilayer!coordinates}%
\index{Layer!index}%
\index{Numbering!layers}%
A sample has $N$ layers,
including the semi-infinite bottom and top layers.
Numbering is from top to bottom,
and from 0 to $N-1$ as imposed by the programming languages C$++$ and Python.
Each layer~$l$
\nomenclature[2l010]{$l$}{Layer index}%
has a constant refractive index $n_l$
\nomenclature[2k022 2l010]{$k_l$}{Wavenumber in layer~$l$}%
\nomenclature[2n020 2l010]{$n_l$}{Refractive index of layer~$l$}%
and a constant wavenumber $k_l\coloneqq K_\text{vac} n_l$.
Any up- or downward travelling solution of the wave equation shall be written
as a sum over partial wavefunctions,
\begin{equation}\label{Epsipmsuml}
  \psi^\pm(\r) = \sum_l \psi_l^\pm(\r),
\end{equation}
with the requirement
\begin{equation}\label{Epsipmloutside}
   \psi_l^\pm(\r) = 0 \text{~for $\r$ outside layer~$l$.}
\end{equation}
The DWBA matrix element~\cref{Edwba} then takes the form
\begin{equation}\label{Edwbal}
  \braket{\psi_\si|\delta v|\psi_\sf}
  = \sum_l \sum_{\pm_\si} \sum_{\pm_\sf}
    \braket{\psi^\pm_{\si l}|\delta v|\psi^\pm_{\sf l}}.
\end{equation}

%--------------------------------------------------------------------------------
\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.4\textwidth]{fig/drawing/multilayer_z_conventions.ps}
\end{center}
\caption{Conventions for layer numbers and interface coordinates.
\index{Coordinates!sample}%
\index{Interface!coordinate}%
\index{Multilayer!numbering}%
\index{Multilayer!coordinates}%
\index{Layer!index}%
\index{Numbering!layers}%
A sample has $N$ layers,
including the semi-infinite bottom and top layers.
\nomenclature[2n120]{$N$}{A multilayer sample has $N$ layers, including the
  semi-infinite bottom and top layers}
Layers are numbered from top to bottom.
The top vacuum (or air) layer (which extends to $z\to+\infty$) has number~0,
the substrate (extending to $z\to-\infty$) is layer~$N-1$.
The parameter $z_l$
\nomenclature[2z020 2l010]{$z_l$}{Vertical
  coordinate at the top of layer~$l$ (at the bottom for $l=0$)}%
is the $z$ coordinate of the \E{top} interface of layer~$l$,
except for $z_0$ which is the coordinate of the \E{bottom} interface
of the air or vacuum layer~0.}
\label{Fdefz}
\end{figure}
%--------------------------------------------------------------------------------

%===============================================================================
\subsection{DWBA for layers with constant mean SLD}\label{SStep}
%===============================================================================

We now specialize to the case that $\mv(z)$ is a step function:
within each layer, $\mv(z)\eqqcolon v_l$ is constant.
Accordingly, within the layer, the directional neutron wavefunction~$\psi^\pm_l$
is a plane wave and factorizes as in~\cref{Ekpar}.
Its amplitude~$A_l^\pm$ is determined recursively
by Fresnel's transmission and reflection coefficients
\index{Fresnel coefficients}%
that are based on continuity conditions at the layer interfaces.
This will be elaborated in \Cref{Sacrolay}.
\index{Multilayer!refractive index profiles}%
\index{Layer!refractive index profiles}%
The vertical wavenumber is determined by \cref{Epsiminus} and~\cref{Ekconst},
\begin{equation}\label{Ekperpl}
  \kappa_l^\pm = \pm\sqrt{k_l^2 - k_\plll^2}.
\end{equation}
In the absence of absorption and above the critical angle,
wavevectors are real
so that we can describe the beam in terms of a glancing angle
\begin{equation}\label{Edef_alpha}
  \alpha_l\coloneqq \arctan(\kappa_l/k_{\plll}).
\end{equation}
Equivalently,
\begin{equation}\label{Ekplllncos}
  k_{\plll}=K n_l \cos\alpha_l.
\end{equation}
Since $k_{\plll}$ is constant across layers,
we have
\begin{equation}\label{ESnell}
  n_l \cos\alpha_l = \text{the same for all }l,
\end{equation}
which is Snell's refraction law.
\index{Refraction!Snell's law}
\index{Snell's law}
In general, however, the vertical wavenumber $\kappa_l$,
determined by $k_l$ and $k_\plll$ as per~\cref{Epsiminus},
can become imaginary (total reflection conditions) or complex (absorbing layer).
\index{Wavevector!complex}%
In these cases, glancing angles are no longer well defined,
and the geometric interpretation of~$\psi_l(\r)$ less obvious.
so that one has to fully rely on the algebraic formalism.

With the indicator function
\nomenclature[1χ032 2l010]{$\chi_l(z)$}{Indicates whether $z$ is in layer~$l$}%
\begin{equation}\label{Echildef}
  \chi_l(\r)\coloneqq\left\{\begin{array}{ll}
  1&\text{~if $z_l\le z \le z_{l+1}$,}\\[.2ex]
  0&\text{~otherwise,} \end{array}\right.
\end{equation}
the vertical wavefunction can be written
\begin{equation}\label{Ephizwj}
  \phi^\pm_l(z)=A^\pm_l\e^{\pm i\kappa_l(z-z_l)}\chi_l(z).
\end{equation}
\nomenclature[2a123 2w010 2l010 \pm]{$A^\pm_{wl}$}{Amplitude
  of the plane wave $\phi^\pm_{wl}(\r)$}%
The offset~$z_l$ has been included in the phase factor for later convenience.
\iffalse See \cref{Snokz} for the case of vanishing~$\kappa$.\fi

The DWBA transition matrix element~\cref{Edwba} is
\index{DWBA!multilayer}%
\begin{equation}\label{Edwba_ml0}
  \braket{\psi_\si|\delta v|\psi_\sf}
  = \sum_l \sum_{\pm_\si} \sum_{\pm_\sf}
    A^{\pm *}_{\si l} A^\pm_{\sf l}
     \delta v_l(\k^\pm_{\sf l}-\k^\pm_{\si l})
\end{equation}
with the Fourier transform of the SLD
restricted to layer~$l$
\begin{equation}\label{Echij}
  \delta v_l(\v{q})
  \coloneqq  \int_{z_l}^{z_{l-1}}\!\d z \int\!\d^2r_\plll\, \e^{i\v{q}\,\r}\delta v(\r)
  = \int\!\d^3r\, \e^{i\v{q}\,\r}\delta v(\r) \chi_l(z).
\end{equation}
\nomenclature[1δ00 2v030 2l010 2q040]{$\delta v_l(\v{q})$}{Fourier transform
of the SLD $\delta v(\r)$, evaluated in one sample layer}%
To alleviate later calculations,
we number the four DWBA terms from 1 to~4 as shown in \cref{Fdwba4terms},
and define the corresponding wavenumbers and amplitude factors and as
\begin{equation}\label{Eudef}
  \begin{array}{l@{\hspace{2em}}l}
    \q^1 \coloneqq  \k^+_\sf - \k^-_\si,& C^1 \coloneqq  A^{-*}_\si A^+_\sf, \\[.6ex]
    \q^2 \coloneqq  \k^-_\sf - \k^-_\si,& C^2 \coloneqq  A^{-*}_\si A^-_\sf, \\[.6ex]
    \q^3 \coloneqq  \k^+_\sf - \k^+_\si,& C^3 \coloneqq  A^{+*}_\si A^+_\sf, \\[.6ex]
    \q^4 \coloneqq  \k^-_\sf - \k^+_\si,& C^4 \coloneqq  A^{+*}_\si A^-_\sf.
  \end{array}
\end{equation}
Accordingly, we can write \cref{Edwba_ml0} as
\begin{equation}\label{Edwba_ml}
  \braket{\psi_\si|\delta v|\psi_\sf}
  = \sum_l \sum_{u} C^u_l \delta v_l(\q_l^u).
\end{equation}
Since $\k_\plll=\text{const}$,
 all wavevectors $\q^u_l$ have the same horizontal component~$\q_\plll$;
they differ only in their vertical component~$q^u_{l\perp}$.

%===============================================================================
\subsection{Wave amplitudes}\label{Sacrolay}
%===============================================================================

\index{Fresnel coefficients}%
\index{Transmission|see{Fresnel coefficients}}%
\index{Reflection|seealso {Fresnel coefficients}}%

The plane-wave amplitudes $A^\pm_{wl}$ need to be computed recursively
from layer to layer.
Since these computations are identical for incident and final waves,
we omit the subscript~$w$ in the remainder of this section.
At layer interfaces, the optical potential changes discontinuously.
From elementary quantum mechanics we know that
piecewise solutions of the Schrödinger equations must be connected
such that the wavefunction $\phi(\r)$ and its first derivative
$\Nabla\phi(\r)$ evolve continuously.

%--------------------------------------------------------------------------------
\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.46\textwidth]{fig/drawing/multilayer_boundary.ps}
\end{center}
\caption{The transfer matrix $M_l$ connects the wavefunctions
\index{Transfer matrix}%
$\Phi_l$, $\Phi_{l-1}$ in adjacent layers.}
\label{Fboundary}
\end{figure}
%--------------------------------------------------------------------------------

To deal with the coordinate offsets introduced in \cref{Ephizwj},
we introduce the function%
\begin{equation}\label{Edldef}
  d_l\coloneqq z_l-z_{l+1},
\end{equation}
which is the thickness of layer~$l$,
except for $l=0$,
where the special definition of $z_0$ (\cref{Fdefz}) implies $d_0=0$.
We consider the interface between layers $l$ and $l-1$,
with~$l=1,\ldots,N-1$, as shown in \cref{Fboundary}.
This interface has the vertical coordinate $z_l=z_{l-1}-d_{l-1}$.
Accordingly, the continuity conditions at the interface are
\begin{equation}\label{Econtcond}
  \begin{array}{lcl}
 \hphantom{\partial_z}\phi_l(z_l) &=& \hphantom{\partial_z}\phi_{l-1}(z_{l-1}-d_{l-1}),\\
           \partial_z \phi_l(z_l) &=&           \partial_z \phi_{l-1}(z_{l-1}-d_{l-1}).
  \end{array}
\end{equation}
We define the phase factor
\begin{equation}\label{Ddell}
   \delta_l \coloneqq  \e^{i\kappa_l d_l}.
\end{equation}
Here and in the following, we will write the downward travelling transmitted
and of the upward travelling reflected amplitude as
\begin{equation}
  t_l \coloneqq A^-_l \quad\text{and}\quad r_l \coloneqq A^+_l.
\end{equation}
For the plane waves \cref{Ephizwj},
the continuity conditions~\cref{Econtcond} take the form
\begin{equation}\label{Econt2}
  \begin{array}{@{}l@{}lcl@{}l}
  \hphantom{+}t_l &\;+\;r_l
  &=&
  \hphantom{+}\delta_{l-1} t_{l-1} &\;+\; \delta_{l-1}^{-1} r_{l-1},
  \\[1.1ex]
  -\kappa_l t_l   &\;+\;  \kappa_l r_l
  &=&
  -\kappa_{l-1} \delta_{l-1} t_{l-1} &\;+\; \kappa_{l-1}\delta_{l-1}^{-1} r_{l-1} .
  \end{array}
\end{equation}
After some lines of linear algebra,
we can rewrite this equation system as
\begin{equation}\label{EcMc}
  \left( \begin{array}{c}t_{l-1}\\ r_{l-1}\end{array} \right)
  = M_l \left( \begin{array}{c}t_l\\r_l\end{array} \right)
\end{equation}
with the transfer matrix\footnote
{This approach is generally attributed to Abelès,
\index{Abelès matrix}%
who elaborated it in his thesis from 1949, published 1950.
The usually cited paper \cite{Abe50a} is no more than a short advertisement.}
\index{Transfer matrix}%
\begin{equation}\label{EMil}
  M_l \coloneqq \Delta_{l-1} S_l,
\end{equation}
which we write using the phase rotation matrix
\begin{equation}\label{DmatD}
  \Delta_l
   \coloneqq
   \left(\begin{array}{cc}
     \delta_{l}^{-1}&0\\
       0 & \delta_{l}
   \end{array}\right)
\end{equation}
and the refraction matrix
\begin{equation}\label{DmatS}
  S_l
   \coloneqq
   \left(\begin{array}{cc}
       s^+_l&s^-_l\\
       s^-_l&s^+_l
   \end{array}\right)
\end{equation}
with coefficients
\begin{equation}\label{Dslpm}
  s^\pm_l \coloneqq \frac{1 \pm \kappa_l/\kappa_{l-1}}{2}.
\end{equation}
Energy conservation can be easily verified for real-valued wave numbers.
The vertical flux is $J=|\Phi|^2\kappa$.
Under the action of either $\Delta$ or S,
\begin{equation}\label{EConservation}
  \kappa_l (|t_l|^2 - |r_l|^2) = \text{const for all~$l$}.
\end{equation}

%===============================================================================
\subsection{Wave amplitudes for X-rays}\label{SmulayX}
%===============================================================================

\def\Ep{\v{\Phi}}
\def\hn{\v{n}}

We shall now translate the above results from unpolarized neutrons to X-rays.
The vectorial amplitude of the electromagnetic field will require
nontrivial modifications.
In place of the factorization~\cref{Ekpar}, we write
\begin{equation}
  \v{E}(\r)=\e^{i\k_\plll\r}\Ep(z).
\end{equation}
In place of~\cref{Ephizwj},
the vertical wavefunction is
\begin{equation}
  \Ep^\pm_l(z) = \v{A}^\pm_l \e^{\pm i\kappa(z-z_l)}\chi_l(z).
\end{equation}

%--------------------------------------------------------------------------------
\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.8\textwidth]{fig/drawing/s-vs-p-polarization.ps}
\end{center}
\caption{Conventions for polarization directions relative to a refracting interface:
\index{Polarization!X-ray@X-ray ($s$ and $p$)}%
\index{X-ray!polarization}%
For $s$ polarization, the electric field vector~$\v{E}$ is perpendicular (\E{senkrecht} in German)
 to the plane spanned by the interface normal~$\hn$ and the incoming~wavevector~$\k$;
for $p$ polarization, it is parallel.
In either case, $\v{E}$ is perpendicular to~$\k$.}
\label{Fsppol}
\end{figure}
%--------------------------------------------------------------------------------

The vectorial character of $\v{A}^\pm_{wl}$ will require changes
with respect to~\cref{Sacrolay}.
For electromagnetic radiation in nonmagnetic media,
the boundary conditions at an interface with normal $\hn$ are \cite[eq. 7.37]{Jac75}
% , Born \& Wolf \cite[ch.~1.1.3]{BoWo99}, or Hecht \cite[ch.~4.6.1]{Hec02}.}
\nomenclature[2n04]{$\hn$}{Normal vector of an interface}
\begin{align}
  &\sum_\pm\,\overline{\epsilon}\,\v{E}^\pm\,\hn = \text{const}, \label{EbcE1}\\[1.4ex]
  &\sum_\pm\,\v{E}^\pm\times\hn = \text{const}, \label{EbcE2}\\[1.4ex]
  &\sum_\pm\,\k^\pm_l\times\v{E}^\pm = \text{const}. \label{EbcE3}
\end{align}
We will only consider the two polarization directions,
\index{Polarization!X-ray@X-ray ($s$ and $p$)|(}%
\index{X-ray!polarization|(}%
conventionally designated as $p$ and~$s$, defined in \Cref{Fsppol}.
As some algebra on \cref{EbcE1,EbcE2,EbcE3} would show,
these are \E{principal axes},
meaning that if both incoming fields $\v{E}^-_{l-1}$ and~$\v{E}^+_l$ are strictly
polarized in either $s$ or $p$ direction,
then the outgoing fields $\v{E}^+_{l-1}$ and~$\v{E}^-_l$
are polarized in the same direction.
Conversely, if the incoming fields are mixtures of $s$ and $p$ polarization,
then the outgoing fields will be, in general, mixed differently.
Therefore if polarization factors are quantitatively important in an experiment,
one should strive to accurately polarize the incident beam in $s$ or $p$ direction
in order to avoid the extra complication of variably mixed polarizations.

Further algebra on \cref{EbcE1,EbcE2,EbcE3} replicates the
reflection law that relates $\k^-$ and $k^+$
and Snell's law~\cref{ESnell}.
Taking these for granted,
we only retain equations that are needed to determine the field amplitudes~$E^\pm$.
For $s$~polarization they yield
\begin{equation}
  \left(\begin{array}{cc}1&1\\
       -\kappa&\kappa\end{array}\right)
  \left(\begin{array}{c}E^-\\
       E^+\end{array}\right) = \text{const}.
\end{equation}
and for $p$~polarization
\begin{equation}
  \left(\begin{array}{cc}n&n\\
       -\kappa/n&\kappa/n\end{array}\right)
  \left(\begin{array}{c}E^-\\
       E^+\end{array}\right) = \text{const},
\end{equation}
The former equation can be brought into the form~\cref{Econt2}.
In consequence,
$s$-polarized X-rays are refracted and reflected in
exactly the same ways as unpolarized neutrons.

For $p$ polarization,
the refraction matrix coefficients become
\begin{equation}
  s_l^\pm = \frac{1}{2}
       \left(\frac{n_l}{n_{l-1}} \pm \frac{\kappa_l}{\kappa_{l-1}}\frac{n_{l-1}}{n_l}\right)
\end{equation}
instead of \cref{Dslpm}.\footnote
{Support for $p$ polarization is not implemented in BornAgain.
It can be added easily if there is need.}
\index{Polarization!X-ray@X-ray ($s$ and $p$)|)}%
\index{X-ray!polarization|)}%

%===============================================================================
\subsection{Scattering of X-rays}\label{SscatterX}
%===============================================================================

The DWBA matrix element is
\begin{equation}\label{Edwba_mlE}
  \braket{\v{E}_\si|\delta v|\v{E}_\sf}
  = \sum_l \sum_{u} C^u_l \delta v_l(\q_l^u).
\end{equation}
in full analogy with~\cref{Edwba_ml},
but the coefficients
$C^1=\v{A}^{-*}_\si\v{A}^{+}_\sf$ etc are now scalar products
of vectorial amplitudes.
For $s$~polarization all amplitudes point in the same direction,
so that we are back to the products of scalar factors of~\cref{Eudef}.
For $p$~polarization,
incident and scattered field amplitudes point in slightly different directions,
which results in correction factors\footnote
{Also currently not implemented in BornAgain.}
\begin{equation}\label{Edwbap}
  \begin{array}{@{}lcl}
    C^1 &=& A^{-*}_\si A^+_\sf\cos(\alpha^{-}_\si+\alpha^+_\sf),\\
    C^2 &=& A^{-*}_\si A^-_\sf\cos(\alpha^{-}_\si+\alpha^-_\sf),\\
    C^3 &=& A^{+*}_\si A^+_\sf\cos(\alpha^{+}_\si+\alpha^+_\sf),\\
    C^4 &=& A^{+*}_\si A^-_\sf\cos(\alpha^{+}_\si+\alpha^-_\sf).
  \end{array}
\end{equation}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Solution of the split boundary problem}\label{Ssolvsplit}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%==================================================================================================%
\subsection{The split boundary problem}\label{Ssplibou}
%==================================================================================================%

We now consider beam propagation through the entire multilayer sample,
from the semiinfinite top layer at $l=0$ to the semiinfinite substrate at $l=N-1$,
which for brevity shall be denoted by $\nu\coloneqq N-1$.

Let us assume that the radiation source or sink is located at $z>0$.
Then in the top layer, $t_0 =1$ is given by the
incident or back-traced final plane wave.
In the substrate, $t_{\nu} =0$ because there is no radiation
coming from $z\to-\infty$.
This leaves us with two unkown amplitudes,
the overall coefficients of transmission~$t_{\nu}$ and reflection~$r_0$.
These two unknowns are connected by a system of two linear equations,
\begin{equation}\label{E1Ap}
  \left( \begin{array}{c}1\\ r_0\end{array} \right)
  = M \left( \begin{array}{c}t_{\nu}\\0\end{array} \right)
\end{equation}
with the matrix product
\begin{equation}\label{DM22}
  M \coloneqq M_1 \cdots M_{\nu}
   \eqqcolon
   \begin{pmatrix} M_{tt} &M_{tr}\\ M_{rt} &M_{rr}\end{pmatrix}.
\end{equation}
To apply this and all the following to the scattered beam in transmission GISAS
\index{Transmission geometry}%
\index{Detector!transmission geometry}%
(sink location $z<0$),
we just reverse the order of layers:
$(0,\ldots, \nu) \mapsto (\nu,\ldots,0)$.

Equation~\cref{E1Ap} is a \emph{split boundary problem} because
the given amplitudes $t_0=1$, $r_\nu=0$ appear on different sides of the equation.
It can be reorganized as
\begin{equation}\label{EWif}
  \begin{pmatrix}{t_\nu}\\ {r_0}\end{pmatrix}
  = W \begin{pmatrix}{1}\\ {0}\end{pmatrix}
\end{equation}
with
\begin{equation}\label{EM2W}
  W = \mathcal{W}(M)
  \coloneqq \begin{pmatrix} {M_{tt}}^{-1} &{M_{tt}}^{-1}M_{tr}\\[.2em]
                     M_{rt}{M_{tt}}^{-1} &
                        (M_{rr}-M_{rt}{M_{tt}}^{-1}M_{tr}) \end{pmatrix}.
\end{equation}
% For later use, we note the inverse function
% \begin{equation}\label{EW2M}
%   M = \mathcal{M}(W)
%   = \begin{pmatrix}  (W_{tt}-W_{tr}{W_{rr}}^{-1}W_{rt}) & W_{tr}{W_{rr}}^{-1}\\[.2em]
%                      {W_{rr}}^{-1}W_{rt} &{W_{rr}}^{-1} \end{pmatrix}.
% \end{equation}
% This formalism, originally developed for dynamic X-ray diffraction \cite{Koh91,StKK98},
% holds also if the matrix components are not commutative under multiplication.
% This will allow us later (for polarized neutrons, \cref{SPol})
% to replace the scalar matrix components by $2\times2$ matrices.
%
From \cref{EWif} and \cref{EM2W}, we can read off
\begin{equation}\label{Etfri0J}
      t_\nu = M_{tt}^{-1} \quad\text{and}\quad r_0 = M_{rt} M_{tt}^{-1}.
\end{equation}
With this, the split boundary problem is formally solved.
However, the matrix product $M$ \cref{DM22} is numerically unstable
\cite[Sects.~III, IV]{StKK98}.
Therefore, the actual computation of $r_0$ and $t_\nu$ is done through
a recursion (\cref{Srt1,SParrattPol}).

If there is one single interface ($\nu=1$),
then $M=S_1$ yields the standard Fresnel results, namely the transmitted amplitude
\begin{equation}\label{EtFresnel}
  t_1 =\frac{2\kappa_0}{\kappa_0+\kappa_1}
\end{equation}
and the reflected amplitude
\begin{equation}\label{ErFresnel}
  r_0 =\frac{\kappa_0-\kappa_1}{\kappa_0+\kappa_1}.
\end{equation}
In connection with roughness models,
we will need to express the coefficients of the refraction matrix~\cref{DmatS}
through $t$ and~$r$,
\begin{equation}
  s_1^+ = \frac{1}{t_1} \quad\text{ and }\quad s_1^- = \frac{r_0}{t_1}.
\end{equation}

%==================================================================================================%
\subsection{Recursive solution}\label{Srt1}
%==================================================================================================%

As mentioned under \cref{Etfri0J},
the matrix product $M$ \cref{DM22} is numerically unstable
\cite[Sects.~III, IV]{StKK98}.
It is therefore preferable to solve the split boundary problem through
the recursion algorithm of Parratt \cite{Par54}.
\index{Parratt recursion!scalar}%
It is based on the insight that one does not need to compute $t_l$ and $r_l$ separately,
but only their ratio $x_l\coloneqq r_l/t_l$.
Spelling out~\cref{EcMc} with $\delta\coloneqq\delta_{l-1}$
and $s^\pm\coloneqq s^\pm_l$, we obtain
\begin{equation}\label{EParratt}
  x_{l-1} = \frac{\delta s^- + \delta s^+ x_l}{\delta^{-1}s^+ + \delta^{-1}s^- x_l}
         = \delta^{2}\frac{R+x_l}{1+Rx_l}.
\end{equation}
The second expression involves the single-interface Fresnel reflection coefficient
\begin{equation}
  R \coloneqq \frac{s^-}{s^+} = \frac{\kappa_{l-1}-\kappa_l}{\kappa_{l-1}+\kappa_l}.
\end{equation}
The recursion starts at the bottom with $x_{\nu}=0$.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Implementation}\label{SimplML}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Last updated to reflect the actual code in May 2023.

%===============================================================================
\subsection{Call chain}
%===============================================================================

{
\footnotesize
\parindent=0pt
\parskip=1ex
\flushleft

All simulations are run through the virtual function \T{ISimulation::runComputation}.

For classes \T{ScatteringSimulation} and \T{OffspecSimulation},\\
most work is done in \T{Compute::scattered\_and\_reflected},

for class \T{SpecularSimulation} in \T{Compute::reflectedIntensity},

whereas class \T{DepthprobeSimulation} performs the computation directly
in \T{runComputation}.

In function \T{Compute::scattered\_and\_reflected},\\
incoming and outgoing fluxes are obtained from functions
\T{ReSample::fluxesIn} and \T{fluxesOut},\\
and stored in instances of class \T{Fluxes},
which incarnates \T{OwningVector<IFlux>}.\\
Following that, scattering is computed by
functions \T{Compute::dwbaContribution} and
 \T{Compute::roughMultiLayerContribution}.\\
Specular intensity is added to the appropriate detector pixel by function
\T{Compute::gisasSpecularContribution}.

In \T{DepthprobeSimulation::runComputation}, incoming fluxes are obtained
from function \T{ReSample::fluxesIn}.

In functions \T{ReSample::fluxesIn} and \T{fluxesOut} call either
\T{Compute::SpecularScalar::fluxes} or \T{Compute::SpecularMagnetic::fluxes}.

For specular simulations, function \T{Compute::reflectedIntensity}
calls either \T{Compute::SpecularScalar::topLayerR} or \T{Compute::SpecularMagnetic::topLayerR}.
These functions only return amplitudes reflected from the top of the sample,
whereas the \T{fluxes} functions called for scattering or depthprobe simulation
compute up and down travelling amplitudes for each sample layer.

Functions \T{fluxes} and \T{topLayerR} are implemented
in files \SRC{Resample/Specular}{ComputeFluxScalar.cpp} and
\SRC{Resample/Specular}{ComputeFluxMagnetic.cpp},
where they share some local functions.
}

%===============================================================================
\subsection{Scalar fluxes}
%===============================================================================

The core numeric algorithm for the scalar flux computation is
implemented in \SRC{Resample/Specular}{ComputeFluxScalar.cpp}.
Here the code is simplified by omitting roughness and transmission geometry.
The code uses class \T{Spinor}, which has components \T{u} and \T{v},
here representing transmitted and reflected amplitude.
Interfaces are numbered as in~\cref{Fdefz}.

\begin{lstlisting}[language=c++, style=eclipseboxed, escapechar=|, style=footnotesize]
std::vector<Spinor>
computeTR(SliceStack& slices, std::vector<cmplx>& kz)
{
    // Parratt algorithm, pass 1:
    //    compute t/t factors and r/t ratios from bottom to top.
    size_t N = slices.size();
    std::vector<cmplx> tfactor(N-1); // transmission damping
    std::vector<cmplx> ratio(N);     // Parratt's x=r/t
    ratio[N-1] = 0;
    for (size_t i = N-1; i > 0; i--) {
        cmplx slp = 1 + kz[i]/kz[i-1];
        cmplx slm = 1 - kz[i]/kz[i-1];
        cmplx delta = exp_I(kz[i-1] * slices[i-1].thicknessOr0());
        cmplx f = delta / (slp + slm * ratio[i]);
        tfactor[i-1] = 2 * f;
        ratio[i-1] = delta * (slm + slp * ratio[i]) * f;
    }

    // Parratt algorithm, pass 2:
    //    compute r and t from top to bottom.
    std::vector<Spinor> TR(N);
    TR[0] = Spinor(1., ratio[0]);
    for (size_t i = 1; i < N; ++i) {
        TR[i].u = TR[i-1].u * tfactor[i-1]; // Spinor.u is t|\label{Lti}|
        TR[i].v = ratio[i] * TR[i].u;       // Spinor.v is r|\label{Lri}|
    }

    return TR;
}
\end{lstlisting}

The are two code blocks, each with a loop over interfaces.
The first loop runs from bottom $l=\nu$ to top $l=1$.
Variables \T{slp} and \T{slm} are the coefficients $s^\pm_l$ of~\cref{Dslpm}.
Variable \T{delta} is $\delta_{l-1}$ as defined in~\cref{Ddell}.
These are used for recursively computing transmission damping factors
\begin{equation}
  h_{l-1} \coloneqq \frac{2\delta_{l-1}}{s^+_l + s^-_l x_{l}}
\end{equation}
and Parratt ratios \cref{EParratt}
\begin{equation}
  x_{l-1} = \delta_{l-1} \frac{s^-_l+ s^+_l x_{l}}{2} h_{l-1}
         =  \delta_{l-1}^2 \frac{s^-_l+ s^+_l x_{l}}{s^+_l + s^-_l x_{l}},
\end{equation}
starting from the bottom value $x_\nu=0$.
The second loop starts from the top where $t_0=1$, $r_0=0$.
From~\cref{EcMc},
\begin{equation}
  t_{l-1} = \delta^{-1}\left(s^+ t_l + s^- r_l\right)
         = \frac{s^+ + s^- x_l}{\delta}t_l
         = h_{l-1}^{-1} t_l.
\end{equation}
Bringing $h_{l-1}$ to the other side, we obtain code~\cref{Lti}.
By definition, $x_l=r_l/t_l$.
Bringing $t_l$ to the other side, we obtain code~\cref{Lri}.

\index{Multilayer|)}%