Detailed Description

Core-polarisation (RPA) corrections to matrix elements of an external field.

External field: Mixed-states + Core Polarisation.

Provides classes and functions for computing all-order core-polarisation corrections to matrix elements of an external field operator.

In the presence of an external field of frequency \( \omega \), the time-dependent operator is

\[ T(t) = t_+ e^{-i\omega t} + t_- e^{+i\omega t}, \]

where \( t_+ = t^k_q(\omega) \) is an irreducible tensor operator of rank \( k \) and \( t_- = t_+^\dag(-\omega) \). Each orbital, including those in the core, acquires a first-order perturbation,

\[ \delta\phi_a(t) = \varphi^a_+ e^{-i\omega t} + \varphi^a_- e^{+i\omega t}. \]

Since the core orbitals are perturbed, the Hartree-Fock potential is also perturbed:

\[ \delta V_\pm \phi_i = \sum_a^{\rm core} \left[ \matel{\phi_a}{Q}{\varphi^a_+}\phi_i - \matel{\phi_a}{Q}{\phi_i}\varphi^a_+ + \matel{\varphi^a_-}{Q}{\phi_a}\phi_i - \matel{\varphi^a_-}{Q}{\phi_i}\phi_a \right]. \]

The resulting core-polarisation corrections to matrix elements are given by

\[ \matel{b}{t_\pm}{a} \to \matel{b}{t_\pm + \delta V_\pm}{a}, \]

where \( \delta V_\pm \) is the correction to the HF potential arising from the perturbed core orbitals \( \{\varphi^a_\pm\} \).

Since \( \delta V \) is solved self-consistently, this accounts for core polarisation to all orders in the Coulomb interaction.

Two equivalent methods are implemented:

TDHF (TDHF, TDHFbasis): solves the TDHF equations
\[ (h_{\rm HF} - \en_a \mp \omega)\varphi^a_\pm = -(t_\pm + \delta V_\pm - \delta\en^a_\pm)\phi_a \]
self-consistently for all core orbitals. The corrections can also be found via the basis expansion
\[ \varphi^a_\pm = \sum_n \frac{\ket{n}\matel{n}{t_\pm + \delta V_\pm}{a}} {\en_a - \en_n \pm \omega}. \]
See Dzuba et al. (1984).
Diagram/Goldstone (DiagramRPA): evaluates the four core-polarisation diagrams directly,
\[ \redmatel{w}{\delta V_\pm}{v} = \sum_{na} \frac{(-1)^{k+w-n}}{[k]} \left( \frac{T_{na}\,W^k_{wavn}}{\en_a - \en_n \pm \omega} + \frac{(-1)^{a-n} T_{an}\,W^k_{wnva}}{\en_a - \en_n \mp \omega} \right), \]
with \( T_{ij} = t_{ij} + \redmatel{i}{\delta V_\pm}{j} \) iterated to convergence. See Johnson et al., Phys. Rev. A 21, 409 (1980).

Use make_rpa to construct the appropriate object from a method string, and calcMatrixElements to compute matrix elements including the correction.

Classes
class	CorePolarisation
	Virtual base class for core-polarisation (RPA); computes dV corrections. More...

class	DiagramRPA
	RPA correction to matrix elements using the diagram technique. More...

struct	MEdata
	Small struct to store a single calculated reduced matrix element with RPA correction. More...

class	TDHF
	Uses TDHF to include core-polarisation (RPA) corrections to matrix elements of an external field operator. More...

class	TDHFbasis
	Like TDHF, but solves the TDHF equations via basis expansion. More...

Enumerations
enum class	Method { TDHF , basis , diagram , none , Error }
	Available RPA/core-polarisation methods. More...

enum class	dPsiType { X , Y }
	Selects the perturbed orbital: X = varphi_+, Y = varphi_-. More...

enum class	StateType { bra , ket }
	Whether the state is a bra or ket. More...

Functions
std::vector< MEdata >	calcMatrixElements (const std::vector< DiracSpinor > &b_orbs, const std::vector< DiracSpinor > &a_orbs, DiracOperator::TensorOperator const h, CorePolarisation const dV=nullptr, double omega=0.0, bool each_freq=false, bool diagonal=true, bool off_diagonal=true, bool calculate_both=false)
	Calculates reduced matrix elements <a\|\|h\|\|b> for lists of orbitals.

Coulomb::meTable< double >	me_table (const std::vector< DiracSpinor > &a_orbs, const std::vector< DiracSpinor > &b_orbs, const DiracOperator::TensorOperator h, const CorePolarisation dV=nullptr, const MBPT::StructureRad *srn=nullptr, std::optional< double > omega={})
	Builds a lookup table of reduced matrix elements <a\|\|h\|\|b>.

std::unique_ptr< CorePolarisation >	make_rpa (const std::string &method, const DiracOperator::TensorOperator h, const HF::HartreeFock vhf, bool print=false, const std::vector< DiracSpinor > &basis={}, const std::string &identity="")
	Factory function to construct a core-polarisation (RPA) object.

std::vector< MEdata >	calcMatrixElements (const std::vector< DiracSpinor > &orbs, DiracOperator::TensorOperator const h, CorePolarisation const dV=nullptr, double omega=0.0, bool each_freq=false, bool diagonal=true, bool off_diagonal=true, bool calculate_both=false)
	Calculates reduced matrix elements for a single list of orbitals.

Coulomb::meTable< double >	me_table (const std::vector< DiracSpinor > &a_orbs, const DiracOperator::TensorOperator h, const CorePolarisation dV=nullptr, const MBPT::StructureRad *srn=nullptr, std::optional< double > omega={})
	Builds a matrix element table for a single set of orbitals.

Method	ParseMethod (std::string_view str)
	Parses method string to Method enum (case-insensitive)

DiracSpinor	solveMixedState (const DiracSpinor &Fa, double omega, const std::vector< double > &vl, double alpha, const std::vector< DiracSpinor > &core, const DiracSpinor &Fs, double eps_target=1.0e-9, const MBPT::CorrelationPotential const Sigma=nullptr, const HF::Breit const VBr=nullptr, const std::vector< double > &H_mag={})
	Solves the inhomogeneous TDHF (mixed-states) equation for perturbed orbital dF.

void	solveMixedState (DiracSpinor &dF, const DiracSpinor &Fa, double omega, const std::vector< double > &vl, double alpha, const std::vector< DiracSpinor > &core, const DiracSpinor &Fs, double eps_target=1.0e-9, const MBPT::CorrelationPotential const Sigma=nullptr, const HF::Breit const VBr=nullptr, const std::vector< double > &H_mag={})
	As solveMixedState(), but updates an existing solution `dF` in place.

DiracSpinor	solveMixedState (const DiracSpinor &Fa, double omega, const DiracSpinor &Fs, const HF::HartreeFock const hf, double eps_target=1.0e-9, const MBPT::CorrelationPotential const Sigma=nullptr)
	Solves Mixed States (TDHF) equation. Overload; takes hf object.

void	solveMixedState (DiracSpinor &dF, const DiracSpinor &Fa, double omega, const DiracSpinor &Fs, const HF::HartreeFock const hf, double eps_target=1.0e-9, const MBPT::CorrelationPotential const Sigma=nullptr)
	Solves Mixed States (TDHF) equation. Overload; takes hf object.

DiracSpinor	solveMixedState_basis (const DiracSpinor &Fa, const DiracSpinor &hFa, double omega, const std::vector< DiracSpinor > &basis)
	Solves for dF via explicit sum over basis; mainly for tests.

Variables
constexpr bool	print_final_eps = false

constexpr bool	print_each_eps = false

Enumeration Type Documentation

◆ Method

enum class ExternalField::Method

strong

Available RPA/core-polarisation methods.

◆ dPsiType

enum class ExternalField::dPsiType

strong

Selects the perturbed orbital: X = varphi_+, Y = varphi_-.

Corresponds to the two first-order corrections to a core orbital,

\[ \delta\phi_a(t) = \varphi^a_+ e^{-i\omega t} + \varphi^a_- e^{+i\omega t}. \]

X selects \( \varphi^a_+ \) (absorption/forward), Y selects \( \varphi^a_- \) (emission/backward).

◆ StateType

enum class ExternalField::StateType

strong

Whether the state is a bra or ket.

Function Documentation

◆ calcMatrixElements() [1/2]

std::vector< MEdata > ExternalField::calcMatrixElements	(	const std::vector< DiracSpinor > &	b_orbs,
		const std::vector< DiracSpinor > &	a_orbs,
		DiracOperator::TensorOperator *const	h,
		CorePolarisation *const	dV = `nullptr`,
		double	omega = `0.0`,
		bool	each_freq = `false`,
		bool	diagonal = `true`,
		bool	off_diagonal = `true`,
		bool	calculate_both = `false`
	)

Calculates reduced matrix elements <a||h||b> for lists of orbitals.

Computes reduced matrix elements \( \redmatel{a}{h}{b} \) for all non-zero combinations of states from a_orbs and b_orbs. Optionally includes RPA corrections via dV.

The transition frequency is set to \( \omega_{ab} = \varepsilon_a - \varepsilon_b \) for off-diagonal elements, and 0 for diagonal elements.

If each_freq is true, the operator and dV are updated at each transition frequency individually. Otherwise, they are updated once using omega.

Selection rules are applied via the operator's isZero() method. For odd-parity operators, only elements with the even-parity state on the right are included (unless calculate_both is true). For even-parity operators, each off-diagonal pair is included once (upper triangle, unless calculate_both is true). When b_orbs != a_orbs, calculate_both is set to true automatically.

Parameters

b_orbs	Bra states (final states, index b).
a_orbs	Ket states (initial states, index a).
h	Pointer to the tensor operator.
dV	Optional RPA/core-polarisation correction. If nullptr, no correction is applied.
omega	Frequency used to update frequency-dependent operators and dV when `each_freq` is false.
each_freq	If true, update the operator and dV at the natural transition frequency for each element.
diagonal	If true, include diagonal elements \( \redmatel{a}{h}{a} \).
off_diagonal	If true, include off-diagonal elements.
calculate_both	If true, include both \( \redmatel{a}{h}{b} \) and \( \redmatel{b}{h}{a} \) for off-diagonal pairs. Forced true when `b_orbs` != `a_orbs`.

Returns: Vector of MEdata structs, one per non-zero matrix element.

◆ me_table() [1/2]

Coulomb::meTable< double > ExternalField::me_table	(	const std::vector< DiracSpinor > &	a_orbs,
		const std::vector< DiracSpinor > &	b_orbs,
		const DiracOperator::TensorOperator *	h,
		const CorePolarisation *	dV = `nullptr`,
		const MBPT::StructureRad *	srn = `nullptr`,
		std::optional< double >	omega = `{}`
	)

Builds a lookup table of reduced matrix elements <a||h||b>.

Fills and returns a Coulomb::meTable<double> with reduced matrix elements

\[ t_{ab} = \redmatel{a}{h}{b} + \delta V_{ab} + \delta_{\rm SR}^{ab} \]

for all non-zero pairs from a_orbs and b_orbs.

Optionally includes RPA correction dV and structure radiation srn. The symmetry-conjugate \( \redmatel{b}{h}{a} \) is also stored via symm_sign(). Table is filled with OpenMP parallelisation.

Parameters

a_orbs	Bra states.
b_orbs	Ket states.
h	Pointer to the (const) tensor operator.
dV	Optional RPA correction. If nullptr, not applied.
srn	Optional structure radiation/normalisation correction. If nullptr, not applied.
omega	Only used for Structure Radiation. If set, use this fixed frequency for all elements; otherwise uses \|eb - ea\|.

Returns: meTable containing t_ab for all non-zero pairs (and conjugates).

Note: For frequency-dependent operators, omega should be set before calling this function.

◆ make_rpa()

std::unique_ptr< CorePolarisation > ExternalField::make_rpa	(	const std::string &	method,
		const DiracOperator::TensorOperator *	h,
		const HF::HartreeFock *	vhf,
		bool	print = `false`,
		const std::vector< DiracSpinor > &	basis = `{}`,
		const std::string &	identity = `""`
	)

Factory function to construct a core-polarisation (RPA) object.

Parses method and returns a std::unique_ptr<CorePolarisation> of the appropriate type. Returns nullptr if method is "none" or "false".

Supported methods (case-insensitive):

"TDHF": time-dependent Hartree-Fock (TDHF).
"basis": TDHF solved in a basis set.
"diagram": diagram RPA.
"none", "false", "": no RPA; returns nullptr.

If the method string is not recognised, prints an error and defaults to none.

Parameters

method	String specifying the RPA method (see above).
h	Pointer to the operator for which RPA is to be solved.
vhf	Pointer to the Hartree-Fock object (provides core potential).
print	If true, print a brief description of the chosen method.
basis	Basis set for basis/diagram methods (ignored for TDHF).
identity	Identifier string passed to DiagramRPA (e.g. for caching).

Returns: Unique pointer to the constructed CorePolarisation object, or nullptr if RPA is disabled.

Warning: An unrecognised method string triggers an error message and falls through to no RPA rather than throwing.

◆ calcMatrixElements() [2/2]

std::vector< MEdata > ExternalField::calcMatrixElements	(	const std::vector< DiracSpinor > &	orbs,
		DiracOperator::TensorOperator *const	h,
		CorePolarisation *const	dV = `nullptr`,
		double	omega = `0.0`,
		bool	each_freq = `false`,
		bool	diagonal = `true`,
		bool	off_diagonal = `true`,
		bool	calculate_both = `false`
	)

inline

Calculates reduced matrix elements for a single list of orbitals.

Convenience overload; calls calcMatrixElements(orbs, orbs, ...) with both bra and ket taken from the same set orbs.

◆ me_table() [2/2]

Coulomb::meTable< double > ExternalField::me_table	(	const std::vector< DiracSpinor > &	a_orbs,
		const DiracOperator::TensorOperator *	h,
		const CorePolarisation *	dV = `nullptr`,
		const MBPT::StructureRad *	srn = `nullptr`,
		std::optional< double >	omega = `{}`
	)

inline

Builds a matrix element table for a single set of orbitals.

Convenience overload; calls me_table(a_orbs, a_orbs, ...) with both bra and ket taken from a_orbs.

◆ ParseMethod()

Method ExternalField::ParseMethod ( std::string_view str )

inline

Parses method string to Method enum (case-insensitive)

◆ solveMixedState() [1/4]

DiracSpinor ExternalField::solveMixedState	(	const DiracSpinor &	Fa,
		double	omega,
		const std::vector< double > &	vl,
		double	alpha,
		const std::vector< DiracSpinor > &	core,
		const DiracSpinor &	Fs,
		double	eps_target = `1.0e-9`,
		const MBPT::CorrelationPotential *const	Sigma = `nullptr`,
		const HF::Breit *const	VBr = `nullptr`,
		const std::vector< double > &	H_mag = `{}`
	)

Solves the inhomogeneous TDHF (mixed-states) equation for perturbed orbital dF.

Solves

\[ (h_{\rm HF} - \en_a \mp \omega)\delta F + F_S = 0 \]

for \( \delta F \), where \( F_S \) is the source term. Typically

\[ F_S = (t_\pm + \delta V_\pm - \delta\en^a_\pm)\phi_a. \]

The angular momentum \( \kappa \) of the solution is that of Fs.
\( t \): Extenral field operator
\( \delta V_\pm \): core polarisation correction (see CorePolarisation)
Solved iteratively using the Green's function method.

Parameters

Fa	Unperturbed orbital \( \phi_a \).
omega	External-field frequency \( \omega \).
vl	Local potential (nuclear + direct).
alpha	Fine-structure constant.
core	Core electrons (for exchange).
Fs	Source term \( F_S \) (note sign: this is \( h\phi_a \), not \( -h\phi_a \)).
eps_target	Convergence goal for the inhomogeneous ODE solver.
Sigma	Optional correlation potential.
VBr	Optional Breit interaction.
H_mag	Magnetic part of QED radiative potential (electric part should be included in `vl`).

Returns: Perturbed orbital \( \delta F \).

◆ solveMixedState() [2/4]

void ExternalField::solveMixedState	(	DiracSpinor &	dF,
		const DiracSpinor &	Fa,
		double	omega,
		const std::vector< double > &	vl,
		double	alpha,
		const std::vector< DiracSpinor > &	core,
		const DiracSpinor &	Fs,
		double	eps_target = `1.0e-9`,
		const MBPT::CorrelationPotential *const	Sigma = `nullptr`,
		const HF::Breit *const	VBr = `nullptr`,
		const std::vector< double > &	H_mag = `{}`
	)

As solveMixedState(), but updates an existing solution dF in place.

Starts from dF as an initial guess rather than zero; converges faster if dF is already an approximate solution (e.g., from a nearby frequency).

◆ solveMixedState() [3/4]

DiracSpinor ExternalField::solveMixedState	(	const DiracSpinor &	Fa,
		double	omega,
		const DiracSpinor &	hFa,
		const HF::HartreeFock *const	hf,
		double	eps_target,
		const MBPT::CorrelationPotential *const	Sigma
	)

Solves Mixed States (TDHF) equation. Overload; takes hf object.

◆ solveMixedState() [4/4]

void ExternalField::solveMixedState	(	DiracSpinor &	dF,
		const DiracSpinor &	Fa,
		double	omega,
		const DiracSpinor &	hFa,
		const HF::HartreeFock *const	hf,
		double	eps_target,
		const MBPT::CorrelationPotential *const	Sigma
	)

Solves Mixed States (TDHF) equation. Overload; takes hf object.

◆ solveMixedState_basis()

DiracSpinor ExternalField::solveMixedState_basis	(	const DiracSpinor &	Fa,
		const DiracSpinor &	hFa,
		double	omega,
		const std::vector< DiracSpinor > &	basis
	)

Solves for dF via explicit sum over basis; mainly for tests.

\[ \delta F = \sum_n \frac{\ket{n}\matel{n}{F_S}{a}}{\en_a - \en_n \pm \omega} \]

where hFa is the already-evaluated source spinor \( F_S \).

Detailed Description

Classes

Enumerations

Functions

Variables

Enumeration Type Documentation

◆ Method

◆ dPsiType

◆ StateType

Function Documentation

◆ calcMatrixElements() [1/2]

◆ me_table() [1/2]

◆ make_rpa()

◆ calcMatrixElements() [2/2]

◆ me_table() [2/2]

◆ ParseMethod()

◆ solveMixedState() [1/4]

◆ solveMixedState() [2/4]

◆ solveMixedState() [3/4]

◆ solveMixedState() [4/4]

◆ solveMixedState_basis()