%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%   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{Polarized wave propagation and scattering}\label{SPol}

In this chapter,
we generalize our treatment of wave propagation and
grazing-incidence small-angle scattering
to polarized neutrons.
\index{Neutron!polarized|(}
\index{Polarization!neutron|(}
We therefore need to study spinor wave equations,
in contrast to the scalar theory of the previous chapters.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Polarized neutrons in 2+1 dimensions}\label{Snpol0}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%===============================================================================
\subsection{Schrödinger equation for neutron spinors}\label{SnSpinor}
%===============================================================================

\index{Neutron!spin|(}%
\index{Spin|(}%

\index{Magnetizing field|see{H field}}%
\index{Magnetic field|see{B field}}%

In presence of a magnetic field,\footnote
{According to Ref.~\cite{ZaTT07},
the magnetic field is usually applied parallel to the sample surface,
but we do not rely on this.}
the propagation of free neutrons becomes spin dependent.
Therefore the scalar wavefunction of \cref{SnScalar}
must be replaced by the spinor\footnote
{Spinors can also be used to describe polarized X-rays \cite{BlKi68}.
Please let us know if there is a use case for BornAgain.}%
\index{Spinor}%
\begin{equation}\label{Dspinor}
    \Psi(\r) = \begin{pmatrix} \psi_{z+}(\r)\\\psi_{z-}(\r) \end{pmatrix}.
\end{equation}
\nomenclature[1ψ150 2r040]{$\Psi(\r)$}{Stationary coherent spinor wavefunction}%
The coupling between the neutron and the $\v{B}$~field
is given by the operator $-\gamma_\text{n}\mu_\text{nucl}\v{B}\Pauli$
with the neutron gyromagnetic factor $\gamma_\text{n}\simeq-1.91$,
the nuclear magnetron $\mu_\text{nucl}$,
and the Pauli vector $\Pauli$, composed of the three Pauli matrices
\nomenclature[1σ04]{$\Pauli$}{Pauli vector,
  composed of the three Pauli matrices: $\bm\sigma=(\sigma_x,\sigma_y,\sigma_z)$}%
\index{Pauli matrices and vector}%
(the \textit{breve} diacritic denotes
an operator in spin space, represented by a complex 2 × 2 matrix).
With the unsigned magnetic moment of the neutron,
$\mu_\text{n}\coloneqq|\gamma_\text{n}\mu_\text{nucl}|$,
\nomenclature[1μ024 2n000]{$\mu_\text{n}$}{Magnetic moment of the neutron}
\nomenclature[2h150 2r040 2t020]{$\v{B}(\r,t)$}{Magnetic induction}%
\index{B field!coupling to neutron moment}%
the Schrödinger equation~\cref{ESchrodi1}
\index{Schrodinger@Schrödinger equation!macroscopic}%
becomes
\begin{equation}\label{EHSchrodi}
  \left\{-\frac{\hbar^2}{2m}\Nabla^2+V(\r)
         +\mu_\text{n}\v{B}(\r)\Pauli-\hbar\omega\right\}
      \Psi(\r) = 0.
\end{equation}
Except near Bragg reflections, $\v{B}$ is an averaged, macroscopic field \cite{Sch78a},
just like $V$ is an averaged potential~\cref{Evrcoarse}.

We abbreviate the nuclear and the magnetic scattering-length density as
\begin{equation}
  \sldN(\r) \coloneqq v_\snuc(\r)\quad\text{and}\quad
  \sldM(\r) \coloneqq \frac{m\mu_\text{n}}{2\pi\hbar^2} B(\r),
\end{equation}
and we write $\eB$ for the unit vector in direction of the magnetic field~$\v{B}$.
\nomenclature[1μ024 00]{$\mu_0$}{Vacuum permeability, $4\pi\cdot10^{-7}$ Vs/Am}%
So the total reduced potential is given by the operator
\begin{equation}
  \OPR v(\r)
  \coloneqq  \sldN(\r)+\sldM(\r)\eB(\r)\Pauli,
\end{equation}
and we can rewrite the Schrödinger equation in analogy to~\cref{ESchrodi2} as
\index{Schrodinger@Schrödinger equation!macroscopic}%
\begin{equation}\label{ESchrodi2}
  \left\{\Nabla^2+K^2-4\pi \OPR v(\r)\right\} \Psi(\r) = 0.
\end{equation}
\index{Neutron!spin|)}%
\index{Spin|)}%

%===============================================================================
\subsection{Propagation in a multilayer}\label{Smulayer2}
%===============================================================================

In the decomposition~\cref{Edecompose_z},
both terms may become operators acting in spin space,
\begin{equation}\label{Edecompose_z2}
  \OPR v(\r) \eqqcolon \OPR\TL(z) + \delta \OPR v(\r).
\end{equation}
The unperturbed distorted wave has the form
\begin{equation}
  \Phi(\r) = \e^{i\k_\parallel\r_\parallel} \Phi(z).
\end{equation}
The horizontal wave vector~$\k_\parallel$ is constant across layers.
This motivates us to introduce the vertical vacuum wavenumber
$\kappa_0\coloneqq\sqrt{K^2-k_\parallel^2}$.
The vertical spinor wave function $\Phi(z)$ obyes the equation
\begin{equation}\label{EWaveZ2}
  \left\{\Nabla^2+\kappa_0^2-4\pi \OPR\mv(z)\right\} \Phi(z) = 0.
\end{equation}
In absence of a magnetic field, $\mv(z)$ is scalar (or proportional to the unit matrix $\OPR 1$),
and each spinor component will propagate exactly as in the scalar case of~\cref{Swave21}.
Conversely, if there is a nonzero magnetic field,
then the neutron spin will undergo Larmor precession,
which in spinor representation shows up as oscillations between the two spinor components.
In consequence, when an incident plane wave hits a magnetic medium
it becomes a superposition of two plane
waves that propagate with two different vertical wavenumbers that correspond
to the two eigenvalues of~\cref{EWaveZ2}.

We now consider a homogeneous layer with constant potential.
Similar to \cite{KeRT03,KeRT08},
we write the formal solution of~\cref{EWaveZ2} as
\begin{equation}
  \Phi(z) = \e^{-i\opkappa z\,} T + \e^{i\opkappa z} R,
\end{equation}
where $T$ and~$R$ are the transmitted and reflected spinor amplitudes.
By comparison with~\cref{EWaveZ2}, we see that the square of the operator~$\opkappa$ is
\begin{equation}\label{Dhp2}
  \opkappa^2 = \kappa_0^2-4\pi \OPR\mv = \kappa_0^2-4\pi(\sldN + \sldM \eB \Pauli).
\end{equation}

%===============================================================================
\subsection{Wavenumber operator $\opkappa$}
%===============================================================================

Without derivation,\footnote
{To verify, use standard properties of Pauli matrices.
Square~\cref{Ehp1} to reproduce~\cref{Dhp2}.
Then confirm that $\evp_\pm^2$ are eigenvalues of $\opkappa^2$.
See also \cite[\S~55, Exercice~1, p.~198]{LL3}.}
we state that the square root of $\opkappa^2$ is the operator
\begin{equation}\label{Ehp1}
  \opkappa
  = \frac{1}{2}\left[(\evp_++\evp_-) + (\evp_+-\evp_-) \eB \Pauli\right],
\end{equation}
expressed through its eigenvalues
\begin{equation}\label{Devp}
\evp_\pm
  \coloneqq \sqrt{ \kappa_0^2 - 4 \pi \sldN \pm 4 \pi \sldM }.
\end{equation}
With the abbreviations
\begin{equation}\label{Dabb}
  \alpha \coloneqq \evp_+ + \evp_-,\quad
  \beta \coloneqq \evp_+ - \evp_-,  \quad\text{and}\quad
  \v{b} \coloneqq \beta \eB,
\end{equation}
we obtain the matrix components\footnote
{Currently (jun23) implemented in function \T{MatrixFlux::computeKappa()}.}
\begin{equation}
\opkappa = \frac{1}{2}  \left( \alpha + \v{b} \Pauli  \right)
         = \frac{1}{2} \begin{pmatrix}
              \alpha + b_z & b_x - i b_y\\
              b_x + i b_y & \alpha - b_z
          \end{pmatrix}.
\end{equation}
For future reference, we note the inverse operator\footnote
{Currently (jun23) implemented in function \T{MatrixFlux::computeInverseKappa()}.}
\begin{align}\label{Ehpi}
  \opkappa^{-1}
  &= \frac{1}{2\evp_+\evp_-}
            \left[(\evp_++\evp_-) - (\evp_+-\evp_-) \eB \Pauli\right]\\
  &= \frac{2}{\alpha^2 - \beta^2} (a-\v{b}\Pauli)\\
  &= \frac{2}{\alpha^2 - \beta^2}
              \begin{pmatrix}
                 \alpha - b_z & -b_x + i b_y\\
                 -b_x - i b_y & \alpha + b_z
            \end{pmatrix} \,.
\end{align}
It does not exist if $\sldN$ is real and $\sldM= \kappa_0^2/(4\pi)-\sldN$.
If $\sldM$ is even larger, then $\opkappa$ becomes pure imaginary,
causing evanescent waves, to be discussed later~(\cref{Seva}).

%==================================================================================================%
\subsection{Eigendecomposition of $\opkappa$}\label{Skeigen}
%==================================================================================================%

To evaluate functions of the operator~$\opkappa$, we will need its eigenvalue decomposition.
We start with the matrix $\eB\Pauli$,
which has the eigenvalues $\pm 1$ and the normalized eigenspinors
\begin{equation}\label{Ev1v2}
  V_1 = \frac{1}{\sqrt{2(1+\eb_z)}} \begin{pmatrix}
 1 +\eb_z \\ \eb_x+i\eb_y \end{pmatrix},\quad
  V_2 = \frac{1}{\sqrt{2(1+\eb_z)}} \begin{pmatrix}
 \eb_x-i\eb_y\\ -1 -\eb_z \end{pmatrix}.
\end{equation}
For readability, we have omitted the subscript $\v{B}$ from the components of~$\eB$.
and the same eigenvectors as $\eB\Pauli$.
We introduce the eigenvector matrix
\begin{equation}\label{D0QofB}
  \OPR Q(\v{B}) \coloneqq \left(V_1, V_2\right)
          = \frac{1}{\sqrt{2(1+\eb_z)}}
            \begin{pmatrix} 1 +\eb_z & \eb_x-i\eb_y \\ \eb_x+i\eb_y & -1 -\eb_z \end{pmatrix}.
\end{equation}
The normalization factor becomes singular for $\eb_z=-1$.
In this case, the matrix $\eB\Pauli$ is just $\OPR\sigma_z$
and has eigenvectors $V_1=(1,0)^\dagger$ and $V_2=(0,1)^\dagger$.
Furthermore, we need to take care of the case $\v{B}=0$.
Altogether, we let
\begin{equation}\label{DQofB}
  \OPR Q(\v{B}) \coloneqq \left\{\begin{array}{ll}
          \OPR 1 &\text{ if } B=0, \\
          \OPR\sigma_x  &\text{ if } B_z=-B, \\
          (\eB+\UVEC{\v{z}})\Pauli/\sqrt{2(1+\eb_z)} &\text{ else.}
          \end{array}\right.
\end{equation}
The matrix $\opkappa$ has the eigenvalues $\evp_\pm$,
and the same eigenvectors as~$\eB\Pauli$.
Accordingly, it has the eigendecomposition
\begin{equation}\label{Ekeigen}
  \opkappa
  = \OPR Q\, \begin{pmatrix}\evp_+ & 0 \\ 0 & \evp_-\end{pmatrix} \OPR Q^\dagger,
\end{equation}
and any holomorphic function~$f(\opkappa)$ can be computed as\footnote
{Currently (jun23) implemented in function \T{MatrixFlux::eigenToMatrix}.}
\begin{equation}\label{Efkeigen}
  f(\opkappa)
  = \OPR Q \, \begin{pmatrix}f(\evp_+) & 0 \\ 0 & f(\evp_-)\end{pmatrix} \OPR Q^\dagger.
\end{equation}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Refraction and reflection at interfaces}\label{Spolif}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%===============================================================================
\subsection{Transfer matrix}
%===============================================================================

To match solutions at layer interfaces,
we use the transfer matrix method introduced in~\cref{Sacrolay}.
That section was formulated in such ways that only minimal modifications are needed now.
Instead of the vertical wave function~$\phi(z)$ and the amplitudes $t$ and~$r$,
we now have the spinors $\Phi(z)$, $T$, and~$R$.
Instead of the vertical wavenumber $\kappa\equiv k_\perp$ \cref{Dkappa},
we have the operator~$\opkappa$.
The phase factor $\delta$ \cref{Ddell} also becomes an operator,
\begin{equation}\label{Dopdel}
\OPR\delta_l\coloneqq\e^{i\opkappa_l d_l}.
\end{equation}
The equation system~\cref{EcMc} becomes
\begin{equation}\label{EcMcP}
  \begin{pmatrix}T_{l-1}\\ R_{l-1}\end{pmatrix}
  = \WW{M}_{l} \begin{pmatrix}T_l\\R_l\end{pmatrix}
\end{equation}
with the $4\times4$ transfer matrix\footnote
{Occasionally called \textit{supermatrix}
\index{Supermatrix}%
for being made of $2\times2$ submatrices \cite{Top02a}.}
\begin{equation}\label{EMDSP}
  \WW{M}_l \coloneqq \WW{D}_{l-1}\, \WW{S}_{l}
\end{equation}
 in place of~\cref{EMil}.
The phase rotation matrix~\cref{DmatD} is replaced by the block matrix
\begin{equation}\label{DDlP}
  \WW{D}_l
   \coloneqq
   \begin{pmatrix}
     \OPR\delta_l^{-1}&0\\
     0&\OPR\delta_l
   \end{pmatrix},
\end{equation}
to be discussed in the next section.
The refraction matrix~\cref{DmatS} also is replaced by a block matrix,
\begin{equation}\label{DSlP}
  \WW{S}_{l}
   \coloneqq
   \left(\begin{array}{cc}
       \OPR s^+_l & \OPR s^-_l\\
       \OPR s^-_l & \OPR s^+_l
   \end{array}\right)
\end{equation}
with the coefficients\footnote
{Currently (jun23), the matrix blocks $\OPR s^+_l$ and $\OPR s^-_l$,
possibly modified by roughness factors (see below),
are computed through local function \T{refractionMatrixBlocks}
in \SRC{Resample/Specular}{ComputeFluxMagnetic.cpp}.}
\begin{equation}\label{Dhslpm}
  \OPR s_l^\pm \coloneqq \frac{1 \pm \opkappa_{l-1}^{-1}\,\opkappa_l}{2}.
\end{equation}

%==================================================================================================%
\subsection{Phase rotation matrix}\label{Sphase}
%==================================================================================================%

With the eigendecomposition~\cref{Efkeigen},
the phase rotation matrix \cref{Dopdel} can be written\footnote
{Currently (jun23) implemented in local function \T{PhaseRotationMatrix}
in file \SRC{Resample/Flux}{MatrixFlux.cpp}.}
\begin{equation}\label{EdP2}
  \OPR\delta
  = \e^{i\opkappa d}
  = \OPR Q \, \begin{pmatrix}\e^{id\evp_+} & 0 \\ 0 & \e^{id\evp_-}\end{pmatrix} \OPR Q^\dagger.
\end{equation}
For the analysis of numerical stability,
the critical factor $\e^{i\alpha d/2}$ may be drawn in front of $\OPR Q$ in~\cref{EdP2},
\begin{equation}\label{EdP2a}
  \OPR\delta
  = \e^{i\alpha d/2}
  \OPR Q \, \begin{pmatrix}\e^{id\beta/2} & 0 \\ 0 & \e^{-id\beta/2}\end{pmatrix} \OPR Q^\dagger.
\end{equation}

%==================================================================================================%
\subsection{Interface with tanh profile}
%==================================================================================================%

\index{Tanh profile!polarized|(}%
\def\RF{\OPR{\mathcal{R}}}%

In the scalar case, the refraction matrix~\cref{DmatS}
has coefficients~\cref{Dslpm} for a sharp interface,
and modified coefficients~\cref{EslpmTanh} for a graded interface with tanh profile.
By analogy,
for polarized neutrons
the refraction matrix of a sharp interface has matrix blocks~\cref{Dhslpm},
which for a graded interface with tanh profile are replaced by
\begin{equation}\label{EhslpmTanh}
  \OPR s^\pm_a = \RF_{ab}^{-1} \pm \RF_{ab}\opkappa_{b}/\opkappa_a
\end{equation}
with the roughness factor
\begin{equation}\label{DhRbaTanh}
\RF_{ab} \coloneqq
          \sqrt{\text{tanhc}\; \pi\tau\opkappa_b}/\sqrt{\text{tanhc}\; \pi\tau\opkappa_a}
\end{equation}
that replaces the scalar factor~\cref{DRba}.
The constant~$\tau$, defined in~\cref{Dtau},
is proportional to the vertical roughness length parameter~$\sigma$.
The  eigendecomposition~\cref{Efkeigen} is applied separately
to $\opkappa_a$ and $\opkappa_b$ dependent factors,\footnote
{Currently (jun23) implemented in function \T{Compute::refractionMatrixBlocksTanh}.}
\begin{equation}\label{EhspmTanh}
  \begin{split}
  \OPR s^\pm_a \,=\,& \OPR Q_b\begin{pmatrix}1/h_b^+&0\\0&1/h_b^-\end{pmatrix}  \OPR Q_b^\dagger
          \;\OPR Q_a\begin{pmatrix}h_a^+&0\\0&h_a^-\end{pmatrix}  \OPR Q_a^\dagger
        \\&\pm\,
                  \OPR Q_b\begin{pmatrix}h_b^+c_b^+&0\\0&h_b^-c_b^-\end{pmatrix}  \OPR Q_b^\dagger
          \;\OPR Q_a\begin{pmatrix}1/(h_a^+c_a^+)&0\\0&1/(h_a^-c_a^-)\end{pmatrix}  \OPR Q_a^\dagger
   \end{split}
\end{equation}
with $h^\pm\coloneqq\sqrt{\text{tanhc}\; \pi\tau c^\pm}$.
\index{Tanh profile!polarized|)}%

%==================================================================================================%
\subsection{Névot-Croce approximation}
%==================================================================================================%

\index{N\'evot-Croce approximation!polarized|(}%
To apply the Névot-Croce approximation to polarized neutrons,
we rewrite~\cref{EslpmNC} in operator form as
\begin{equation}\label{hEslpmNC}
  \OPR s^\pm_a = \frac{1 \pm \opkappa_b/\opkappa_a}{2} \exp(-(\opkappa_b\mp\opkappa_a)^2\sigma^2/2).
\end{equation}
In contrast to the tanh roughness factor~\cref{DhRbaTanh},
the Gaussian factor here does not factorize
into separate functions of $\opkappa_{l-1}$ and $\opkappa_l$.
Therefore we need a dedicated eigendecomposition for the operator
\begin{equation}
  \opkappa^\pm \coloneqq (\opkappa_b\mp\opkappa_a)^2
       = \frac{1}{2}\left((\alpha^\pm \mp \v{b}^\pm\Pauli)^2\right)
       = \frac{1}{4}\left((\alpha^\pm)^2 + (\v{b}^\pm)^2 \mp 2 \alpha^\pm\v{b}^\pm\Pauli\right),
\end{equation}
where
\begin{equation}
  \alpha^\pm \coloneqq \alpha_b \mp \alpha_a \quad\text{ and }\quad
                   \v{b}^\pm \coloneqq \v{b}_a \mp \v{b}_b.
\end{equation}
Note that we take non-conjugated squares of complex vectors, not squared norms.
The eigendecomposition of~$\opkappa^\pm$ is computed
exactly as for the operator $\opkappa$ in~\cref{Skeigen}.\footnote
{Currently (jun23) implemented in function \T{Compute::refractionMatrixBlocksNevot}.}
\index{N\'evot-Croce approximation!polarized|)}%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Implementation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%===============================================================================
\subsection{Generalized Parratt recursion}\label{SParrattPol}
%===============================================================================

\index{Parratt recursion!polarized|(}%
We now describe the currently implemented solution
of the split boundary problem for polarized neutrons.
We start from the transfer matrix equation~\cref{EcMcP},
which we apply simultaneously to different polarization states.
To this end, we replace the spinors $T$ and $R$ by $2\times2$ matrices $\OPR t$ and $\OPR r$:
\begin{equation}\label{EotorM}
  \begin{pmatrix}\OPR t_{l-1}\\ \OPR r_{l-1}\end{pmatrix}
  = \WW{M}_{l} \begin{pmatrix}\OPR t_l\\\OPR r_l\end{pmatrix}.
\end{equation}
To generalize the Parratt recursion,
we define
\begin{equation}\label{DxP}
  \OPR x_l \coloneqq \OPR r_l \OPR {t_l}^{-1}.
\end{equation}
With \cref{EMDSP,DDlP,DSlP}, we find \cite[Eq~65]{Top02a}
\begin{align}
  \OPR x_{l-1}
   &= {\OPR\delta_{l-1}}^2 \left(\OPR s^- \OPR t + \OPR s^+ \OPR r\right)_l
                           \left(\OPR s^+ \OPR t + \OPR s^- \OPR r\right)_l^{-1}
                            \\
   &= {\OPR\delta_{l-1}}^2 \left(\OPR s^- \OPR t + \OPR s^+ \OPR r\right)_l
                           {\OPR t_l}^{-1} \OPR t_l
                            \left(\OPR s^+ \OPR t + \OPR s^- \OPR r\right)_l^{-1}
                           \\
   &= {\OPR\delta_{l-1}}^2 \left(\OPR s^- + \OPR s^+ \OPR x\right)_l
                            \left(\OPR s^+ + \OPR s^- \OPR x\right)_l^{-1},\label{ExPd2}
\end{align}
which indeed generalizes the scalar recursion~\cref{EParratt}.
\index{Parratt recursion!polarized|)}%
This recursive solution is numerically stable,
in contrast to the \textit{supermatrix
\index{Supermatrix}%
formalism} \cite{RuTD99}
that solves the split boundary value problem by inversion of~$\WW{M}$.
When modelling specular reflectivity,
then it is sufficient to compute the reflected intensity emanating from the top layer.
For incident $\OPR t_0$, the corresponding reflected matrix flux is\footnote
{Currently (jun23) implemented in function \T{Compute::polarizedReflectivity}.}
\begin{equation}
  \OPR r_0 = \OPR x_0 \OPR t_0.
\end{equation}
For a given incident spinor amplitude~$T_0$,
the reflected spinor amplitude is
\begin{equation}\label{ERxT}
  R_0 = \OPR x_0 T_0.
\end{equation}

%===============================================================================
\subsection{Fluxes inside the sample}
%===============================================================================

For modelling GISAS, we need the transmitted and reflected fluxes
in all layers of the sample.
In a first loop we compute the~$\OPR x_l$ from bottom to top as before,
and store them all in an array.
Then in a second loop we compute the $\OPR t_l$ and~$\OPR r_l$ from top to bottom.\footnote
{Currently (jun23) implemented in function \T{Compute::polarizedFluxes} and below.}

From~\cref{EotorM,DxP} we have
\begin{equation}
  \OPR t_{l-1} = {\OPR\delta_{l-1}}^{-1} \left(\OPR s^+ + \OPR s^- \OPR x\right)_l \OPR t_l.
\end{equation}
Inverting this, we obtain our recipee for computing transmitted intensities,\footnote
{Alternative expressions, involving $\OPR x_{l-1}$ rather than $\OPR s^\pm$,
can be found in \cite[Eq A.3]{KeRT03} and \cite[Eq 68]{ZaTT07}.}
\begin{equation}\label{EtPolRec}
  \OPR t_l =  {\OPR\delta_{l-1}} \left(\OPR s^+ + \OPR s^- \OPR x\right)_l^{-1} \OPR t_{l-1}
\end{equation}
The reflected intensities are then simply
\begin{equation}
  \OPR r_l = \OPR x_l \OPR t_l.
\end{equation}
For an efficient implementation, we rearrange \cref{ExPd2} from the first loop as
\begin{equation}
   \OPR x_{l-1} = {\OPR\delta_{l-1}} \left(\OPR s^- + \OPR s^+ \OPR x\right)_l \OPR F_l
\end{equation}
with matrices
\begin{equation}
  \OPR F_l \coloneqq \OPR\delta_{l-1}\left(\OPR s^+ + \OPR s^- \OPR x\right)_l^{-1},
\end{equation}
which we store to reuse them in the second loop in computing~\cref{EtPolRec},
which is just
\begin{equation}
  \OPR t_l = \OPR F_l \OPR t_{l-1}.
\end{equation}

%===============================================================================
\subsection{Numeric stability}
%===============================================================================

\iftodo

\cite{KeRT03} mention the numerical stability of this algorithm due to the strictly positive imaginary parts in the phase factors

%-------------------------------------------------------------------------------
\subsubsection{Limiting case $\kappa \rightarrow 0$}
%-------------------------------------------------------------------------------

From old document ``PolarizedImplementation'':
This case is implemented in the same way as for the scalar case, that is described in Sec.~2.2.2 of the BornAgain manual version 1.7.2.
For clarity, we briefly summarize the treatment here.
\begin{itemize}
\item One single layer: This is a trivial case, nothing needs to be calculated here as the outgoing wave is equal to the incoming one.
As a consequence, it means that we have $\UU{T_0} = \UU{1}$ and $\UU{R_0} = 0$
\item More than one layer: In that case the limit $\kappa \rightarrow 0$ is well defined.
For $\kappa = 0$, we have $\UU{R_0} = - \UU{T_0} = - \UU{1}$ and $\UU{T_j} = \UU{R_j} = 0$ for $j > 0$.
\item $\kappa = 0$ in intermediate layer: This case is not treated separately but automatically covered by the solution also present for scalar computations.
In \T{KzComputation::checkForUnderflow} a tiny imaginary part is added if the resulting value for $\kappa^2$ is getting very small.
\end{itemize}
For a single layer, the correct computation of these conditions is checked in
\T{SpecularMagneticTest::test\_degenerate}

\else...\fi

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Magnetic field}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\iftodo

From old document ``PolarizedImplementation''.\footnote
{See also issues \#654 and other issues linked there.}

%===============================================================================
\subsection{Magnetic field $z$-component conservation}
%===============================================================================

In BornAgain, the sample is assumed to be infinite along the $x$ and $y$ axes,
all parameters being constant inside each layer.
This is equivalent to the requirement of translational invariance along these axes.
On the other hand, magnetic field is known to be divergence-free,
\begin{equation*}
\nabla \cdot \boldsymbol{B} = 0 \,.
\end{equation*}
Both of these conditions result in the $z$-component of the magnetic field
(that is, the component normal to the sample surface), $B_z$,
being preserved in the whole sample and fronting medium:
\begin{equation*}
\frac{\partial B_z}{\partial z} \equiv 0 \,.
\end{equation*}

%===============================================================================
\subsection{Magnetic field in the fronting medium}
%===============================================================================

As previously described in \cite[Sec.~4.2, Fig.~13]{MaOB06},
it is reasonable to assume that the incoming beam
penetrates the fronting medium of the sample assembly from a side.
This results in $k_z$ being preserved
even when there is a non-zero magnetic field in the fronting medium.
To account for that in the calculations,
one needs to replace $k_{0z}^2$ with $k_{0z}^2 + 4 \pi \OPR\rho_{front}$,
 with $\OPR\rho_{front}$ being the SLD matrix for the fronting medium.
It is also equivalent
to subtracting the magnetic field of the fronting medium, $\boldsymbol{B}_{front}$,
from the magnetic field of each layer,
thus amending $\OPR\rho_{M}$:
\begin{equation*}
  \OPR\rho_{M}'
     = -\frac{m}{2 \pi \hbar^2} \boldsymbol{\OPR\mu}
     \left( \boldsymbol{B} - \boldsymbol{B}_{front} \right).
\end{equation*}
This amendment also concerns the nuclear (non-magnetic) scattering length density:
\begin{equation*}
\rho_n' = \rho_n - \rho_{n, front},
\end{equation*}
where $\rho_{n, front}$ is the nuclear SLD of the fronting medium.

TODO: Check this in the code

Further in the text we will omit the primes and handling of the fronting medium's properties,
however, implying that both magnetic fields and nuclear SLDs are amended in the way mentioned above.

\else...\fi

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Density operator formalism}\label{SPolDens}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%===============================================================================
\subsection{Reflected flux}
%===============================================================================

The density matrix is defined as
\index{Density operator}%
\begin{equation}
  \OPR\rho \coloneqq \displaystyle\sum_A A\, p_{\!A} A^+,
\end{equation}
where the spinors~$A$ are normalized, but not necessarily orthogonal,
and the statistical weights $p_{\!A}$ add up to~1.
When an operator~$\OPR o$ transforms the~$A$ into $\OPR o A$,
then the density operator is transformed into $\OPR o \OPR\rho \OPR o^+$.

A neutron polarizer is described an operator~$\Pi$
that shall not be further specified because
it affects observables only through a density operator
to be defined below.
Under the action of~$\OPR\Pi$,
the density matrix of the unpolarized source beam
\begin{equation}\label{Drho0}
  \OPR\rho_0 \coloneqq \begin{pmatrix}1/2&0\\0&1/2\end{pmatrix} = \frac{\OPR 1}{2}
\end{equation}
becomes transformed into the density matrix of the polarized incident beam
\begin{equation}\label{Erho1}
  \OPR\rho_1
  = \OPR\Pi_\si\,\OPR\rho_0\,\OPR\Pi^+_\si
  = \OPR\Pi_\si\,\OPR\Pi^+_\si\,\OPR\rho_0
  \equiv \OPR\rho_\si\,\OPR\rho_0.
\end{equation}
In the second equality we used the fact that $\OPR\rho_0$ is proportional to the
unity matrix and therefore commutes with any other matrix.
The product of $\OPR\Pi_\si$ and its conjugate transpose are then combined
into the polarizer density operator
\begin{equation}\label{Drhoi}
  \OPR\rho_\si
  \coloneqq \OPR\Pi_\si\,\OPR\Pi^+_\si.
\end{equation}


In~\cref{ERxT} we found for a given incident spinor amplitude $T$
a reflected spinor amplitude $\OPR x_0 T$,
where $\OPR x_0$ is a matrix obtained from the generalized Parratt recursion.
Accordingly, the density matrix of the incident beam is transformed into the
density matrix of the reflected beam
\begin{equation}
  \rho_2 = \OPR x_0\, \OPR\rho_1\, \OPR x^+.
\end{equation}
Finally, the beam is passed through a polarization analyzer and
the density matrix becomes%
\begin{equation}
  \rho_3 = \OPR\Pi_\sf\, \OPR\rho_2\, \OPR\Pi^+_\sf.
\end{equation}
At this point, the flux is given by the trace
\index{Flux!polarized}%
\begin{equation}
  I_3 = \Tr \OPR\rho_3
  = \Tr \OPR\Pi_\sf\, \OPR\rho_2\, \OPR\Pi^+_\sf
  = \Tr \OPR\Pi^+_\sf\,\OPR\Pi_\sf\, \OPR\rho_2
  \equiv \Tr \OPR\rho_\sf\,\OPR\rho_2.
\end{equation}
In the second equality we used
the invariance of a trace under rotation of matrix factors.
In the final identity,
we introduced the polarizer density operator
\begin{equation}\label{Drhof}
  \OPR\rho_\sf \coloneqq \OPR\Pi^+_\sf\,\OPR\Pi_\sf.
\end{equation}
Collecting everything, we obtain\footnote
{The leading factor $1/2$, which comes from the density matrix \cref{Drho0}
of the unpolarized source beam, is ignored in BornAgain.
Besides that, \cref{EI3} is currently (June 2023) implemented in
function \T{Compute::magneticR} in file \SRC{Sim/Computation}{SpecularComputation.cpp}.}
\begin{equation}\label{EI3}
  I_3 = \frac{1}{2} \Tr \OPR\rho_\sf\, \OPR x_0\, \OPR\rho_\si\, \OPR x_0^+.
\end{equation}

%===============================================================================
\subsection{Parameterization of the polarizer density operator}
%===============================================================================

As any other 2$\times$2 matrix,
the polarization operator
\index{Polarization!density operator}%
can be written as
\begin{equation}\label{EdecomposePi}
  \OPR\Pi = p_0\OPR 1 + \v{p}\Pauli,
\end{equation}
and the polarizer density operator
\index{Polarization!density operator}%
as
\begin{equation}\label{Edecompose1rho}
  \OPR\rho = r_0\OPR 1 + \v{r}\Pauli.
\end{equation}
From \cref{Drhoi} or~\cref{Drhof}, we know that $\OPR\rho=\OPR\Pi\,\OPR\Pi^+$.
Inserting \cref{EdecomposePi}, we can conclude that $\OPR\rho$ is Hermitean,
that $r_0$ and~$\v{r}$ are real,
and that $|\v{r}|\le|r_0|$.
This allow up to replace \cref{Edecompose1rho} by
\begin{equation}\label{Edecompose2rho}
  \OPR\rho = \left(\OPR 1 + \v{P}\Pauli\right) \tau.
\end{equation}
We identify $\v{P}$ as the \textit{polarization vector},
\index{Polarization!vector}%
and $\tau$ as the \textit{mean transmission} of an unpolarized beam;
it can take values between 0 and~1/2,
whereas the polarization strength $P\coloneqq|\v{P}|$ may take
values between 0 and~1.
For a source flux~$I_0$, the flux after a beam polarizer has the components
\begin{equation}
  I_\pm \coloneqq \Tr (\pm \v{\hat P} \Pauli) \OPR\rho_\si \OPR\rho_0 I_0
   = \frac{1}{2}(1\pm P)\tau I_0.
\end{equation}
The polarization ratio is
\begin{equation}
  \frac{I_+ - I_-}{I_+ + I_-}
  = P
\end{equation}
in accord with the conventional definition of the polarization degree~$P$
\cite{Wil88}.

\index{Neutron!polarized|)}
\index{Polarization!neutron|)}