Algorithms

Algorithms#

We need to compute

\[ I \propto \frac{1}{\mathcal{Z}(T)}\sum_{i}e^{- E_{i}/(k_\mathrm{B}T)} \sum_f |M_{fi}|^2 \delta(E_f + \hbar\omega_{\boldsymbol{k}^\prime} - E_i - \hbar\omega_{\boldsymbol{k}} ) \]

\[ M_{fi} = \sum_n \frac{\bra{f} {\cal D}^\dagger_{\boldsymbol{k}^\prime\hat{\epsilon}^\prime}\ket{n}\bra{n} {\cal D}^{\phantom\dagger}_{\boldsymbol{k}\hat{\epsilon}}\ket{i}}{E_n - E_i - \hbar\omega_{\boldsymbol{k}}+\mathrm{i}\Gamma_c}, \]

where \(\delta\) is simulated with a Lorentzian of half width at half maximum \(\Gamma\).

Direct solution#

EDRIXS hosts a pure Python implementation of the cross-section that is suitable for performance insensitive cases and testing. This was what we used in the atomic model section. See the diagonalization step in edrixs.solvers.ed_1v1c_py and spectrum construction in edrixs.solvers.rixs_1v1c_py. The following is a crude illustration of what happens

Build emat and umat in \(Y_l^m\) single-particle basis

import edrixs
import scipy
import numpy as np

ten_dq = 1.0
emat_i = np.zeros((14, 14), dtype=complex)
emat_i[:10, :10] = edrixs.cf_cubic_d(ten_dq) # plus SOC and hopping

F0, F2, F4 = 4.0, 1.0, 0.3
umat_i = np.zeros((14, 14, 14, 14), dtype=complex)
umat_i = edrixs.get_umat_slater('d', F0, F2, F4) # plus core hole terms 

# For L3-edge RIXS of d-shell material
print(f"emat is 14x14")
print(f"umat is 14x14x14x14")

emat is 14x14
umat is 14x14x14x14

Construct Fock basis for \(\ket{i}\) and \(\ket{n}\).

v_norb = 10
v_noccu = 8
c_norb = 4
basis_i = edrixs.get_fock_bin_by_N(v_norb, v_noccu, c_norb, c_norb)
basis_n = edrixs.get_fock_bin_by_N(v_norb, v_noccu+1, c_norb, c_norb - 1)
print(np.array(basis_i)[:3, :10])

[[1 1 1 1 1 1 1 1 0 0]
 [1 1 1 1 1 1 1 0 1 0]
 [1 1 1 1 1 1 1 0 0 1]]

Build Hamiltonian for initial and intermediate state using edrixs.two_fermion and edrixs.four_fermion.

ncfg_i, ncfg_n = len(basis_i), len(basis_n)
print(f"{ncfg_i} ground state vectors")
print(f"{ncfg_n} intermediate state vectors")
hmat_i = np.zeros((ncfg_i, ncfg_i), dtype=complex)
hmat_n = np.zeros((ncfg_n, ncfg_n), dtype=complex)
hmat_i[:, :] += edrixs.two_fermion(emat_i, basis_i, basis_i)
hmat_i[:, :] += edrixs.four_fermion(umat_i, basis_i)
# hmat_n left out for brevity. 

45 ground state vectors
40 intermediate state vectors

Diagonalize Hamiltonian for initial and intermediate state using scipy, which calls LAPACK under the hood.

eval_i, evec_i = scipy.linalg.eigh(hmat_i)
eval_n, evec_n = scipy.linalg.eigh(hmat_n)

Load relevant dipole matrix operators

trans_op_Ylm = edrixs.get_trans_oper('dp32')

Build the many-body transition operators for each polarization and express them in the eigenbasis

npol = 3
trans_op = np.zeros((npol, ncfg_n, ncfg_i), dtype=complex)
for i in range(npol):
    trans_op[i] = edrixs.cb_op2(edrixs.two_fermion(trans_op_Ylm[i], basis_n, basis_i),
                                evec_n, evec_i)
print(trans_op.shape)

(3, 40, 45)

RIXS spectra can be built from trans_op, eval_i, and eval_n.

Current EDRIXS solver#

Currently, EDRIXS calls into its Fortran layer. This is sparsely documented and suffers from serious limitations in portability, maintainability, and flexibility. See edrixs.ed_siam_fort and edrixs.rixs_siam_fort, which calls ed_driver.f90 and rixs_driver.f90. We used this in the Anderson impurity model section.

Mathematics#

To do the calculation, it is useful to reformulate the cross-section compared to the most common Kramers-Heisenberg form. We first note that one can re-express \(M_{fi}\) as

\[ M_{fi} = \bra{f} {\cal D}^\dagger_{\boldsymbol{k}^\prime\hat{\epsilon}^\prime} \frac{1}{{\cal \widetilde{H}} - E_i - \hbar\omega_{\boldsymbol{k}}+i\Gamma_c} {\cal D}^{\phantom\dagger}_{\boldsymbol{k}\hat{\epsilon}}\ket{i}, \]

where \({\cal \widetilde{H}}\) is the intermediate state Hamiltonian with \(\widetilde{\mathcal H}\ket{n} = E_n\ket{n}\) and \(\sum_n \ket{n}\bra{n} = \mathbb{I}.\) To see this, note that using a spectral decomposition we can write

\[ \frac{1}{{\cal \widetilde{H}} - E_i - \hbar\omega_{\boldsymbol{k}}+i\Gamma_c } = \sum_n |n\rangle \frac{1}{E_n - E_i - \hbar\omega_{\boldsymbol{k}}+i\Gamma_c }\langle n|. \]

Inserting this operator between the dipole operators, we get

\[ M_{fi} = \bra{f}{\cal D}^\dagger_{\mathbf{k}'\hat\epsilon'} \Biggl(\sum_n |n\rangle \frac{1}{E_n-E_i-\hbar\omega_{\mathbf{k}}+i\Gamma_c}\langle n|\Biggl) {\cal D}_{\mathbf{k}\hat\epsilon}\ket{i}. \]

By rearranging, we can show that this is equivalent to the original equation with

\[ M_{fi} =\sum_n \frac{\langle f|{\cal D}^\dagger_{\mathbf{k}'\hat\epsilon'}|n\rangle \langle n|{\cal D}_{\mathbf{k}\hat\epsilon}|i\rangle}{E_n-E_i-\hbar\omega{\mathbf{k}}+i\Gamma_c}. \]

The next step is to define an effective scattering operator

\[ \mathcal O_i \equiv \mathcal D^\dagger_{\boldsymbol{k}^\prime\hat\epsilon^\prime} \frac{1}{\widetilde{\mathcal H}-E_i-\hbar\omega_{\boldsymbol k}+i\Gamma_c} \mathcal D_{\boldsymbol{k}\hat\epsilon}, \]

so that we have

\[ \sum_f |M_{fi}|^2 \delta(E_f-E_i-\hbar\omega_{\boldsymbol q})=\sum_f \langle i|\mathcal O_i^\dagger|f\rangle\langle f|\mathcal O_i|i\rangle \delta(E_f-E_i-\hbar\omega_{\boldsymbol q}). \]

Since the set of \(|f\rangle\) is mutually orthonormal, they have the property that

\[ \sum_f |f\rangle \delta(E_f-E)\langle f|=\delta(\mathcal H-E) \]

such that

\[ \sum_f |M_{fi}|^2 \delta(E_f-E_i-\hbar\omega_{\boldsymbol q})=\langle i|\mathcal O_i^\dagger \delta({\mathcal H} -E_i-\hbar\omega_{\boldsymbol q}) \mathcal O_i|i\rangle . \]

The Dirac \(\delta\) function can be represented as a limiting case of a Lorentzian with width \(\Gamma\)

\[ \delta({\mathcal H} -E_i-\hbar\omega_{\boldsymbol q})\to -\frac{1}{\pi} \Im\left(\frac{1}{{\mathcal H} -E_i-\hbar\omega_{\boldsymbol q} +i \Gamma}\right). \]

If we define the prepared state using \(|F_i\rangle \equiv \mathcal O_i|i\rangle\), and leave out a factor \(\pi\), we can assemable the full equation as

\[ I \propto -\sum_i e^{-E_i/(k_BT)} \Im\langle F_i|\frac{1}{\mathcal H-E_i-\hbar\omega_{\boldsymbol q}+i\Gamma}|F_i\rangle. \]

In this form, it is possible to effciently compute the spectrum using Lanczos methods. Let’s define

\[ \sigma_i = \langle F_i|\frac{1}{\mathcal H-E_i-\hbar\omega_{\boldsymbol q}+i\Gamma}|F_i\rangle \]

and

\[ z = E_i + \hbar\omega_{\boldsymbol q} - i\Gamma. \]

The key thing to compute is

\[ \sigma_i = |F_i|^2 \langle f_i|\frac{1}{\mathcal H- z}|f_i\rangle, \]

where \(\ket{F_i} = |F_i| \ket{f_i}\) splits the vectors into a nomalized component and a normalization factor.

Krylov’s methods give us a procedure for constructing a new basis \(V_m = [\ket{q_0}, \ket{q_1}, ... \ket{q_m}]\), whose first element is \(\ket{q_0} = \ket{F_i}\), which allows us to effciently represent the Hamiltonanian in a tridiagonal form

\[ {\cal H} \approx V_m T_m V_m^\dagger \]

\[\begin{split} T_m = \begin{pmatrix} \alpha_0 & \beta_1 & 0 & 0 & \cdots & 0 \\ \beta_1 & \alpha_1 & \beta_2 & 0 & \cdots & 0 \\ 0 & \beta_2 & \alpha_2 & \beta_3 & \cdots & 0 \\ 0 & 0 & \beta_3 & \alpha_3 & \ddots & \vdots \\ \vdots & \vdots & \vdots & \ddots & \ddots & \beta_{m-1} \\ 0 & 0 & 0 & \cdots & \beta_{m-1} & \alpha_{m-1} \end{pmatrix}. \end{split}\]

Once this form is constructed, the components of the spectrum can be built via.

\[ \sigma_i = |F_i|^2\,\langle f_i|(\mathcal H - z)^{-1}|f_i\rangle \approx |F_i|^2\, \cfrac{1}{ \alpha_1 - z - \cfrac{\beta_1^2}{ \alpha_2 - z - \cfrac{\beta_2^2}{ \alpha_3 - z - \ddots - \cfrac{\beta_{m-1}^2}{ \alpha_m - z } } }} \]

The recursive continued fraction arises because the inverse of a tridiagonal matrix can be constructed by repeatedly eliminating deeper blocks, and Lanczos produces exactly the tridiagonal structure required for this elimination to proceed step by step.

Computational methods#

The general process is:

Fock basis is generated in integer representation

for binary_string in ["011", "101", "110"]:
    print(f"{binary_string} is {int(binary_string, 2)}")

is 3
is 5
is 6

Hamiltonian is constructed in CSR format. For example, the matrix

\[\begin{split} \begin{bmatrix} 10 & 20 & 0 & 0 & 0 & 0 \\ 0 & 30 & 0 & 40 & 0 & 0 \\ 0 & 0 & 50 & 60 & 70 & 0 \\ 0 & 0 & 0 & 0 & 0 & 80 \end{bmatrix} \end{split}\]

is encoded as

V         = [10, 20, 30, 40, 50, 60, 70, 80]
COL_INDEX = [0,  1,  1,  3,  2,  3,  4,  5]
ROW_INDEX = [0,  2,  4,  7,  8]

A handful of the eigenvectors of the ground state Hamiltonian \({\cal H}\) are obtained via the Lanczos method within the parallel ARPACK library.

\[ {\cal H} \ket{i} = E_i \ket{i} \]

Apply the absorption transition operator to the ground state

\[ \ket{b_i} = {\cal D}_{\boldsymbol{k},\hat{\epsilon}} \ket{i}. \]

Solve the following linear equation, involving the intermediate state Hamiltonian \({\cal \widetilde{H}}\) via sparse MINRES methods

\[ \left({\cal \widetilde{H}} - E_i - \hbar\omega_{\boldsymbol{k}}+i\Gamma_c\right) \ket{x_i} = \ket{b_i} \]

Apply emission operator

\[ \ket{F_i} = {\cal D}^\dagger_{\boldsymbol{k}^\prime,\hat{\epsilon}^\prime} \ket{x_i} \]

The spectrum can then be represented as

\[ I \propto -\sum_{i}e^{- E_{i}/(k_\mathrm{B}T)} \Im \bra{F_i} \frac{1}{{\cal H} - E_i - \hbar\omega_{\boldsymbol{q}}+i \Gamma} \ket{F_i} \]

The continued fraction technique is then used to construct the spectrum

Plans#

We are still finalizing the planned approach.

Algorithm changes#

Improved basis lookup methods for constructing Hamiltonian based on Ref. [1].
Implement methods to re-use the shared structure of Krylov subspaces for different incident x-ray energies based on Ref. [2]. Kipton recommends this preprint [3].

Infrastructure changes#

The current plan is to call into more modern platform-agnostic CPU/GPU packages PETSc+SLEPc. These, which are primarily written in C, have Python interfaces.

Algorithms

Contents