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 Amplitudes::matrix_elements 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...

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...

class	TDHFcntm
	TDHF (RPA) core polarisation above ionisation threshold(s): the open channels are solved with the outgoing-wave boundary condition, so the corrections and the corrected matrix elements are complex. 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::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="", const DiracOperator::TensorOperator *h_minus=nullptr)
	Factory function to construct a core-polarisation (RPA) object.

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

std::vector< const DiracSpinor * >	conditioning_states (const std::vector< DiracSpinor > &core, const DiracSpinor &Fa, int kappa, double e0)
	Find bound states of the solve channel that make (h_l - e0) near-singular.

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.

std::vector< double >	anderson_coefficients (const std::vector< std::vector< double > > &gram)
	Anderson mixing coefficients for a fixed-point iteration x -> G(x).

std::vector< double >	anderson_coefficients (const std::vector< DiracSpinor > &residuals)
	As above, from the residual spinors directly (Gram matrix of their inner products).

void	solveMixedState_cntm (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 HF::Breit *const VBr=nullptr, const std::vector< double > &H_mag={})
	Bound mixed-state solve with Anderson acceleration; the bound channels of TDHFcntm.

void	solveContinuumMixedState (DiracSpinor phi, DiracSpinor Freg, DiracSpinor Firr, double K, 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 DiracSpinor *const Fhole=nullptr)
	Continuum (en_+ > 0) mixed-state solve with the standing-wave boundary condition and the non-local exchange iterated in the source: the open channels of TDHFcntm.

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

◆ 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 = `""`,
		const DiracOperator::TensorOperator *	h_minus = `nullptr`
	)

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 forward operator ( \( t_+ \)).
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).
h_minus	Pointer to the backward operator ( \( t_- \)); if nullptr (default), `h` is used for both. See TDHF constructor for when this is needed.

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.

◆ ParseMethod()

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

inline

Parses method string to Method enum (case-insensitive)

◆ conditioning_states()

std::vector< const DiracSpinor * > ExternalField::conditioning_states	(	const std::vector< DiracSpinor > &	core,
		const DiracSpinor &	Fa,
		int	kappa,
		double	e0
	)

Find bound states of the solve channel that make (h_l - e0) near-singular.

For the channel of kappa kappa and energy \( e_0 = \en_a \pm \omega \), the radial operator \( (h_l - e_0) \) is (near-)singular for components along any bound core state \( \phi_m \) of the same kappa with \( \en_m \approx e_0 \): exactly singular for the diagonal ( \( \phi_a \), \( \omega = 0 \)) case, near-singular for e.g. fine-structure partners. Those components cannot be resolved reliably by the Green's-function solve, so they are projected out of the source (forces the solution orthogonal to them), they should be restore the off-diagonal ones analytically afterwards.

The set is the same-kappa core states satisfying a relative nearness criterion \( |e_0 - \en_m| < \eta\,|e_0 + \en_m| \) ( \( \eta = 0.2 \)), plus Fa itself when it shares the channel kappa (the \( \matel{a}{\delta F}{} = 0 \) / left-orthogonality constraint).

◆ 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).

Note: Near-resonant channels are handled automatically: \( (h_{\rm HF} - \en_a \mp \omega) \) is (near-)singular for components along any same-kappa bound state with \( \en_m \approx \en_a \pm \omega \) (the diagonal Fa, fine-structure partners, etc.). These are projected out of the source, the solution is forced orthogonal to them, and the off-diagonal components are restored analytically. The caller need not pre-condition Fs.

◆ anderson_coefficients() [1/2]

std::vector< double > ExternalField::anderson_coefficients ( const std::vector< std::vector< double > > & gram )

Anderson mixing coefficients for a fixed-point iteration x -> G(x).

Given the residuals \( r_k = G(x_k) - x_k \) of the stored iterates (oldest first) through their Gram matrix \( B_{kl} = \braket{r_k}{r_l} \), returns the \( c_k \) minimising

\[ \Big\|\sum_k c_k r_k\Big\|^2 \quad\text{subject to}\quad \sum_k c_k = 1, \]

via the bordered system [B 1; 1^T 0][c; lambda] = [0; 1]. The next iterate is \( x = \sum_k c_k G(x_k) \). For a linear map this is equivalent to (truncated) GMRES: it converges wherever the linear problem is non-singular, where damped iteration may not.

A saturated history makes B ill-conditioned: the oldest entries are then dropped, and the returned vector holds one coefficient per KEPT entry, the last c.size() entries of the history (the caller drops the same entries from its own history). Empty if no entry is usable (take a plain step).

◆ anderson_coefficients() [2/2]

std::vector< double > ExternalField::anderson_coefficients ( const std::vector< DiracSpinor > & residuals )

As above, from the residual spinors directly (Gram matrix of their inner products).

◆ solveMixedState_cntm()

void ExternalField::solveMixedState_cntm	(	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 HF::Breit *const	VBr = `nullptr`,
		const std::vector< double > &	H_mag = `{}`
	)

Bound mixed-state solve with Anderson acceleration; the bound channels of TDHFcntm.

Same physics and conditioning as solveMixedState (in-place overload), but the linear equation \( (h_{\rm HF} - \en_0)\,\delta F = -F_S \) is solved by Anderson mixing (anderson_coefficients) of the preconditioned fixed-point map instead of damped iteration. At the high frequencies of ionisation a bound channel's \( \en_0 = \en_a + \omega \) can sit near a spurious eigenvalue of the local preconditioner \( (h_{\rm local} + U_x - \en_0) \), where damped iteration diverges although \( (h_{\rm HF} - \en_0) \) is non-singular; Anderson mixing converges on the conditioning of the true operator. Parameters as solveMixedState.

◆ solveContinuumMixedState()

void ExternalField::solveContinuumMixedState	(	DiracSpinor *	phi,
		DiracSpinor *	Freg,
		DiracSpinor *	Firr,
		double *	K,
		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 DiracSpinor *const	Fhole = `nullptr`
	)

Continuum (en_+ > 0) mixed-state solve with the standing-wave boundary condition and the non-local exchange iterated in the source: the open channels of TDHFcntm.

Solves

\[ (h_r^{(\kappa)} + V^{\rm nl} - \en_+)\,\varphi = -F_S , \qquad \en_+ = \en_a + \omega > 0 , \]

with \( \varphi \to K\,F_{\rm irr} \) at large r, by outward integration plus F_reg subtraction (DiracODE::solveContinuumForward) on the two real homogeneous solutions at en_+: the regular, energy-normalised continuum orbital F_reg and its irregular partner F_irr (DiracODE::solveContinuumIrregular), in the local potential \( v = v_l + U_x \), with \( U_x \) the orbital-independent Kohn-Sham exchange (HF::vex_KS) for conditioning (cancelled in the source at the fixed point). The non-local exchange remainder is iterated in the source with Anderson mixing (anderson_coefficients); K is linear in phi and is mixed with the same coefficients. After each solve, phi is orthogonalised to Fa when the channel shares its kappa (norm conservation, as the bound solver); other occupied components are kept, since their dV contributions cancel pairwise across channels.

The hole-particle treatment is the caller's: pass vl already adjusted ( \( v_l - y^0_{aa} \)) and the hole orbital as Fhole, so that the iterated exchange is \( V^{\rm exch} - X_a \) (HF::vexFa_1el); together these put the ejected electron in the V^{N-1} Hamiltonian (do both or neither).

If Freg already holds the pair at en_+ (Freg->en() == en_+, nonzero), the homogeneous solutions are reused (they depend only on the channel and omega); otherwise they are built and written.

Parameters

phi	In/out: the correction orbital (kappa = channel kappa, as `Fs`). Used as the starting guess if nonzero.
Freg	In/out: regular energy-normalised continuum orbital at en_+.
Firr	In/out: irregular partner.
K	Output (if non-null): standing-wave amplitude, phi -> K F_irr.
Fa	Core orbital phi_a (energy and grid).
omega	Frequency (en_+ = Fa.en() + omega must be > 0).
vl	Local potential (including any hole-particle term).
alpha	Fine-structure constant.
core	Core orbitals (exchange).
Fs	Source spinor F_S (sign as solveMixedState: +h phi_a).
eps_target	Convergence goal of the exchange iteration.
Fhole	Optional hole orbital (V^{N-1} exchange part).

Warning: Requires en_+ > 0 and a grid dense enough at large r; a reused pair must have been built in the same local potential (not checked).

◆ 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

◆ make_rpa()

◆ ParseMethod()

◆ conditioning_states()

◆ solveMixedState() [1/4]

◆ solveMixedState() [2/4]

◆ anderson_coefficients() [1/2]

◆ anderson_coefficients() [2/2]

◆ solveMixedState_cntm()

◆ solveContinuumMixedState()

◆ solveMixedState() [3/4]

◆ solveMixedState() [4/4]

◆ solveMixedState_basis()