AXWP: Attosecond X-Ray Wave Packet - v1.0

This code simulates the electronic dynamics in a neutral molecule triggered by an attosecond X-Ray pulse. The code wraps around the Qchem quantum chemistry package, which is used to calculate the molecular electronic structure and its properties.

Theoretical model

The molecular state $\ket{\Psi(t)}$ is expanded as

$$\ket{\Psi(t)} = (1+ag^{(2)}(t))\ket{\psig} + \sumc ac^{(1)}(t)e^{-i(\omegac-i\frac{\Gammav}{2})t}\ket{\psic} + \sumv av^{(2)}(t)e^{-i(\omegav-i\frac{\Gammav}{2})t}\ket{\psiv}$$.

Here, $\ket{\psig}$, $\ket{\psic}$ and $\ket{\psi_v}$ represent the ground, core-excited and valence-excited states, respectively. The $\Gamma$ terms in the exponents correspond to the states' decay rates, which are inversely proportional to their lifetime. The expressions for the probability amplitudes 'a(t)' are derived within time-dependent perturbation theory. As indicated by the superscripts, the core amplitudes are approximated at the first perturbative order, while the valence amplitudes at the second perturbative order. The probability amplitudes are combined to obtained the density matrix, which is the final output of the code.

Numerical implementation

The probability amplitudes $a_c^{(1)}(t)$ are calculated by computing the expression

Equation 1: core-excited state probability amplitudes.

The excitation energy, $\omega_{cg}$, and the transition dipole moment, $\mu_{cg}$, are calculated by the quantum chemistry package. The integration is performed numerically on a frequency array, choosing $\Gamma_c$ according to the selected excitation edge. The probability amplitudes $a_k^{(2)}(t)$, with $k\in [g,v]$, are calculated according to the expression

Equation 2: valence-/ground-state probability amplitudes.

Here, the term $M_{kg}(\omega_p,\omega_d)$ represents the RIXS transition moment, obtained at the EOM-CC level of theory. Qchem considers only the $\omega_p$ dependence $M_{kg}(\omega_p,\omega_d)$, neglecting its linear dependence on $\omega_d$. Qchem calculates $M_{kg}(\omega_p,\omega_d)$ over an interval of $\omega_p$, fixing $\omega_d$ to the first value of the $\omega_p$ vector. The code interpolates the $\omega_p$-dependent data from Qchem; the interpolating function is used to generate $M_{kg}(\omega_p,\omega_d)$ over a square frequency grid of choice, with the $\omega_p$-dependent array being 'repeated' for each value of $\omega_d$.

Figure 1: frequency grid division in terms of integration style: grey = numerical; blue = analytical. The red line indicates the strip's center, where the resonance condition $\omega_{kg}-\omega_p+\omega_d = 0$ is satisfied. The green segment indicates the discrete "strip element" the strip is divided into.

Performing the double integral in eq.2 numerically is very difficult, either because, when k=v, and $\Gamma_v$ is small and requires a large grid, or because, when k=g, $\Gamma_g=0$ and the integral diverges. The code tackles this problem by integrating numerically only in the 'grey' areas of the frequency grid in fig.1 - where $\omega-\omega_p\ll 0$ or $\gg 0$ - where the properties of the integrand allow to use a computationally affordable grid. The regions where $\omega-\omega_p\approx 0$ (the blue 'strip' in fig.1), are divided into discrete "strip elements" (the green segments in fig. 1), where the integral in eq.2 is approximated in the limit $\Gamma_k\to 0$ as

Eq. 3: valence-/ground-state probability amplitude approximation on strip element.

Here, the function $f(\omega _0)$ corresponds to the terms highlighted in red in equation 2, which is Taylor-expanded up to first order. $\omega_0$ corresponds to the centre of the strip element (i.e. intersection of green and red lines in fig. 1), while $\delta$ to the strip's half width. The code calculates the contributions from all the discrete strip elements, summing them together to obtain an approximation of the integral in equation 2 on the "blue strip". The contributions from the strip are then summed to those obtained from the numerical integration over the "grey parts" of the grid, leading to the approximate value of $a_k^{(2)}(t)$.

Implementation of the "strip" approximation

The solution in eq.3 is implemented on each strip element, whose number depends on the frequency grid chosen by the user. The derivative at $\omega0$ is estimated numerically as the difference quotient involving the upper and lower limits of the strip. To converge to the correct qualitative behaviour of $ak^{(2)}(t)$, the parameter $\delta$, which controls the width of "strip", and the finesse of the frequency grid need to be optimised. The correct behaviour of $ak^{(2)}(t)$ can be estimated by monitoring the corresponding populations, $|ak^{(2)}(t)|^2$. An example of "non-converged" behaviour is shown in figure 2, where in some cases $|ak^{(2)}(t)|^2$ diverges to infinity is are non-zero before the pulse.

Figure 2: $|ak^{(2)}(t)|^2$ in a "non-converged" calculation for a pulse centred around 0 as. Here, the chosen $\delta$ and frequency grid finesse are not ideal.

A "converged" behaviour can be obtained by fine-tuning the $\delta$ and frequency grid size until a result similar to fig. 3 is obtained. Here, $|ak^{(2)}(t)|^2$ shows the qualitatively correct behaviour, not diverging to infinity and being zero before the pulse arrives at 0 as.

Figure 2: $|ak^{(2)}(t)|^2$ in a "converged" calculation for a pulse centred around 0 as. Here, the chosen $\delta$ and frequency grid finesse are the ideal ones.

The current implementation imposes the same $\delta$ (set inside the code) and frequency grid size (input to the script) to approximate eq. 2 for all the considered ground/valence-excited states. In testing, this led in some case to a slight overestimation/underestimation of $|a_k^{(2)}(t)|^2$. However, the correct qualitative behaviour is achieved across all the considered molecules, excitation edges and pulse parameters examined.

Improvements - in progress

1. Implement selection of $\delta$ and grid step size for groung and each valence-excited state.
2. Add higher orders to Taylor expansion in eq. 3.
3. Implement option for analytical derivative (for Gaussian pulses).

Dependencies

The following packages are required:

- multiprocessing
- cmath, math, time
- numpy, scipy, matplotlib
- itertools
- numexpr

Calculation setup

A Qchem calculation with the sample SRIXS and XAS input files must be performed. The output files must end with '-XAS.out' and '-SRIXS.out', respectively. The wave packet calculation can be started by performing

python3 wp_calc.py 'path_Qchem_out' '#colors' C1_en C2_en C1_bw C2_bw 'C1_pol' 'C2_pol' #freq_grid 'output_array'

The keys have the corresponding meaning

1. path_Qchem_out: name of output files before ending suffix
2. #colors: '1C'= 1 color Gaussian, '2C' = 2 color Gaussian
3. C1_en/C2_en: color1/color2 central energy (eV)
4. C1_bw/C2_bw: color1/color2 bandwidth (eV)
5. C1_pol/C2_pol: color1/color2 polarisation. Options are:
   • linear polarised: 'x', 'z', 'xz'
   • circularly polarised: 'Rxy', 'Lxy'
6. #freq_grid: number of points frequecy grid with shape (#freq_grid,#freq_grid)
7. output_array: path+name output

Output contains python dictionary with entries:

• 'Density Matrix': density matrix of shape=(#states,#states), where #states is ordered as gs, valence-excited, core-excited.
• '1PDM': 1PDMs array of shape=(#states,#states,#MO,#MO); state 1PDM diagonal axes 0,1; transition 1PDM off-diagonal axes 0,1.
• 'time_array': array of considered time points.
• 'pulse_time': pulse in the time domain.
• '#_val_states': number valence-excited states + 1 (gs).
• '#_core_states': number of core-excited states.

Improvements - in progress

• Parameters input via input file:
  • Section for pulse parameters
  • Section for integration and strip parameters
• Separate module for pulses