Detailed Description

Functions for atomic ionisation form factors.

Bethe-ridge (plane-wave) completion of the truncated multipole sum.

Classes
struct	FormFactorsRPA
	Bare and RPA form factors of every core orbital, from calculate_formFactors_RPA() More...

struct	IonisationChannel
	One ionisation channel: a hole in a core orbital, and the ejected (continuum) state. More...

struct	IonisedOrbital
	A core orbital ionised at one energy transfer, with the continuum states of its ejected electron. More...

struct	MomentumOrbital
	Momentum-space (Fourier transformed) bound orbital, for the all-K plane-wave response. More...

struct	RPAOptions
	Options for the RPA (core polarisation) form factors. More...

Typedefs
using	FormFactorSet = std::array< LinAlg::Matrix< double >, 13 >
	The 13 form factors of one bound orbital (or their sum), in the fixed order: {V_T, V_E, V_M, V_L, X, A_T, A_E, A_M, A_L, Y, Z, S, P}. A factor that is not calculated is left empty (0x0).

using	ChannelAmplitudes = std::array< std::complex< double >, 10 >
	Reduced matrix elements <e\|\|h\|\|a> of one (bound, continuum) channel for each operator of the multipole set, in the order of multipole_operators(): {t, E, M, L, t5, E5, M5, L5, S, S5}. Real (imaginary part zero) without RPA; complex (outgoing-wave amplitudes) with RPA. Zero if not calculated.

Enumerations
enum class	Coupling { Vector , Scalar , AxialVector , PseudoScalar , Error }
	DM-electron couplings. More...

enum class	OutputFormat { matrix , xyz , Error }
	Format for output file. (All new code should use xyz; matrix kept for legacy code) More...

enum class	Units { Atomic , Particle , Error }
	Units used in output file. More...

enum class	AtomicMethod { HF , Zeff , ZeffAnalytic , RPA }
	Method used to solve bound/continuum states for form factors. More...

Functions
AtomicMethod	parseStatesMethod (const std::string &in_method)
	Parses string (HF, Zeff, ZeffAnalytic, RPA) to AtomicMethod (case-insensitive).

std::string	parseStatesMethod (const AtomicMethod &in_method)
	AtomicMethod to string (HF, Zeff, ZeffAnalytic, RPA)

LinAlg::Matrix< double >	calculateK_nk (const HF::HartreeFock vHF, const DiracSpinor &Fnk, int max_L, const Grid &Egrid, const DiracOperator::jL jl, bool force_rescale, bool hole_particle, bool force_orthog, bool zeff_cont, bool zeff_bound, double ec_cut=1.0e99)
	Calculates the ionisation factor K(E,q) for one core state, using the standard (single multipole operator) method. New code should use calculate_formFactors.

FormFactorSet	allocate_formFactors (std::size_t E_steps, std::size_t q_steps, bool vectorQ, bool axialQ, bool scalarQ, bool pseudoscalarQ, bool spatialQ)
	Small helper to allocate/size the requested factors.

std::array< std::unique_ptr< DiracOperator::TensorOperator >, 10 >	multipole_operators (const Grid &grid, bool low_q, const SphericalBessel::JL_table *jK_tab, bool vectorQ, bool axialQ, bool scalarQ, bool pseudoscalarQ, bool spatialQ)
	Constructs the required multipole operator set for the form factors (at w=q=0)

void	accumulate_formFactors (FormFactorSet *K_factors, std::size_t iE, std::size_t iq, double tkp1_x, const ChannelAmplitudes &A)
	Adds the contribution of one (bound, continuum) channel to the form factors at (iE, iq).

std::vector< DiracSpinor >	model_bound_states (const HF::HartreeFock &vHF, AtomicMethod method)
	The bound state of each core orbital as used in the matrix elements: the real (HF) orbital, or its H-like (Zeff) version.

void	accumulate_multipole_sum (FormFactorSet *K_factors, std::size_t iE, const DiracSpinor &Fa, double occ_frac, const std::vector< DiracSpinor > &continuum, std::array< std::unique_ptr< DiracOperator::TensorOperator >, 10 > &multipoles, int Kmin, int Kmax, const std::vector< std::pair< std::size_t, double > > &q_columns)
	Adds the multipole sum of one bound orbital at one energy transfer to the form factors: every rank K, momentum transfer, and continuum state.

void	add_formFactors (std::vector< FormFactorSet > *K_nk, const std::vector< FormFactorSet > &dK)
	Adds the factors of each orbital of dK to those of K_nk (same orbitals, same allocation pattern; factors not allocated in either are skipped)

std::vector< FormFactorSet >	calculate_formFactors (const HF::HartreeFock *vHF, const std::vector< DiracSpinor > &bound_states, const std::optional< std::array< int, 2 > > &lc_minmax, double ec_min, double ec_max, bool force_rescale, bool hole_particle, bool force_orthog, const std::vector< double > &Egrid, const std::vector< double > &qgrid, bool diagonal_Eq, bool low_q, const SphericalBessel::JL_table &jK_tab, int Kmin, int Kmax, bool vectorQ, bool axialQ, bool scalarQ, bool pseudoscalarQ, bool spatialQ, AtomicMethod method=AtomicMethod::HF)
	Calculates all 13 form factors (V, A, S, P) for every core orbital.

std::pair< int, int >	continuum_l_range (const DiracSpinor &Fa, int Kmax, const std::optional< std::array< int, 2 > > &lc_minmax)
	Range of continuum l reached from a bound orbital by multipoles of rank up to Kmax.

std::vector< std::vector< IonisedOrbital > >	solve_ionised_orbitals (const HF::HartreeFock *vHF, const std::vector< double > &omegas, double ec_min, double ec_max, int Kmax, const std::optional< std::array< int, 2 > > &lc_minmax, bool force_rescale, bool hole_particle, bool force_orthog)
	The core orbitals ionised by each of a set of energy transfers, each with the Hartree-Fock continuum states of its ejected electron.

std::vector< IonisationChannel >	construct_channels (const std::vector< DiracSpinor > &core, const std::vector< IonisedOrbital > &ionised)
	Every (hole orbital, ejected state) channel at one energy transfer.

std::vector< std::size_t >	active_operators (const std::array< std::unique_ptr< DiracOperator::TensorOperator >, 10 > &operators)
	Indices of the requested (non-null) operators of a multipole_operators() set, in its order.

std::optional< double >	solve_channel_amplitudes (const DiracOperator::TensorOperator &h, std::size_t i_op, ExternalField::TDHFcntm rpa, double omega, const RPAOptions &rpa_options, const std::vector< IonisationChannel > &channels, std::vector< ChannelAmplitudes > A_bare, std::vector< ChannelAmplitudes > *A_rpa, bool print=false)
	One RPA solve: the bare and RPA amplitudes of every channel, for one operator at one (E, K, q).

FormFactorsRPA	solve_formFactors_RPA (const HF::HartreeFock *vHF, const std::optional< std::array< int, 2 > > &lc_minmax, double ec_min, double ec_max, bool force_rescale, bool hole_particle, bool force_orthog, const std::vector< double > &Egrid, const std::vector< double > &qgrid, bool diagonal_Eq, bool low_q, const SphericalBessel::JL_table &jK_tab, int Kmin, int Kmax, bool vectorQ, bool axialQ, bool scalarQ, bool pseudoscalarQ, bool spatialQ, const RPAOptions &rpa_options={})
	Bare and RPA (core polarisation) form factors of every core orbital, with the RPA solved at every point of the grids, from the outgoing-wave TDHF.

FormFactorsRPA	calculate_formFactors_RPA (const HF::HartreeFock *vHF, const std::optional< std::array< int, 2 > > &lc_minmax, double ec_min, double ec_max, bool force_rescale, bool hole_particle, bool force_orthog, const std::vector< double > &Egrid, const std::vector< double > &qgrid, bool diagonal_Eq, bool low_q, const SphericalBessel::JL_table &jK_tab, int Kmin, int Kmax, bool vectorQ, bool axialQ, bool scalarQ, bool pseudoscalarQ, bool spatialQ, const RPAOptions &rpa_options={})
	Calculates the form factors for every core orbital, without and with RPA (core polarisation), the RPA solved only within the limits of RPAOptions.

std::pair< std::size_t, std::size_t >	count_failed_rpa (const LinAlg::Matrix< double > &eps, double eps_fail)
	Counts the points of an RPA eps table at which the solve failed, and those of them with a converged neighbour along the columns.

void	interpolate_failed_rpa (const LinAlg::Matrix< double > &eps, double eps_fail, const LinAlg::Matrix< double > &K_bare, LinAlg::Matrix< double > *K_rpa)
	Corrects the points at which the RPA solve failed, by interpolating the relative RPA shift from the neighbouring columns.

void	interpolate_failed_rpa_2d (const LinAlg::Matrix< double > &eps, double eps_fail, const LinAlg::Matrix< double > &K_bare, LinAlg::Matrix< double > *K_rpa)
	Corrects the (E, q) points at which the RPA solve failed, by interpolating the relative RPA shift from the neighbours in E and q.

std::vector< std::size_t >	factor_operators (std::size_t i_factor)
	The operators (indices into ChannelAmplitudes) that a form factor is built from.

bool	check_radial_grid_q (double qmax, const Grid &rgrid)
	The q part of check_radial_grid only: grid spacing near r ~ a0 must resolve exp(iq.r) at qmax (au). Returns false (with advice) if not.

bool	check_radial_grid_E (double Emax, const Grid &rgrid, double alpha=PhysConst::alpha)
	The E part of check_radial_grid only: grid must resolve the continuum oscillations out to rmax at Emax (au). Returns false (with advice) if not.

bool	check_radial_grid (double Emax, double qmax, const Grid &rgrid, double alpha=PhysConst::alpha)
	Checks the radial grid is dense enough for the continuum states and momentum transfers required.

void	write_to_file_xyz (const std::string &filename, const std::vector< double > &E_grid, const std::vector< double > &q_grid, const std::vector< std::string > &titles, const std::vector< std::string > &descriptions, std::vector< LinAlg::Matrix_view< const double > > factors, Units units=Units::Particle, int num_digits=6, bool diagonal=false)
	Writes output file in 'xyz' form: for easy 2D interpolation.

void	write_to_file_xyz_13 (const std::string &filename, const std::vector< double > &E_grid, const std::vector< double > &q_grid, const std::vector< std::string > &titles, const std::vector< std::string > &descriptions, const FormFactorSet &K_factors, Units units=Units::Particle, int num_digits=6, bool diagonal=false)
	As write_to_file_xyz, for a FormFactorSet (empty factors are skipped)

void	write_to_file_matrix (const LinAlg::Matrix< double > &K, const std::vector< double > &E_grid, const std::vector< double > &q_grid, const std::string &filename, int num_digits=5, Units units=Units::Particle)
	Writes output file in matrix form.

double	Zeff_nonrel (double en, int n)
	Effective charge of an orbital, from its binding energy.

bool	rpa_failed (double eps, double eps_fail)
	A failed RPA solve: final eps above eps_fail, or nan.

MomentumOrbital	momentum_orbital (const DiracSpinor &Fa, double en, double p_min=1.0e-3, double p_max=1.0e4, std::size_t points_per_decade=150, double density_cut=1.0e-12)
	Fourier transforms a bound orbital onto a logarithmic momentum grid (MomentumOrbital).

std::array< double, 11 >	planewave_traces (double p, double F, double G, int kappa, double pf, double q, double ef, double alpha)
	Integrand of the all-K plane-wave response at bound-electron momentum p: the Dirac traces of every factor.

std::array< double, 13 >	planewave_formFactors (const MomentumOrbital &orb, double E, double q, double alpha)
	The 13 form factors of one orbital with a free (plane-wave) ejected electron, summed over all multipoles, at energy transfer E and momentum transfer q.

double	captured_fraction (const DiracSpinor &Fa, int Kmax, std::size_t iq, const SphericalBessel::JL_table &jK_tab)
	Fraction of the norm of \( e^{i\vb{q}\cdot\vb{r}}\psi_a \) carried by the multipoles K <= Kmax.

std::vector< FormFactorSet >	calculate_ridge_correction (const HF::HartreeFock *vHF, const std::vector< DiracSpinor > &bound_states, double ec_min, double ec_max, const std::vector< double > &Egrid, const std::vector< double > &qgrid, const SphericalBessel::JL_table &jK_tab, int Kmax, bool vectorQ, bool axialQ, bool scalarQ, bool pseudoscalarQ, bool spatialQ, double ridge_eps)
	The Bethe-ridge correction: completes the multipole sum of every core orbital to all K with plane waves.

Class Documentation

◆ Kion::FormFactorsRPA

struct Kion::FormFactorsRPA

Class Members
vector< FormFactorSet >	bare {}	Bare (Hartree-Fock) factors of each core orbital, indexed as the core.
vector< FormFactorSet >	rpa {}	Factors including RPA (the total, not the correction), indexed as core.
Matrix< double >	eps {}	Worst RPA eps over K and operators at each (E, q); zero where no orbital is ionised, or where the RPA was not solved (outside the RPAOptions limits)
size_t	n_solves {0}	Number of RPA solves made (one per (E, K, operator, q) with a non-zero amplitude, within the limits)
size_t	n_failed {0}	How many of those did not converge (eps above RPAOptions::eps_fail, or nan); zero for a first-order solve, whose eps is not tested.
size_t	n_interpolated {0}	How many of the failed solves had a converged neighbour (same K and operator, adjacent E or q) to interpolate the RPA shift from.

◆ Kion::IonisationChannel

struct Kion::IonisationChannel

Class Members
size_t	hole_index {}	Index of the hole orbital in the core.
const DiracSpinor *	hole {nullptr}	The hole (core) orbital.
const DiracSpinor *	ejected {nullptr}	The ejected (continuum) state.

◆ Kion::IonisedOrbital

struct Kion::IonisedOrbital

Class Members
size_t	core_index {}	Index of the orbital in the core.
ContinuumOrbitals	ejected	Continuum states of the ejected electron.

◆ Kion::MomentumOrbital

struct Kion::MomentumOrbital

Class Members
int	kappa {}	Dirac quantum number of the orbital.
double	en {}	Binding energy (au); that of the real (HF) orbital.
double	num_electrons {}	Number of electrons in the orbital: (2j+1) times the occupation fraction.
vector< double >	p {}	Momentum grid (au), logarithmic.
vector< double >	F {}	Transformed large component, f~_a, at each p.
vector< double >	G {}	Transformed small component, g~_a, at each p.

◆ Kion::RPAOptions

struct Kion::RPAOptions

Class Members
int	max_its {60}	Maximum RPA iterations per solve; 1 gives the first-order correction.
double	eps {1.0e-10}	RPA convergence target.
double	eps_fail {1.0e-3}	An RPA solve whose final eps is above this (or nan) is discarded: the bare (no-RPA) amplitude is used for that (E, q, K, operator)
double	E_max {1.0e99}	The RPA is solved only for energy transfers E up to this (au); above it the factors are the bare ones (see calculate_formFactors_RPA)
double	q_max {1.0e99}	The RPA is solved only for momentum transfers q up to this (au); in the diagonal (q = E/c) case this is a limit on E.
int	K_max {999}	The RPA is solved only for multipoles of rank K up to this; the higher multipoles contribute their bare factors.

Typedef Documentation

◆ FormFactorSet

using Kion::FormFactorSet = typedef std::array<LinAlg::Matrix<double>, 13>

The 13 form factors of one bound orbital (or their sum), in the fixed order: {V_T, V_E, V_M, V_L, X, A_T, A_E, A_M, A_L, Y, Z, S, P}. A factor that is not calculated is left empty (0x0).

◆ ChannelAmplitudes

using Kion::ChannelAmplitudes = typedef std::array<std::complex<double>, 10>

Reduced matrix elements <e||h||a> of one (bound, continuum) channel for each operator of the multipole set, in the order of multipole_operators(): {t, E, M, L, t5, E5, M5, L5, S, S5}. Real (imaginary part zero) without RPA; complex (outgoing-wave amplitudes) with RPA. Zero if not calculated.

Enumeration Type Documentation

◆ Coupling

enum class Kion::Coupling

strong

DM-electron couplings.

◆ OutputFormat

enum class Kion::OutputFormat

strong

Format for output file. (All new code should use xyz; matrix kept for legacy code)

xyz : For easy 2D interpolation. List formatted with each row 'E q K(E,q)'
matrix : Outputs entire matrix in table form. E and q grids printed prior.

◆ Units

enum class Kion::Units

strong

Units used in output file.

Atomic : [q] = [1/a_0], [E] = Hartree
Particle : [q] = eV, [E] = eV
Form factors are defined to be dimensionless

◆ AtomicMethod

enum class Kion::AtomicMethod

strong

Method used to solve bound/continuum states for form factors.

HF : Real (Hartree-Fock) bound and continuum states (standard method).
Zeff : H-like (Zeff) bound and continuum states, solved numerically with DiracODE.
ZeffAnalytic : H-like (Zeff) bound and continuum states, using exact analytic Dirac-Coulomb functions. Relativistic continuum requires FLINT (see DiracContinuum::available).
RPA : Hartree-Fock states, with core-polarisation (RPA) corrections to every amplitude, from the outgoing-wave TDHF. A method for the amplitudes, not the states: the states are those of HF (see calculate_formFactors_RPA).

Function Documentation

◆ parseStatesMethod() [1/2]

AtomicMethod Kion::parseStatesMethod ( const std::string & in_method )

Parses string (HF, Zeff, ZeffAnalytic, RPA) to AtomicMethod (case-insensitive).

Parameters

in_method Method name; unknown input warns and defaults to HF.

Returns: The corresponding AtomicMethod.

◆ parseStatesMethod() [2/2]

std::string Kion::parseStatesMethod ( const AtomicMethod & in_method )

AtomicMethod to string (HF, Zeff, ZeffAnalytic, RPA)

◆ calculateK_nk()

LinAlg::Matrix< double > Kion::calculateK_nk	(	const HF::HartreeFock *	vHF,
		const DiracSpinor &	Fnk,
		int	max_L,
		const Grid &	Egrid,
		const DiracOperator::jL *	jl,
		bool	force_rescale,
		bool	hole_particle,
		bool	force_orthog,
		bool	zeff_cont,
		bool	zeff_bound,
		double	ec_cut = `1.0e99`
	)

Calculates the ionisation factor K(E,q) for one core state, using the standard (single multipole operator) method. New code should use calculate_formFactors.

\[ K(E,q) = \sum_{L,e} (2L+1)\,x_{\rm occ}\,|\redmatel{e}{j_L}{a}|^2 \]

summed over the multipoles L up to max_L and the continuum states e with \( l_e \) within max_L of \( l_a \). The continuum energy is \( \en_c = E + \en_a \); grid points at which this is not positive (or exceeds ec_cut) are left zero. Parallelised over E or q, whichever is larger.

Note: Should be equivilant to temporal component of calculate_formFactors. Prefer calculate_formFactors() or calculate_formFactors_RPA() for new code.

Parameters

vHF	Hartree-Fock potential.
Fnk	Bound (core) orbital being ionised.
max_L	Maximum multipolarity L.
Egrid	Energy transfer grid, in au.
jl	Operator providing the reduced matrix elements; its q grid sets the columns.
force_rescale	Rescale V(r) at large r for the continuum states.
hole_particle	Solve the continuum in the V^(N-1) potential.
force_orthog	Orthogonalise the continuum states to the core.
zeff_cont	Use H-like (Zeff) continuum states.
zeff_bound	Use an H-like (Zeff) bound state in the matrix elements (the HF energy is still used).
ec_cut	Maximum continuum energy, in au.

Returns: K(E,q) as a matrix: each row a new E, each column a new q.

◆ allocate_formFactors()

FormFactorSet Kion::allocate_formFactors	(	std::size_t	E_steps,
		std::size_t	q_steps,
		bool	vectorQ,
		bool	axialQ,
		bool	scalarQ,
		bool	pseudoscalarQ,
		bool	spatialQ
	)

Small helper to allocate/size the requested factors.

Each requested factor is allocated (E_steps x q_steps) and zeroed; the others are left empty (0x0), and are skipped by every function that takes a FormFactorSet. The interference terms X and Y follow their vector and axial parts (spatial only); Z requires vector, axial, and spatial.

Parameters

E_steps	Number of energy-transfer grid points (rows).
q_steps	Number of momentum-transfer grid points (columns).
vectorQ	Include the vector factors.
axialQ	Include the axial-vector factors.
scalarQ	Include the scalar factor.
pseudoscalarQ	Include the pseudoscalar factor.
spatialQ	Include the spatial (E, M, L) components; if false, only the temporal components are allocated.

Returns: The allocated set, in the fixed FormFactorSet order.

◆ multipole_operators()

std::array< std::unique_ptr< DiracOperator::TensorOperator >, 10 > Kion::multipole_operators	(	const Grid &	grid,
		bool	low_q,
		const SphericalBessel::JL_table *	jK_tab,
		bool	vectorQ,
		bool	axialQ,
		bool	scalarQ,
		bool	pseudoscalarQ,
		bool	spatialQ
	)

Constructs the required multipole operator set for the form factors (at w=q=0)

Returned in the fixed order {t, E, M, L, t5, E5, M5, L5, S, S5}, with nullptr for those not requested: vector temporal t, electric E, magnetic M, longitudinal L; axial, their \( \gamma^5 \) partners; scalar S; pseudoscalar S5.

Each is constructed at rank 0 and zero frequency: updateRank() then updateFrequency() must be called before use (see DiracOperator::MultipoleOperator for the units of the frequency: qc, or E in the diagonal case).

Parameters

grid	Radial grid on which the operators act.
low_q	Use the low-momentum (long-wavelength) form.
jK_tab	Precomputed spherical Bessel table (may be nullptr).
vectorQ	Build the vector operators.
axialQ	Build the axial-vector operators.
scalarQ	Build the scalar operator.
pseudoscalarQ	Build the pseudoscalar operator.
spatialQ	Build the spatial (E, M, L) components.

Returns: The operator set; entries not requested are nullptr.

◆ accumulate_formFactors()

void Kion::accumulate_formFactors	(	FormFactorSet *	K_factors,
		std::size_t	iE,
		std::size_t	iq,
		double	tkp1_x,
		const ChannelAmplitudes &	A
	)

Adds the contribution of one (bound, continuum) channel to the form factors at (iE, iq).

With weight \( w = (2K+1)x_{\rm occ} \) (tkp1_x), the squared terms get \( w|A|^2 \) and the interference terms \( w\,{\rm Re}(A\,A'^*) \): X from (t, L), Y from (t5, L5), and Z from (E5, M) - (E, M5). For real (bare) amplitudes this is the plain product.

Parameters

K_factors	Factors to add to; empty (not allocated) ones are skipped.
iE	Energy-transfer grid index.
iq	Momentum-transfer grid index.
tkp1_x	Weight (2K+1) times the occupation fraction.
A	Channel amplitudes, in multipole_operators() order.

◆ model_bound_states()

std::vector< DiracSpinor > Kion::model_bound_states	(	const HF::HartreeFock &	vHF,
		AtomicMethod	method
	)

The bound state of each core orbital as used in the matrix elements: the real (HF) orbital, or its H-like (Zeff) version.

For AtomicMethod::Zeff the H-like state is solved numerically (DiracODE) in the pointlike -Zeff/r potential, for AtomicMethod::ZeffAnalytic it is the exact Dirac-Coulomb function; otherwise the HF orbital itself. Zeff is that of each orbital from its binding energy (Zeff_nonrel(), as DarkARC). Every state carries the occupation of the real orbital. The Zeff states have the H-like energy, not the HF one: continuum energies must be taken from the real orbital.

Parameters

vHF	Hartree-Fock; its core defines the orbitals.
method	States method (see AtomicMethod).

Returns: One bound state per core orbital, indexed as the core.

◆ accumulate_multipole_sum()

void Kion::accumulate_multipole_sum	(	FormFactorSet *	K_factors,
		std::size_t	iE,
		const DiracSpinor &	Fa,
		double	occ_frac,
		const std::vector< DiracSpinor > &	continuum,
		std::array< std::unique_ptr< DiracOperator::TensorOperator >, 10 > &	multipoles,
		int	Kmin,
		int	Kmax,
		const std::vector< std::pair< std::size_t, double > > &	q_columns
	)

Adds the multipole sum of one bound orbital at one energy transfer to the form factors: every rank K, momentum transfer, and continuum state.

For each K in [Kmin, Kmax] and each (column, qc) of q_columns, sets the rank and frequency of every (non-null) operator, forms the channel amplitudes \( \redmatel{e}{h}{a} \) with each continuum state e, and accumulates them with weight \( (2K+1)x_{\rm occ} \) (see accumulate_formFactors()) into row iE of K_factors. The operators take qc (or E itself in the diagonal, massless-absorption, case) as their frequency; see DiracOperator::MultipoleOperator.

Parameters

K_factors	Factors to add to; empty (not allocated) ones are skipped.
iE	Row (energy-transfer index) of `K_factors` to add to.
Fa	Bound orbital (as used in the matrix elements).
occ_frac	Its occupation fraction.
continuum	Continuum states of the ejected electron (all l, kappa).
multipoles	Operator set (multipole_operators()); modified in place (rank and frequency), so each thread needs its own.
Kmin	Minimum multipolarity K.
Kmax	Maximum multipolarity K.
q_columns	The (column index, frequency qc) pairs to evaluate.

◆ add_formFactors()

void Kion::add_formFactors	(	std::vector< FormFactorSet > *	K_nk,
		const std::vector< FormFactorSet > &	dK
	)

Adds the factors of each orbital of dK to those of K_nk (same orbitals, same allocation pattern; factors not allocated in either are skipped)

◆ calculate_formFactors()

std::vector< FormFactorSet > Kion::calculate_formFactors	(	const HF::HartreeFock *	vHF,
		const std::vector< DiracSpinor > &	bound_states,
		const std::optional< std::array< int, 2 > > &	lc_minmax,
		double	ec_min,
		double	ec_max,
		bool	force_rescale,
		bool	hole_particle,
		bool	force_orthog,
		const std::vector< double > &	Egrid,
		const std::vector< double > &	qgrid,
		bool	diagonal_Eq,
		bool	low_q,
		const SphericalBessel::JL_table &	jK_tab,
		int	Kmin,
		int	Kmax,
		bool	vectorQ,
		bool	axialQ,
		bool	scalarQ,
		bool	pseudoscalarQ,
		bool	spatialQ,
		AtomicMethod	method = `AtomicMethod::HF`
	)

Calculates all 13 form factors (V, A, S, P) for every core orbital.

The continuum states of each orbital are solved at every E for which the ejected electron energy \( \en_c = E + \en_a \) lies in (ec_min, ec_max]. The continuum l are those reached from the orbital by multipoles of rank up to Kmax (both parities), \( j_e = j_a \pm K_{\rm max} \), \( l_e = j_e \pm 1/2 \), clipped to lc_minmax if given. Parallel over the (orbital, E) pairs.

The bound state used in the matrix elements of each orbital is bound_states[ia]: the Hartree-Fock orbital itself, or (for the H-like methods) its Zeff version, from model_bound_states(). With method Zeff or ZeffAnalytic the continuum is also H-like, in the pointlike potential of the same charge as that model state (Zeff_nonrel() of the HF binding energy), solved numerically (DiracODE) or with the exact analytic Dirac-Coulomb functions. The occupation x_i is that of the bound state used; continuum energies always use the real (HF) orbital.

Parameters

vHF	Hartree-Fock potential; its core defines the orbitals.
bound_states	Bound state of each core orbital as used in the matrix elements, indexed as the core.
lc_minmax	Optional limits on the continuum orbital l.
ec_min	Minimum ejected electron energy, in au.
ec_max	Maximum ejected electron energy, in au.
force_rescale	Rescale V(r) at large r for the continuum states.
hole_particle	Solve the continuum in the V^(N-1) potential of the hole (include the hole-particle interaction).
force_orthog	Orthogonalise the continuum states to the core.
Egrid	Energy transfer grid, in au.
qgrid	Momentum transfer grid, in au (size 1 if diagonal).
diagonal_Eq	Momentum transfer set equal to the energy transfer (absorption of a massless particle).
low_q	Use the low-q form of the operators.
jK_tab	Precomputed spherical Bessel table.
Kmin	Minimum multipolarity K.
Kmax	Maximum multipolarity K.
vectorQ	Calculate the vector factors.
axialQ	Calculate the axial-vector factors.
scalarQ	Calculate the scalar factor.
pseudoscalarQ	Calculate the pseudoscalar factor.
spatialQ	Calculate the spatial (E, M, L) components.
method	Continuum states: AtomicMethod::HF : Hartree-Fock AtomicMethod::Zeff : H-like, DiracODE AtomicMethod::ZeffAnalytic : H-like, analytic (AtomicMethod::RPA is not a states method, and is treated as HF here; see calculate_formFactors_RPA.)

Returns: One FormFactorSet per core orbital, indexed as the core (zero for an orbital that is not ionised anywhere on the E grid). See allocate_formFactors() for which factors are calculated.

Note: force_rescale and hole_particle have no effect for Zeff states.

◆ continuum_l_range()

std::pair< int, int > Kion::continuum_l_range	(	const DiracSpinor &	Fa,
		int	Kmax,
		const std::optional< std::array< int, 2 > > &	lc_minmax
	)

Range of continuum l reached from a bound orbital by multipoles of rank up to Kmax.

Either parity: \( j_e = j_a \pm K_{\rm max} \), \( l_e = j_e \pm 1/2 \), clipped to lc_minmax if given.

Parameters

Fa	Bound (core) orbital.
Kmax	Maximum multipolarity.
lc_minmax	Optional limits {lc_min, lc_max} on the continuum l.

Returns: {lc_min, lc_max}.

◆ solve_ionised_orbitals()

std::vector< std::vector< IonisedOrbital > > Kion::solve_ionised_orbitals	(	const HF::HartreeFock *	vHF,
		const std::vector< double > &	omegas,
		double	ec_min,
		double	ec_max,
		int	Kmax,
		const std::optional< std::array< int, 2 > > &	lc_minmax,
		bool	force_rescale,
		bool	hole_particle,
		bool	force_orthog
	)

The core orbitals ionised by each of a set of energy transfers, each with the Hartree-Fock continuum states of its ejected electron.

An orbital is ionised if its ejected electron energy \( \en_c = \omega + \en_a \) lies in (ec_min, ec_max]. Its continuum states are solved at that energy, with l from continuum_l_range(), the unresolved high-r tail locally averaged (DiracODE::averageTail). In core order at each energy; parallel over (energy, orbital, l).

Parameters

vHF	Hartree-Fock potential; its core defines the orbitals.
omegas	Energy transfers, in au.
ec_min	Minimum ejected electron energy, in au.
ec_max	Maximum ejected electron energy, in au.
Kmax	Maximum multipolarity (sets the continuum l range).
lc_minmax	Optional limits on the continuum orbital l.
force_rescale	Rescale V(r) at large r for the continuum states.
hole_particle	Solve the continuum in the V^(N-1) potential of the hole.
force_orthog	Orthogonalise the continuum states to the core.

Returns: At each energy (as omegas), the ionised orbitals and their continuum states; empty if none is ionised.

◆ construct_channels()

std::vector< IonisationChannel > Kion::construct_channels	(	const std::vector< DiracSpinor > &	core,
		const std::vector< IonisedOrbital > &	ionised
	)

Every (hole orbital, ejected state) channel at one energy transfer.

Ordered by orbital (as ionised), then by continuum state. The channels point into core and ionised, which must outlive them.

Parameters

core	Core orbitals.
ionised	The ionised orbitals and their continuum states at one energy (solve_ionised_orbitals()).

Returns: The channel list.

◆ active_operators()

std::vector< std::size_t > Kion::active_operators ( const std::array< std::unique_ptr< DiracOperator::TensorOperator >, 10 > & operators )

Indices of the requested (non-null) operators of a multipole_operators() set, in its order.

◆ solve_channel_amplitudes()

std::optional< double > Kion::solve_channel_amplitudes	(	const DiracOperator::TensorOperator &	h,
		std::size_t	i_op,
		ExternalField::TDHFcntm *	rpa,
		double	omega,
		const RPAOptions &	rpa_options,
		const std::vector< IonisationChannel > &	channels,
		std::vector< ChannelAmplitudes > *	A_bare,
		std::vector< ChannelAmplitudes > *	A_rpa,
		bool	print = `false`
	)

One RPA solve: the bare and RPA amplitudes of every channel, for one operator at one (E, K, q).

h is the operator rpa was constructed with, with its rank and frequency already set. Solves the TDHF at omega (warm starting from the solver's current state), then fills amplitude i_op of each channel: bare \( \redmatel{e}{h}{a} \) in A_bare, and RPA \( \redmatel{e}{h + \delta V}{a} \) (ExternalField::TDHFcntm::dV_complex) in A_rpa. Channels for which the operator is zero by selection rules are left as they are.

If the solve did not converge (eps above RPAOptions::eps_fail, or nan), the RPA amplitude is the bare one, and the solver is cleared so that the next solve does not warm start from the failed state. Convergence is not tested for a first-order solve (RPAOptions::max_its of 1): its eps is the size of the correction, not a convergence measure. An operator with no non-zero bare amplitude in any channel (a transverse multipole at K = 0) is not solved: the RPA amplitudes are the bare ones, and nothing is returned.

Parameters

h	Multipole operator (rank and frequency set).
i_op	Its index in the multipole_operators() set.
rpa	TDHF solver for `h`.
omega	Energy transfer, in au.
rpa_options	Iterations and convergence limits (see RPAOptions).
channels	The (hole, ejected) channels (construct_channels()).
A_bare	Bare amplitudes, one per channel; entry `i_op` set.
A_rpa	RPA amplitudes, one per channel; entry `i_op` set.
print	Print the RPA iterations of the solve.

Returns: The eps of the solve (see ExternalField::TDHF::last_eps); nullopt if the RPA was not solved.

◆ solve_formFactors_RPA()

FormFactorsRPA Kion::solve_formFactors_RPA	(	const HF::HartreeFock *	vHF,
		const std::optional< std::array< int, 2 > > &	lc_minmax,
		double	ec_min,
		double	ec_max,
		bool	force_rescale,
		bool	hole_particle,
		bool	force_orthog,
		const std::vector< double > &	Egrid,
		const std::vector< double > &	qgrid,
		bool	diagonal_Eq,
		bool	low_q,
		const SphericalBessel::JL_table &	jK_tab,
		int	Kmin,
		int	Kmax,
		bool	vectorQ,
		bool	axialQ,
		bool	scalarQ,
		bool	pseudoscalarQ,
		bool	spatialQ,
		const RPAOptions &	rpa_options = `{}`
	)

Bare and RPA (core polarisation) form factors of every core orbital, with the RPA solved at every point of the grids, from the outgoing-wave TDHF.

The bare factors are exactly those of calculate_formFactors() (Hartree-Fock states only). For the RPA, each channel amplitude \( \redmatel{e}{h}{a} \) is replaced by the complex outgoing-wave amplitude \( \redmatel{e}{h + \delta V}{a} \) (see ExternalField::TDHFcntm::dV_complex), and the factors are accumulated as in accumulate_formFactors().

One RPA solve is required per (E, K, operator, q); it serves every core orbital at once, which is why the orbital loop is inside. The RPA includes every open channel; the ec limits apply to the output only.

Parallel over the coarsest axis that keeps every thread busy. Over the energies at which some orbital is ionised, when there are at least as many as threads: each thread solves its energies in full, with its own continuum states and solver, and the q points in one warm-start chain. Otherwise, within each (E, K, operator) block, over q when there are at least as many as threads (each thread owns a solver, and its solves warm start from the neighbouring q); else q runs serially with the solver's own parallelism. Prints one summary line at the start, then a progress bar.

The RPA is solved at every (E, q) given: the E and q limits of rpa_options are not applied here (see calculate_formFactors_RPA(), which restricts the grids to the region within them; what the formFactors module uses). The K limit is: multipoles above RPAOptions::K_max contribute their bare amplitudes, with no solve.

Parameters

vHF	Hartree-Fock potential; its core defines the orbitals.
lc_minmax	Optional limits on the continuum orbital l.
ec_min	Minimum ejected electron energy, in au.
ec_max	Maximum ejected electron energy, in au.
force_rescale	Rescale V(r) at large r for the continuum states.
hole_particle	Solve the continuum in the V^(N-1) potential of the hole (required here; see warning).
force_orthog	Orthogonalise the continuum states to the core.
Egrid	Energy transfer grid, in au.
qgrid	Momentum transfer grid, in au (size 1 if diagonal).
diagonal_Eq	Momentum transfer set equal to the energy transfer.
low_q	Use the low-q form of the operators.
jK_tab	Precomputed spherical Bessel table.
Kmin	Minimum multipolarity K.
Kmax	Maximum multipolarity K.
vectorQ	Calculate the vector factors.
axialQ	Calculate the axial-vector factors.
scalarQ	Calculate the scalar factor.
pseudoscalarQ	Calculate the pseudoscalar factor.
spatialQ	Calculate the spatial (E, M, L) components.
rpa_options	RPA iterations, convergence target, the eps above which a solve is discarded, and the K limit (see RPAOptions); the E and q limits are not applied.

Returns: The bare and RPA factors of each core orbital, and the worst RPA eps at each (E, q).

Note: An unconverged solve (eps above RPAOptions::eps_fail, or nan) is not used: the bare amplitude is taken for that (E, K, operator, q), and the solver is cleared so the next q does not warm start from it (see solve_channel_amplitudes()). The factors are accumulated per K for the RPA-solved K, and each (K, factor) block is corrected at its failed (E, q) points from the converged neighbours in E and q (interpolate_failed_rpa_2d()) before the sum over K; not after a first-order solve. FormFactorsRPA::eps records the worst eps at each (E, q), and n_solves, n_failed, n_interpolated the counts.; Parallel over chains, each one (E, K, operator) and a run of consecutive q (warm starts along it). The energies are taken in batches of the thread count, the continuum states of a batch held together, and the chains of a batch are one dynamically scheduled pool; the run of q is chosen to give a few chains per thread. The multipoles above the K limit are bare, one task per (E, K).

Warning: The continuum states must be those of the residual ion (hole_particle = true) for the RPA amplitude to be consistent; see ExternalField::TDHFcntm::dV_complex.

◆ calculate_formFactors_RPA()

FormFactorsRPA Kion::calculate_formFactors_RPA	(	const HF::HartreeFock *	vHF,
		const std::optional< std::array< int, 2 > > &	lc_minmax,
		double	ec_min,
		double	ec_max,
		bool	force_rescale,
		bool	hole_particle,
		bool	force_orthog,
		const std::vector< double > &	Egrid,
		const std::vector< double > &	qgrid,
		bool	diagonal_Eq,
		bool	low_q,
		const SphericalBessel::JL_table &	jK_tab,
		int	Kmin,
		int	Kmax,
		bool	vectorQ,
		bool	axialQ,
		bool	scalarQ,
		bool	pseudoscalarQ,
		bool	spatialQ,
		const RPAOptions &	rpa_options = `{}`
	)

Calculates the form factors for every core orbital, without and with RPA (core polarisation), the RPA solved only within the limits of RPAOptions.

The bare factors are those of calculate_formFactors() (Hartree-Fock states) over the full grids. The RPA is solved (solve_formFactors_RPA()) only on the region of the grids with \( E \le E_{\rm max} \) and \( q \le q_{\rm max} \) (RPAOptions::E_max, q_max), and its factors replace the bare ones there; elsewhere the RPA factors are the bare ones. Multipoles above RPAOptions::K_max are bare everywhere (applied inside solve_formFactors_RPA()). In the diagonal case (q = E/c) the q limit is a limit on E. With no E or q limit in effect the RPA is solved over the full grids directly, and the bare factors are calculated once only.

The limits exist because the RPA is expensive, and hard to converge, at high E, q, and K (many open channels; rapidly oscillating operators) where its effect on the factors is small: they confine the work to the low-E, low-q, low-K region where the core polarisation matters.

Parameters

vHF	Hartree-Fock potential; its core defines the orbitals.
lc_minmax	Optional limits on the continuum orbital l.
ec_min	Minimum ejected electron energy, in au.
ec_max	Maximum ejected electron energy, in au.
force_rescale	Rescale V(r) at large r for the continuum states.
hole_particle	Solve the continuum in the V^(N-1) potential of the hole (required here; see warning).
force_orthog	Orthogonalise the continuum states to the core.
Egrid	Energy transfer grid, in au.
qgrid	Momentum transfer grid, in au (size 1 if diagonal).
diagonal_Eq	Momentum transfer set equal to the energy transfer.
low_q	Use the low-q form of the operators.
jK_tab	Precomputed spherical Bessel table.
Kmin	Minimum multipolarity K.
Kmax	Maximum multipolarity K.
vectorQ	Calculate the vector factors.
axialQ	Calculate the axial-vector factors.
scalarQ	Calculate the scalar factor.
pseudoscalarQ	Calculate the pseudoscalar factor.
spatialQ	Calculate the spatial (E, M, L) components.
rpa_options	RPA iterations, convergence target, the eps above which a solve is discarded, and the E, q, K limits (see RPAOptions).

Returns: The bare and RPA factors of each core orbital, and the worst RPA eps at each (E, q): zero outside the RPA region.

Note: An unconverged solve within the region is treated as in solve_formFactors_RPA(): the bare amplitude is used there, and FormFactorsRPA::eps records it.

Warning: The continuum states must be those of the residual ion (hole_particle = true) for the RPA amplitude to be consistent; see ExternalField::TDHFcntm::dV_complex.

◆ count_failed_rpa()

std::pair< std::size_t, std::size_t > Kion::count_failed_rpa	(	const LinAlg::Matrix< double > &	eps,
		double	eps_fail
	)

Counts the points of an RPA eps table at which the solve failed, and those of them with a converged neighbour along the columns.

A point fails if rpa_failed(); its neighbours are the adjacent columns of the same row (the q grid for the form factors, the energy grid for photoionisation). Those with at least one converged neighbour are the points interpolate_failed_rpa() corrects.

Parameters

eps	Final RPA eps of each point (rows x columns).
eps_fail	Failure threshold; as RPAOptions::eps_fail.

Returns: {number of failed points, number of those with a converged neighbour}.

◆ interpolate_failed_rpa()

void Kion::interpolate_failed_rpa	(	const LinAlg::Matrix< double > &	eps,
		double	eps_fail,
		const LinAlg::Matrix< double > &	K_bare,
		LinAlg::Matrix< double > *	K_rpa
	)

Corrects the points at which the RPA solve failed, by interpolating the relative RPA shift from the neighbouring columns.

Where a solve fails the bare value is used, which leaves a step in an otherwise smooth function. The relative shift \( R = K_{\rm RPA}/K_{\rm bare} \) varies slowly along the columns (q for the form factors, the photon energy for photoionisation), so it is taken from the converged neighbours and applied to the bare value, \( K_{\rm RPA}(i,j) = R\,K_{\rm bare}(i,j) \):

\[ \begin{align} R &= \frac{1}{2}(R_{j-1} + R_{j+1}), \\ R &= 1 + \frac{1}{2}(R_{j\mp1} - 1), \end{align} \]

the mean of the two when both have converged; half the relative correction when only one has, since the shift is then unconstrained on the other side. A point with no converged neighbour (in a run of adjacent failures, a resonance say) is left as it is, as is one whose converged neighbours have a zero bare value. Points that fail are never used as neighbours, so the result does not depend on the order of the corrections.

Parameters

eps	Final RPA eps of each point (rows x columns).
eps_fail	Failure threshold (see rpa_failed()).
K_bare	Bare (no-RPA) values, same shape as `eps`.
K_rpa	RPA values, same shape; corrected in place.

Note: Corrects nothing with a single column, since there is no neighbour to interpolate from. Only meaningful for an iterated solve: after a first-order solve (RPAOptions::max_its of 1) eps is the size of the correction, not a convergence measure.

◆ interpolate_failed_rpa_2d()

void Kion::interpolate_failed_rpa_2d	(	const LinAlg::Matrix< double > &	eps,
		double	eps_fail,
		const LinAlg::Matrix< double > &	K_bare,
		LinAlg::Matrix< double > *	K_rpa
	)

Corrects the (E, q) points at which the RPA solve failed, by interpolating the relative RPA shift from the neighbours in E and q.

As interpolate_failed_rpa(), with the neighbours on both axes: the adjacent rows (E) and columns (q). The shift \( R = K_{\rm RPA}/K_{\rm bare} \) is the mean over the converged neighbours with a non-zero bare value (up to four); with a single one, half its relative correction. A point with none is left as it is. Failed points are never used as neighbours. With a single column (the diagonal, q = E/c, case) the neighbours are along E.

Parameters

eps	Final RPA eps of each point (E rows x q columns); a point that was not solved must not count as failed (eps of 0, or negative).
eps_fail	Failure threshold (see rpa_failed()).
K_bare	Bare (no-RPA) values, same shape as `eps`.
K_rpa	RPA values, same shape; corrected in place.

◆ factor_operators()

std::vector< std::size_t > Kion::factor_operators ( std::size_t i_factor )

The operators (indices into ChannelAmplitudes) that a form factor is built from.

The direct factors depend on one amplitude, the cross terms on two (X: t and L; Y: t5 and L5) or four (Z: E, M, E5, M5); see accumulate_formFactors(). A factor is affected by a failed RPA solve of any of its operators.

Parameters

i_factor Index into FormFactorSet.

Returns: The operator indices.

◆ check_radial_grid_q()

bool Kion::check_radial_grid_q	(	double	qmax_au,
		const Grid &	rgrid
	)

The q part of check_radial_grid only: grid spacing near r ~ a0 must resolve exp(iq.r) at qmax (au). Returns false (with advice) if not.

◆ check_radial_grid_E()

bool Kion::check_radial_grid_E	(	double	Emax_au,
		const Grid &	rgrid,
		double	alpha
	)

The E part of check_radial_grid only: grid must resolve the continuum oscillations out to rmax at Emax (au). Returns false (with advice) if not.

◆ check_radial_grid()

bool Kion::check_radial_grid	(	double	Emax,
		double	qmax,
		const Grid &	rgrid,
		double	alpha = `PhysConst::alpha`
	)

Checks the radial grid is dense enough for the continuum states and momentum transfers required.

Two independent checks, both of which print advice (larger num_points, or a different loglinear b) when they fail:

q: the grid spacing near \( r\sim a_0 \) must resolve the oscillations of \( e^{i\vb{q}\cdot\vb{r}} \) at qmax. Only a rough guide: high q may contribute negligibly, in which case error there does not matter.
E: the grid must resolve the continuum oscillations out to rmax at Emax (see DiracODE::RequiredContinuumGrid).

Parameters

Emax	Maximum continuum state energy, in au.
qmax	Maximum momentum transfer, in au.
rgrid	Radial grid to be checked (loglinear expected).
alpha	Fine-structure constant (as used by the Hartree-Fock).

Returns: False if either check fails; calculations may then be inaccurate.

◆ write_to_file_xyz()

void Kion::write_to_file_xyz	(	const std::string &	filename,
		const std::vector< double > &	E_grid,
		const std::vector< double > &	q_grid,
		const std::vector< std::string > &	titles,
		const std::vector< std::string > &	descriptions,
		std::vector< LinAlg::Matrix_view< const double > >	factors,
		Units	units = `Units::Particle`,
		int	num_digits = `6`,
		bool	diagonal = `false`
	)

Writes output file in 'xyz' form: for easy 2D interpolation.

A header of column descriptions, then one row per (E, q) point: 'E q K_1(E,q) ... K_n(E,q)'. A blank line separates each energy (which gnuplot requires, and pyplot ignores). Factors that are empty are skipped, along with their column.

Parameters

filename	Output file name.
E_grid	Energy transfer grid, in au.
q_grid	Momentum transfer grid, in au.
titles	Short column header of each factor.
descriptions	Longer description of each factor, for the header.
factors	Factors to write; must match `titles` in size.
units	Units for E and q (K is dimensionless).
num_digits	Digits printed for each value (clamped to 3 - 16).
diagonal	Momentum transfer equal to the energy transfer: q is written as alpha*E, and `q_grid` is not used.

◆ write_to_file_xyz_13()

void Kion::write_to_file_xyz_13	(	const std::string &	filename,
		const std::vector< double > &	E_grid,
		const std::vector< double > &	q_grid,
		const std::vector< std::string > &	titles,
		const std::vector< std::string > &	descriptions,
		const FormFactorSet &	K_factors,
		Units	units,
		int	num_digits,
		bool	diagonal
	)

As write_to_file_xyz, for a FormFactorSet (empty factors are skipped)

◆ write_to_file_matrix()

void Kion::write_to_file_matrix	(	const LinAlg::Matrix< double > &	K,
		const std::vector< double > &	E_grid,
		const std::vector< double > &	q_grid,
		const std::string &	filename,
		int	num_digits = `5`,
		Units	units = `Units::Particle`
	)

Writes output file in matrix form.

The E and q grids are printed first, then the entire matrix in table form: each row is a new E, each column a new q.

Note: Kept for legacy: new code should use XYZ format

Parameters

K	Factor to write, K(E,q).
E_grid	Energy transfer grid, in au.
q_grid	Momentum transfer grid, in au.
filename	Output file name.
num_digits	Digits printed for each value.
units	Units for E and q in the output (K is dimensionless).

◆ Zeff_nonrel()

double Kion::Zeff_nonrel	(	double	en,
		int	n
	)

inline

Effective charge of an orbital, from its binding energy.

\[ Z_{\rm eff} = n\sqrt{-2\en} \]

Same Zeff as used by DarkARC (see arXiv:1912.08204).

Parameters

en	Orbital energy (binding energy, negative), in au.
n	Principal quantum number.

Returns: Effective charge.

◆ rpa_failed()

bool Kion::rpa_failed	(	double	eps,
		double	eps_fail
	)

inline

A failed RPA solve: final eps above eps_fail, or nan.

◆ momentum_orbital()

MomentumOrbital Kion::momentum_orbital	(	const DiracSpinor &	Fa,
		double	en,
		double	p_min = `1.0e-3`,
		double	p_max = `1.0e4`,
		std::size_t	points_per_decade = `150`,
		double	density_cut = `1.0e-12`
	)

Fourier transforms a bound orbital onto a logarithmic momentum grid (MomentumOrbital).

The grid is logarithmic, points_per_decade points per factor of 10 in p, running upwards from p_min, and stops (at most at p_max) once the density has fallen below density_cut times its peak, beyond the peak. The radial integrals run over the extent of the orbital (its max_pt()), with the cell-averaged Bessel functions of SphericalBessel::fillBesselVec_kr, so the high-p tail (from small r) is integrated correctly on a coarse grid.

Parameters

Fa	Bound orbital (as used in the matrix elements); its occupation sets the electron count.
en	Binding energy to record (au): that of the real orbital, which may differ from Fa.en() for a Zeff model state.
p_min	First momentum (au).
p_max	Largest momentum considered (au).
points_per_decade	Grid points per factor of 10 in p.
density_cut	Stop once density < density_cut times the peak.

Returns: The momentum-space orbital.

◆ planewave_traces()

std::array< double, 11 > Kion::planewave_traces	(	double	p,
		double	F,
		double	G,
		int	kappa,
		double	pf,
		double	q,
		double	ef,
		double	alpha
	)

Integrand of the all-K plane-wave response at bound-electron momentum p: the Dirac traces of every factor.

With z along q, \( \vb{p}_f = \vb{p} + \vb{q} \), and the angle between p and q fixed by energy conservation,

\[ \cos\theta_{pq} = \frac{p_f^2 - p^2 - q^2}{2pq}, \]

returns the traces over the Dirac indices

\[ T = {\rm Tr}\big[(E_f + c\,\vb{\alpha}\cdot\vb{p}_f + \beta mc^2)\, \Gamma\,\rho\,\Gamma'^\dagger\big] \]

for the momentum-space density matrix of the shell without its \( N_a/(8\pi) \) prefactor (2x2 blocks; see planewave_formFactors()),

\[ \rho = \begin{pmatrix} \tilde f_a^2 & s_\kappa\,\tilde f_a\tilde g_a\,\vb{\sigma}\cdot\hat{\vb{p}} \\ s_\kappa\,\tilde f_a\tilde g_a\,\vb{\sigma}\cdot\hat{\vb{p}} & \tilde g_a^2 \end{pmatrix}, \]

where \( s_\kappa = \kappa/|\kappa| \). The traces are returned in the order {V: 00, 33, 11+22, 03; A: 00, 33, 11+22, 03; Im VA 12; S; P}. See planewave_formFactors() for the vertices and the use.

Parameters

p	Bound-electron momentum (au).
F	Transformed large component, f~_a, at p.
G	Transformed small component, g~_a, at p.
kappa	Dirac quantum number of the orbital.
pf	Ejected-electron momentum, p_f (au).
q	Momentum transfer (au).
ef	Ejected-electron kinetic energy (au).
alpha	Fine-structure constant.

Returns: The 11 traces.

◆ planewave_formFactors()

std::array< double, 13 > Kion::planewave_formFactors	(	const MomentumOrbital &	orb,
		double	E,
		double	q,
		double	alpha
	)

The 13 form factors of one orbital with a free (plane-wave) ejected electron, summed over all multipoles, at energy transfer E and momentum transfer q.

The ejected electron is free (the potential does not act during the collision), with kinetic energy, momentum, and total energy

\[ \begin{align} \en_f &= E + \en_a, \\ p_f &= \sqrt{\en_f(2 + \alpha^2\en_f)}, \\ E_f &= \en_f + mc^2, \end{align} \]

and the bound electron has the momentum distribution of the orbital (MomentumOrbital). The factors of the orbital are the Cartesian components (z along q) of its response tensor,

\[ R^{\mu\nu}_a = x_a \sum_{m_a}\sum_f J^\mu_{fa}\, J^{\nu *}_{fa}, \]

summed over the final plane waves at energy \( \en_f \) and over the closed shell. The sum reduces to a single integral over the bound-electron momentum p, the angle between p and q being fixed by energy conservation (planewave_traces()):

\[ R^{\mu\nu}_a = \frac{1}{8\pi^2 c^2 q}\int_{|p_f-q|}^{p_f+q} p\,dp\; {\rm Tr}\big[(E_f + c\,\vb{\alpha}\cdot\vb{p}_f + \beta mc^2)\, \Gamma^\mu \rho_a(p)\, \Gamma^{\nu\dagger}\big], \]

with \( \vb{p}_f = \vb{p} + \vb{q} \), the trace over the Dirac indices, and

\[ \rho_a = \frac{N_a}{8\pi}\begin{pmatrix} \tilde f_a^2 & s_\kappa\,\tilde f_a\tilde g_a\,\vb{\sigma}\cdot\hat{\vb{p}} \\ s_\kappa\,\tilde f_a\tilde g_a\,\vb{\sigma}\cdot\hat{\vb{p}} & \tilde g_a^2 \end{pmatrix} \]

the momentum-space density matrix of the shell ( \( N_a \) electrons). The vertex \( \Gamma^\mu = \gamma^0\gamma^\mu\tilde\gamma \) is that of the operator:

vector: \( (1, \vb{\alpha}) \)
axial: \( (\gamma^5, \vb{\Sigma}) \)
scalar: \( \beta \)
pseudoscalar: \( i\beta\gamma^5 \)

The factors are the components

\[ \begin{align} V_T &= R^{00}, \\ V_L &= R^{33}, \\ V_E + V_M &= R^{11} + R^{22}, \\ X &= -{\rm Re}\,R^{03}, \\ A_T &= R^{00}_A, \\ A_L &= R^{33}_A, \\ A_E + A_M &= R^{11}_A + R^{22}_A, \\ Y &= {\rm Re}\,R^{03}_A, \\ Z &= 2\,{\rm Im}\,R^{12}_{VA}, \\ S &= R_{SS}, \\ P &= R_{PP}, \end{align} \]

where Z is the antisymmetric transverse vector-axial interference. With the couplings in the vertex, \( \gamma^0\gamma^\mu(c_V - c_A\gamma^5) \), the tensor is \( c_V^2 R_V + c_A^2 R_A - c_V c_A (R_{VA} + R_{AV}) \), and its components are the coefficients of the spin-summed cross section:

\[ \begin{align} R^{00} &= c_V^2 V_T + c_A^2 A_T, \\ R^{33} &= c_V^2 V_L + c_A^2 A_L, \\ R^{11} + R^{22} &= c_V^2 (V_E + V_M) + c_A^2 (A_E + A_M), \\ -{\rm Re}\,R^{03} &= c_V^2 X - c_A^2 Y, \\ -{\rm Im}\,R^{12} &= c_V c_A Z, \end{align} \]

(for an isotropic target \( R^{21}_{VA} = -R^{12}_{VA} \), so the VA and AV terms contribute \( {\rm Im}\,R^{12}_{VA} \) each), and \( R = c_S^2 S + c_P^2 P \) for the scalar-pseudoscalar vertex \( c_S\gamma^0 + i c_P\gamma^0\gamma^5 \). For an electron at rest the factors reduce to the free-electron values, e.g.

\[ \begin{align} V_L &= (E/qc)^2\, V_T, \\ X &= -(E/qc)\, V_T, \\ A_L &= V_T. \end{align} \]

The electric and magnetic multipoles are not separately defined by the Cartesian tensor (only their sum enters a spin-summed cross section): the whole transverse response is returned as the electric factor, and the magnetic factors are zero.

Parameters

orb	Momentum-space orbital (momentum_orbital()).
E	Energy transfer (au).
q	Momentum transfer (au).
alpha	Fine-structure constant.

Returns: The factors, in FormFactorSet order; all zero if the orbital is not ionised ( \( \en_f \le 0 \)).

◆ captured_fraction()

double Kion::captured_fraction	(	const DiracSpinor &	Fa,
		int	Kmax,
		std::size_t	iq,
		const SphericalBessel::JL_table &	jK_tab
	)

Fraction of the norm of \( e^{i\vb{q}\cdot\vb{r}}\psi_a \) carried by the multipoles K <= Kmax.

\[ c_a(q) = \sum_{K=0}^{K_{\rm max}} (2K+1)\int (f_a^2 + g_a^2)\, j_K(qr)^2\, dr, \]

which tends to 1 by the identity

\[ \sum_{K=0}^{\infty} (2K+1)\, j_K(x)^2 = 1. \]

The remainder 1 - c is the norm of the part of the state that the multipoles above Kmax carry, and bounds what their omission can cost: the plane-wave completion is skipped where it is negligible.

Parameters

Fa	Bound orbital.
Kmax	Maximum multipolarity of the sum.
iq	Momentum-transfer index in the Bessel table.
jK_tab	Spherical Bessel table, on the orbital's grid.

Returns: The captured fraction.

◆ calculate_ridge_correction()

std::vector< FormFactorSet > Kion::calculate_ridge_correction	(	const HF::HartreeFock *	vHF,
		const std::vector< DiracSpinor > &	bound_states,
		double	ec_min,
		double	ec_max,
		const std::vector< double > &	Egrid,
		const std::vector< double > &	qgrid,
		const SphericalBessel::JL_table &	jK_tab,
		int	Kmax,
		bool	vectorQ,
		bool	axialQ,
		bool	scalarQ,
		bool	pseudoscalarQ,
		bool	spatialQ,
		double	ridge_eps
	)

The Bethe-ridge correction: completes the multipole sum of every core orbital to all K with plane waves.

The multipole sum truncated at Kmax misses the response near the quasi-free (Bethe) ridge,

\[ E \approx \frac{q^2}{2m}, \]

where multipoles up to \( K \sim q\,r_a \) contribute. There the ejected electron is fast and its high partial waves are undistorted, so the missing part is supplied by plane waves. Per orbital and factor, the returned correction is

\[ \Delta K_a = K_a^{\rm PW}({\rm all}\ K) - K_a^{\rm PW}(K \le K_{\rm max}), \]

the all-K plane-wave response (planewave_formFactors()) minus the same multipole sum evaluated with free Dirac spherical waves (DiracODE::freeDirac()), i.e. exactly the plane-wave multipoles above Kmax. Added to the distorted-wave factors of calculate_formFactors(), the total is

\[ K = K^{\rm DW}(K \le K_{\rm max}) + \Delta K. \]

Off the ridge the correction vanishes by itself; on the ridge the result is independent of Kmax once the low multipoles, which carry the distortion, are included.

The plane-wave sum always runs K = 0..Kmax over the full continuum l range, with the same operators, Bessel table, bound states, and continuum truncation (max_pt) as the distorted-wave sum, so that the K <= Kmax parts cancel to grid accuracy. It has no hole-particle or orthogonalisation (distorted-wave physics only). Bound-bound (Pauli) leakage of the plane waves is confined to K <= j_a + j_b and cancels as long as Kmax >= 2 j_max of the core (checked). The squared factors are clamped at zero (they are sums of squares; a negative value is noise from the subtraction). The whole transverse correction goes into the electric factors (see planewave_formFactors()).

The correction is computed only where it can matter: for each orbital, at the momentum transfers where the multipoles above Kmax carry more than ridge_eps of the norm (captured_fraction()), at the energies where the orbital is ionised (ejected energy within ( ec_min, ec_max ), and, at each energy, only where |p_f - q| lies within the momentum grid of the orbital (elsewhere, far off the ridge, the plane-wave response is zero by construction). Parallel over the (orbital, E) pairs. Prints, per orbital, the momentum transfer above which the correction is active.

Parameters

vHF	Hartree-Fock; its core defines the orbitals.
bound_states	Bound state of each orbital as used in the matrix elements (model_bound_states()); must be the same as used for the distorted-wave factors.
ec_min	Minimum ejected electron energy, in au.
ec_max	Maximum ejected electron energy, in au.
Egrid	Energy transfer grid, in au.
qgrid	Momentum transfer grid, in au.
jK_tab	Spherical Bessel table on `qgrid` (as used for the distorted-wave factors).
Kmax	Maximum multipolarity of the distorted-wave sum.
vectorQ	Vector factors.
axialQ	Axial-vector factors.
scalarQ	Scalar factor.
pseudoscalarQ	Pseudoscalar factor.
spatialQ	Spatial (E, M, L) components.
ridge_eps	Skip (orbital, q) where the multipoles above Kmax carry less than this fraction of the norm.

Returns: The correction of each core orbital, indexed as the core, with the allocation pattern of allocate_formFactors().

Warning: Requires Kmax >= 2 j_max(core): if not, prints a warning and returns a zero correction. Not defined for the diagonal (q = E/c) case, which never approaches the ridge.

Detailed Description

Classes

Typedefs

Enumerations

Functions

Class Documentation

◆ Kion::FormFactorsRPA

◆ Kion::IonisationChannel

◆ Kion::IonisedOrbital

◆ Kion::MomentumOrbital

◆ Kion::RPAOptions

Typedef Documentation

◆ FormFactorSet

◆ ChannelAmplitudes

Enumeration Type Documentation

◆ Coupling

◆ OutputFormat

◆ Units

◆ AtomicMethod

Function Documentation

◆ parseStatesMethod() [1/2]

◆ parseStatesMethod() [2/2]

◆ calculateK_nk()

◆ allocate_formFactors()

◆ multipole_operators()

◆ accumulate_formFactors()

◆ model_bound_states()

◆ accumulate_multipole_sum()

◆ add_formFactors()

◆ calculate_formFactors()

◆ continuum_l_range()

◆ solve_ionised_orbitals()

◆ construct_channels()

◆ active_operators()

◆ solve_channel_amplitudes()

◆ solve_formFactors_RPA()

◆ calculate_formFactors_RPA()

◆ count_failed_rpa()

◆ interpolate_failed_rpa()

◆ interpolate_failed_rpa_2d()

◆ factor_operators()

◆ check_radial_grid_q()

◆ check_radial_grid_E()

◆ check_radial_grid()

◆ write_to_file_xyz()

◆ write_to_file_xyz_13()

◆ write_to_file_matrix()

◆ Zeff_nonrel()

◆ rpa_failed()

◆ momentum_orbital()

◆ planewave_traces()

◆ planewave_formFactors()

◆ captured_fraction()

◆ calculate_ridge_correction()