Level Matrix Construction and Inversion
This page describes how AZURE2 constructs and inverts the level matrix (A-matrix) or channel matrix (R-matrix) to compute cross sections. The two formalisms are mathematically equivalent but differ in computational efficiency depending on the number of levels vs. channels.
Class Hierarchy
The matrix calculation is organized through an inheritance hierarchy:
GenMatrixFunc(abstract base)Defines the interface and implements the cross section calculation from T-matrix elements.
AMatrixFunc(level matrix formalism)Works with the A-matrix, a level-by-level matrix. Preferred when the number of channels exceeds the number of levels.
RMatrixFunc(channel matrix formalism)Works with the R-matrix, a channel-by-channel matrix. Preferred when the number of levels exceeds the number of channels.
Key files:
include/GenMatrixFunc.h,src/GenMatrixFunc.cppinclude/AMatrixFunc.h,src/AMatrixFunc.cppinclude/RMatrixFunc.h,src/RMatrixFunc.cppinclude/MatrixInv.h,src/MatrixInv.cpp
Computation Pipeline
For every data point, the following sequence is executed:
ClearMatrices() – allocate and zero the matrices for each \(J^\pi\) group.
FillMatrices(EPoint*) – compute the matrix elements from the R-matrix parameters.
InvertMatrices() – invert the matrix using LU decomposition.
CalculateTMatrix(EPoint*) – compute the T-matrix (transition matrix) from the inverted matrix.
CalculateCrossSection(EPoint*) – compute the cross section from the T-matrix.
A-Matrix Formalism
Matrix Construction: AMatrixFunc::FillMatrices()
The code constructs the inverse A-matrix (\(A^{-1}\)) directly, then inverts it in the next step. The matrix is indexed by level pairs \((\lambda, \lambda')\) within each \(J^\pi\) group.
The input energy in the center-of-mass frame is:
where \(E_\text{sep}\) and \(E_\text{exc}\) are the separation and excitation energies of the entrance pair.
Standard R-matrix parameterization:
where:
\(E_\lambda\) is the resonance energy of level \(\lambda\)
\(\gamma_{\lambda c}\) is the reduced width amplitude of level \(\lambda\) in channel \(c\)
\(L_c = L_0(c)\) is the shift function (Blatt-Biedenharn) evaluated at the data point energy
The sum runs over all channels in the \(J^\pi\) group
Brune parameterization (USE_BRUNE_FORMALISM):
When the Brune parameterization is enabled, additional terms modify the
channel sum for penetrability-type channels (radType == 'P'):
where:
\(B_c\) is the boundary condition for channel \(c\)
For diagonal elements (\(\lambda = \lambda'\)): \(S_c(\lambda, \lambda) = S_c^\lambda\) (the shift function at level \(\lambda\))
For off-diagonal elements (\(\lambda \neq \lambda'\)):
\[S_c(\lambda, \lambda') = \frac{S_c^\lambda (E_\text{in} - E_{\lambda'}) - S_c^{\lambda'} (E_\text{in} - E_\lambda)}{E_\lambda - E_{\lambda'}}\]
RMC formalism (USE_RMC_FORMALISM):
For radiative channels (radType == 'M' or 'E'), an imaginary term is
added to the diagonal elements:
Implementation detail: Only active levels are stored in the matrix.
An index map (level_active_index_) translates original level indices to
compact, zero-based indices for dense matrix storage. Gamma values and shift
functions are pre-cached in buffers to avoid repeated lookups. Channel
contributions with \(|\gamma| < 10^{-12}\) are skipped for efficiency.
Matrix Inversion: MatrixInv
The matrix inversion is performed by the MatrixInv class using GSL
complex LU decomposition:
The \(A^{-1}\) matrix (stored as
std::vector<std::vector<complex>>) is copied into agsl_matrix_complex.gsl_linalg_complex_LU_decomp()performs the LU factorization with partial pivoting.gsl_linalg_complex_LU_invert()computes the inverse from the LU factors.The result is copied back into a
matrix_c(std::vector<std::vector<complex>>).
This is done independently for each \(J^\pi\) group. The inversion runs sequentially (not parallelized) to avoid race conditions on the shared matrix data structures. Move semantics are used to avoid copying the result matrices.
T-Matrix Calculation: AMatrixFunc::CalculateTMatrix()
The collision matrix element \(U\) is computed from the A-matrix:
where \(P_c\) is the penetrability of channel \(c\).
The T-matrix is then:
Elastic scattering (\(c = c'\)):
\[T_{cc} = e^{2i\omega_c} - e^{i\omega_c}\,e^{i\phi_c}\,(1 + U_{cc})\]Transfer reactions (\(c \neq c'\)):
\[T_{cc'} = -e^{i\omega_c}\,e^{i\phi_{c'}}\,U_{cc'}\]
where \(\omega_c = \sigma_l\) is the Coulomb phase shift and \(\phi_c\) is the hard-sphere phase shift.
Phase factors, penetrabilities, and Coulomb amplitudes are pre-computed and
cached on the EPoint object.
R-Matrix Formalism
The R-matrix formalism (RMatrixFunc) works with the channel matrix instead.
R-matrix construction:
The [1 - RL] matrix is then constructed:
This matrix is inverted using the same GSL LU decomposition, and then multiplied by the R-matrix to obtain \([1-RL]^{-1}R\), which enters the T-matrix calculation.
When the Brune formalism is enabled, an intermediate Q-matrix inversion is performed at the level stage before building the R-matrix from the result.
Cross Section Calculation
The cross section is computed in GenMatrixFunc::CalculateCrossSection()
from the T-matrix elements. The calculation depends on the data type:
Angle-integrated cross section:
T-matrix elements from pathways with the same \((J, l, l')\) values are
coherently summed (via the TempTMatrix structure), enabling interference
between internal and external capture pathways:
where \(g\) is the geometrical factor and \(I_{12}\) is the channel spin statistical factor.
Differential cross section (with Coulomb):
When the entrance and exit channels are the same particle pair, the full Coulomb + nuclear + interference cross section is computed:
Type Definitions
The matrix types used throughout the code are defined in Constants.h:
typedef std::complex<double> complex;
typedef std::vector<std::vector<std::complex<double>>> matrix_c;
typedef std::vector<std::vector<std::vector<std::complex<double>>>> vector_matrix_c;
matrix_c– a 2D complex matrix (dense storage as nested vectors)vector_matrix_c– a collection of matrices, one per \(J^\pi\) group