# -*- coding: utf-8 -*-
r"""
ePSproc MFPAD functions
-----------------------
29/10/20 v2 Updated with Wigner delay routine.
Tidied up To Do list items.
Tidied up docstring.
05/08/19 v1 Initial python version.
Based on original Matlab code ePS_MFPAD.m
Structure
---------
Calculate MFPAD on a grid from input ePS matrix elements.
Use fast functions, pre-calculate if possible.
Data in Xarray, use selection functions and multiplications based on relevant quantum numbers, other axes summed over.
Choices for functions...
* `Moble's spherical functions (quaternion based) <https://github.com/moble/spherical_functions>`_
Provides fast wignerD, 3j and Ylm functions, uses Numba.
Install with:
>>> conda install -c conda-forge spherical_functions
* `Scipy special.sph_harm <https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.sph_harm.html#scipy.special.sph_harm>`_
To Do
-----
* More efficient computation? Use Xarray group by?
Formalism
---------
MF-PADs numerical expansion
============================
This is implemented numerically in :py:func:`epsproc.MFPAD.mfpad()`. For direct computation of $\beta_{LM}$ parameters (rather than full observable expansion) see :py:func:`epsproc.MFBLM` (original looped version) or :py:func:`epsproc.geomFunc.mfblmGeom` (geometric tensor version).
From `ePSproc: Post-processing suite for ePolyScat electron-molecule scattering calculations <https://www.authorea.com/users/71114/articles/122402/_show_article>`_.
.. math::
I_{\mu_{0}}(\theta_{\hat{k}},\phi_{\hat{k}},\theta_{\hat{n}},\phi_{\hat{n}})=\frac{4\pi^{2}E}{cg_{p_{i}}}\sum_{\mu_{i},\mu_{f}}|T_{\mu_{0}}^{p_{i}\mu_{i},p_{f}\mu_{f}}(\theta_{\hat{k}},\phi_{\hat{k}},\theta_{\hat{n}},\phi_{\hat{n}})|^{2}\label{eq:MFPAD}
T_{\mu_{0}}^{p_{i}\mu_{i},p_{f}\mu_{f}}(\theta_{\hat{k}},\phi_{\hat{k}},\theta_{\hat{n}},\phi_{\hat{n}})=\sum_{l,m,\mu}I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)Y_{lm}^{*}(\theta_{\hat{k}},\phi_{\hat{k}})D_{\mu,\mu_{0}}^{1}(R_{\hat{n}})\label{eq:TMF}
I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)=\langle\Psi_{i}^{p_{i},\mu_{i}}|\hat{d_{\mu}}|\Psi_{f}^{p_{f},\mu_{f}}\varphi_{klm}^{(-)}\rangle\label{eq:I}
In this formalism:
* :math:`I_{l,m,\mu}^{p_{i}\mu_{i},p_{f}\mu_{f}}(E)` is the radial part of the dipole matrix element, determined from the initial and final state electronic wavefunctions :math:`\Psi_{i}^{p_{i},\mu_{i}}` and :math:`\Psi_{f}^{p_{f},\mu_{f}}`, photoelectron wavefunction :math:`\varphi_{klm}^{(-)}` and dipole operator :math:`\hat{d_{\mu}}`. Here the wavefunctions are indexed by irreducible representation (i.e. symmetry) by the labels :math:`p_{i}` and :math:`p_{f}`, with components :math:`\mu_{i}` and :math:`\mu_{f}` respectively; :math:`l,m` are angular momentum components, :math:`\mu` is the projection of the polarization into the MF (from a value :math:`\mu_{0}` in the LF). Each energy and irreducible representation corresponds to a calculation in ePolyScat.
* :math:`T_{\mu_{0}}^{p_{i}\mu_{i},p_{f}\mu_{f}}(\theta_{\hat{k}},\phi_{\hat{k}},\theta_{\hat{n}},\phi_{\hat{n}})` is the full matrix element (expanded in polar coordinates) in the MF, where :math:`\hat{k}` denotes the direction of the photoelectron :math:`\mathbf{k}`-vector, and :math:`\hat{n}` the direction of the polarization vector :math:`\mathbf{n}` of the ionizing light. Note that the summation over components :math:`\{l,m,\mu\}` is coherent, and hence phase sensitive.
* :math:`Y_{lm}^{*}(\theta_{\hat{k}},\phi_{\hat{k}})` is a spherical harmonic.
* :math:`D_{\mu,\mu_{0}}^{1}(R_{\hat{n}})` is a Wigner rotation matrix element, with a set of Euler angles :math:`R_{\hat{n}}=(\phi_{\hat{n}},\theta_{\hat{n}},\chi_{\hat{n}})`, which rotates/projects the polarization into the MF.
* :math:`I_{\mu_{0}}(\theta_{\hat{k}},\phi_{\hat{k}},\theta_{\hat{n}},\phi_{\hat{n}})` is the final (observable) MFPAD, for a polarization :math:`\mu_{0}` and summed over all symmetry components of the initial and final states, :math:`\mu_{i}` and :math:`\mu_{f}`. Note that this sum can be expressed as an incoherent summation, since these components are (by definition) orthogonal.
* :math:`g_{p_{i}}` is the degeneracy of the state :math:`p_{i}`.
For more details see: Toffoli, D., Lucchese, R. R., Lebech, M., Houver, J. C., & Dowek, D. (2007). Molecular frame and recoil frame photoelectron angular distributions from dissociative photoionization of NO2. The Journal of Chemical Physics, 126(5), 054307. https://doi.org/10.1063/1.2432124
(MF) Wigner Delays
==================
This is implemented numerically in :py:func:`epsproc.MFPAD.mfWignerDelay()`.
The partial wave, or channel resolved, Wigner delay can be given as the energy-derivative of the continuum wavefunction phase, $\arg(\psi_{lm}^{*})$
.. math::
\begin{equation}
\tau_{w}(k,\theta,\phi)=\hbar\frac{d\arg(\psi_{lm}^{*})}{d\epsilon}
\end{equation}
And the total, or group, delay from the (coherent) sum over these channels.
.. math::
\begin{equation}
\tau_{w}^{g}(k,\theta,\phi)=\hbar\frac{d\arg(\sum_{l,m}\psi_{lm}^{*})}{d\epsilon}\label{eq:tauW_mol_sum}
\end{equation}
Here the full wavefunction $\sum_{l,m}\psi_{lm}^{*}$ is essentially equivalent to the phase of the (conjugate) $T_{\mu_{0}}^{p_{i}\mu_{i},p_{f}\mu_{f}*}(\theta_{\hat{k}},\phi_{\hat{k}},\theta_{\hat{n}},\phi_{\hat{n}})$ functions given above (and by :py:func:`epsproc.MFPAD.mfpad()`).
For more details see, for example:
- Hockett, Paul, Eugene Frumker, David M Villeneuve, and Paul B Corkum. “Time Delay in Molecular Photoionization.” Journal of Physics B: Atomic, Molecular and Optical Physics 49, no. 9 (May 2016): 095602. https://doi.org/10.1088/0953-4075/49/9/095602.
- Carvalho, C.A.A. de, and H.M. Nussenzveig. “Time Delay.” Physics Reports 364, no. 2 (June 2002): 83–174. https://doi.org/10.1016/S0370-1573(01)00092-8.
- Wigner, Eugene. “Lower Limit for the Energy Derivative of the Scattering Phase Shift.” Physical Review 98, no. 1 (April 1955): 145–47. https://doi.org/10.1103/PhysRev.98.145.
"""
# Imports
import numpy as np
import pandas as pd
import xarray as xr
# Special functions
# from scipy.special import sph_harm
import spherical_functions as sf
import quaternion
import scipy.constants
# Package fns.
# TODO: tidy this up!
from epsproc.util import matEleSelector
from epsproc.sphCalc import sphCalc, setPolGeoms
from epsproc.sphPlot import plotTypeSelector
[docs]def mfpad(dataIn, thres = 1e-2, inds = {'Type':'L','it':1}, res = 50, R = None, p = 0):
"""
Parameters
----------
dataIn : Xarray
Contains set(s) of matrix elements to use, as output by epsproc.readMatEle().
thres : float, optional, default 1e-2
Threshold value for matrix elements to use in calculation.
ind : dictionary, optional.
Used for sub-selection of matrix elements from Xarrays.
Default set for length gauage, single it component only, inds = {'Type':'L','it':'1'}.
res : int, optional, default 50
Resolution for output (theta,phi) grids.
R : list or Xarray of Euler angles or quaternions, optional.
Define LF > MF polarization geometry/rotations.
For default case (R = None), calls :py:func:`epsproc.setPolGeoms()`, where 3 geometries are calculated, corresponding to z-pol, x-pol and y-pol cases.
Defined by Euler angles (p,t,c) = [0 0 0] for z-pol, [0 pi/2 0] for x-pol, [pi/2 pi/2 0] for y-pol.
p : int, optional.
Defines LF polarization state, p = -1...1, default p = 0 (linearly pol light along z-axis).
TODO: add summation over p for multiple pol states in LF.
Returns
-------
Ta
Xarray (theta, phi, E, Sym) of MFPADs, summed over (l,m)
Tlm
Xarray (theta, phi, E, Sym, lm) of MFPAD components, expanded over all (l,m)
"""
# Define reduced data from selection over all data
daRed = matEleSelector(dataIn, thres = thres, inds = inds)
# Generate spherical harmonics
Lmax = daRed.l.max()
YlmX = sphCalc(Lmax, res = res)
# Reindex to match data (should happen automagically, but not always!)
# YlmXre = YlmX.reindex_like(daRed)
# Set rotation angles for LF > MF
if R is None:
# Set (x,y,z) projection terms only
# Nangs = 10
# pRot = np.linspace(0,180,Nangs)
# tRot = np.linspace(0,90,Nangs)
# cRot = np.linspace(0,180,Nangs)
# eAngs = np.array([pRot, tRot, cRot,])*np.pi/180
# Convert to quaternions
# R = quaternion.from_euler_angles(pRot*np.pi/180, tRot*np.pi/180, cRot*np.pi/180)
# # Eugler angles for rotation of LF->MF, set as [0 0 0] for z-pol, [0 pi/2 0] for x-pol, [pi/2 pi/2 0] for y-pol
# pRot = [0, 0, np.pi/2]
# tRot = [0, np.pi/2, np.pi/2]
# cRot = [0, 0, 0]
# eAngs = np.array([pRot, tRot, cRot]) # List form to use later
# Euler = pd.MultiIndex.from_arrays(eAngs, names = ['P','T','C'])
#
# # Convert to quaternions
# R = quaternion.from_euler_angles(pRot, tRot, cRot)
# 02/10/20 - updated to use sphCalc.setPolGeoms()
R = setPolGeoms()
#**************** Calculate MFPADs
Tlm = []
Ta = []
# Loop over pol geoms R
for n, Rcalc in enumerate(R):
T = []
# Loop over mu terms and multiply
for mu in np.arange(-1,2):
# Set by element replacement (preserves whole structure)
# daTemp = daRed.copy() # Set explicit copy for rotation.
# daTemp.loc[{'mu':mu}].values = daTemp.loc[{'mu':mu}].values * sf.Wigner_D_element(Rcalc, 1, mu, 0).conj()
# Issues with reindexing to extra coords at the moment, so reindex and multiply for specific mu only
# daTemp = daTemp.sel({'mu':mu})
# YlmXre = YlmX.reindex_like(daTemp)
# T.append(YlmXre.conj() * daTemp) # Output full (l,m,mu) expansion
# Set by looping and selection
daTemp = daRed.sel({'mu':mu}) * sf.Wigner_D_element(Rcalc.item(), 1, mu, 0).conj()
YlmXre = YlmX.reindex_like(daTemp)
T.append(YlmXre.conj() * daTemp) # Output full (l,m,mu) expansion
# Concat & sum over symmetries
Ts = xr.combine_nested([T[0], T[1], T[2]], concat_dim=['LM'])
# Add dims - currently set for Euler angles only.
# Can't seem to add mutiindex as a single element, so set dummy coord here and replace below.
Ts = Ts.expand_dims({'Euler':[n]}) # Set as index
# Ts = Ts.expand_dims({'p':[eAngs[0,n]], 't':[eAngs[1,n]], 'c':[eAngs[2,n]]})
Tlm.append(Ts)
Ta.append(Ts.sum(dim = 'LM'))
TlmX = xr.combine_nested(Tlm, concat_dim=['Euler'])
TaX = xr.combine_nested(Ta, concat_dim=['Euler'])
# Assign Euler angles to dummy dim
# TlmX = TlmX.assign_coords(Euler = Euler)
# TaX = TaX.assign_coords(Euler = Euler)
TlmX = TlmX.assign_coords(Euler = R['Euler']) # Coords from Euler Xarray.
TaX = TaX.assign_coords(Euler = R['Euler'])
# Set generic data type as TX
TlmX.attrs['dataType'] = 'TX'
TaX.attrs['dataType'] = 'TX'
return TaX, TlmX # , Ta, Tlm # For debug also return lists
[docs]def mfWignerDelay(dataIn, pType = 'phaseUW'):
"""
Compute MF Wigner Delay as phase derivative of continuum wavefunction.
Parameters
----------
dataIn : Xarray
Contains set(s) of wavefunctions to use, as output by :py:func:`epsproc.MFPAD.mfpad()`.
pType : str, optional, default = 'phaseUW'
Used to set data conversion options, as implemented in :py:func:`epsproc.plotTypeSelector()`
- 'phase' use np.angle()
- 'phaseUW' unwapped with np.unwrap(np.angle())
"""
# Set data
Tw = dataIn.copy()
Tw['EJ'] = Tw['Eke']*scipy.constants.e # Convert to J units
unitConv = scipy.constants.hbar/1e-18 # Convert to attoseconds
# TODO: move this to util functions.
# Use existing functionality to set phases - easier if unwrap required (uses lambdas)
Tw = plotTypeSelector(Tw.conj(), pType = pType).differentiate('EJ')* unitConv
# Propagate attrs
Tw.attrs = dataIn.attrs
Tw.attrs['dataType'] = 'Wigner'
Tw.attrs['units'] = 'attosecond'
return Tw