Transition operators for photon spectroscopy

In order to calculate excitation spectra one needs to define a transition operator. In principle any operator can be written in second quantization, making it possible to create many different kinds of excitation spectroscopy. In practice we very often want to calculate absorption spectroscopy due to dipole transitions. In this case one can use the crystal-field operator to quickly create the transition operator wanted.

The Hamiltonian including the interaction between the light and the electrons is given as: \begin{equation} H = H_0 + \frac{e}{m} (p \cdot A ) + \frac{e^2}{2m} A \cdot A, \end{equation} with $H_0$ the Hamiltonian of the system without the light and $A$ the vector field of the photons. Here we used the Coulomb Gauge to simplify $p \cdot A + A \cdot p = 2 p \cdot A$. There are higher order terms, there is the magnetic field of the photon and at high flux there need to be some other terms. (To be looked at !!!!). Below we discuss the two different terms separately.

The $p \cdot A$ term simplifies in the dipole approximation to $\varepsilon \cdot \hat{r}$ with $\varepsilon$ the polarization of the light and $\hat{r}$ a unit vector. The operator $\hat{r}$ can be created in Quanty by expanding this on renormalized spherical harmonics: \begin{eqnarray} \varepsilon \cdot \hat{r} &=& (\varepsilon_x + \imath \varepsilon_y) C_{k=1,m=-1}\\ \nonumber &+& \varepsilon_z C_{k=1,m=0}\\ \nonumber &+& (-\varepsilon_x + \imath \varepsilon_y) C_{k=1,m=1}. \end{eqnarray} The transition operator in Quanty for light polarized in the $\varepsilon=\{\varepsilon_x,\varepsilon_y,\varepsilon_z\}$ direction can be created as:

Example.Quanty

-- dipole operator for 
-- light polarized in the
-- {ex,ey,ez} direction
NF = 16
NB = 0
IndexDn_2p = {0,2,4}
IndexUp_2p = {1,3,5}
IndexDn_3d = {6,8,10,12,14}
IndexUp_3d = {7,9,11,13,15}
ex = 1
ey = 2
ez = 3
exn = ex / math.sqrt(ex^2 + ey^2 + ez^2)
eyn = ey / math.sqrt(ex^2 + ey^2 + ez^2)
ezn = ez / math.sqrt(ex^2 + ey^2 + ez^2)
 
Akm = {{1, -1,( exn + I * eyn)/math.sqrt(2)},
        {1,  0, ezn},
        {1,  1,(-exn + I * eyn)/math.sqrt(2)}}
OppT = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm)

A spectrum calculated using the above defined transition operator returns the complex quantity: \begin{equation} G(\omega) = \bigg\langle \psi_i \bigg| \varepsilon^* \cdot \hat{r} \frac{1}{\omega - H + \imath \Gamma/2} \varepsilon \cdot \hat{r} \bigg| \psi_i \bigg\rangle \end{equation}

For non-resonant inelastic x-ray scattering (nIXS or Raman or Compton scattering) we follow the notation and derivations as given in Phys. Rev. Lett \textbf{99}, 257401 (2007) and rewrite the intensity in terms of the $\vec{q}$ dependent structure factor, which relates to the dynamical susceptibility defined as $G(\omega)$ with $T=e^{\imath \vec{q} \cdot \vec{r}}$.\cite{Haverkort:2007bv} In order to include this we expand $T=e^{\imath \vec{q} \cdot \vec{r}}$ on renormalized spherical harmonics and spherical Bessel functions: \begin{eqnarray} \nonumber e^{\imath \vec{q}\cdot \vec{r}} &=& \sum_{k=0}^{\infty} \sum_{m=-k}^{m=k} \imath^k \, (2k+1) \, j_k(qr) \, C_{k,m}^*(\theta_q,\phi_q)\\ &&\times C_{k,m}(\theta_r,\phi_r). \end{eqnarray} If we write our one particle orbital basis functions as a product of a radial wavefunction times a spherical harmonic ($\phi_{n,l,m}(\vec{r})=R_{n,l}(r)Y_{l,m}(\theta,\phi)$) we can express the transition operator as: \begin{equation} T = \sum_{\tau_1,\tau_2} \sum_{k,m} A_{k,m} \big\langle Y_{l_1,m_1} \big| C_{k,m} \big| Y_{l_2,m_2} \big\rangle a^{\dagger}_{\tau_1} a^{\phantom{\dagger}}_{\tau_2}, \end{equation} with $A_{k,m}$ given as: \begin{equation} A_{k,m} = \imath^k \, (2k+1)\, \big\langle R_{n_1,l_1} \big| j_k(kr) C_{k,m}^*(\theta_q,\phi_q) \big| R_{n_2,l_2} \big\rangle \end{equation}

In Quanty one can define the appropriate functions to calculate these prefactors. Given the radial wavefunctions stored in a file the $q$ dependent $d-d$ excitation spectrum for Ni in NiO could be calculated as:

Example.Quanty

-- read the radial wave functions
-- they are stored in the file "NiO_Radial"
-- order of functions
-- r	1S	2S	2P	3S	3P	3D
file = io.open("NiO_Radial", "r")
Rnl = {}
for line in file:lines() do
  RnlLine={}
  for i in string.gmatch(line, "%S+") do
    table.insert(RnlLine,i)
  end
  table.insert(Rnl,RnlLine)
end
 
-- some constants
a0      =  0.52917721092
Rydberg = 13.60569253
Hartree = 2*Rydberg
 
-- dd transitions from 3d (index 7 in Rnl) to 
-- 3d (index 7 in Rnl)
-- define a function that calculates
-- <R(r) | j_k(q r) | R(r)>
function RjRdd (q)
  Rj0R = 0
  Rj2R = 0
  Rj4R = 0
  dr = Rnl[3][1]-Rnl[2][1]
  r0 = Rnl[2][1]-2*dr
  for ir = 2, #Rnl, 1 do
    r = r0 + ir * dr
    Rj0R = Rj0R + Rnl[ir][7] * 
      math.SphericalBesselJ(0,q*r) * Rnl[ir][7] * dr
    Rj2R = Rj2R + Rnl[ir][7] * 
      math.SphericalBesselJ(2,q*r) * Rnl[ir][7] * dr
    Rj4R = Rj4R + Rnl[ir][7] * 
      math.SphericalBesselJ(4,q*r) * Rnl[ir][7] * dr
  end
  return Rj0R, Rj2R, Rj4R
end
 
-- the angular part is given as 
-- C(theta_q, phi_q)^* C(theta_r, phi_r)
-- which is a potential expanded on 
-- spherical harmonics
function ExpandOnClm(k,theta,phi,scale)
  ret={}
  for m=-k, k, 1 do
    table.insert(ret,{k,m,scale * math.SphericalHarmonicC(k,m,theta,phi)})
  end
  return ret
end
 
-- define nIXS transition operators
function TnIXS_dd(q, theta, phi)
  Rj0R, Rj2R, Rj4R = RjRdd(q)
  k=0
  A0 = ExpandOnClm(k, theta, phi, I^k*(2*k+1)*Rj0R)
  T0 = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, A0)
  k=2
  A2 = ExpandOnClm(k, theta, phi, I^k*(2*k+1)*Rj2R)
  T2 = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, A2)
  k=4
  A4 = ExpandOnClm(k, theta, phi, I^k*(2*k+1)*Rj4R)
  T4 = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, A4)
  T = T0+T2+T4
  T.Chop()
  return T
end
 
-- q in units per a0 (if you want in units per A 
-- take 5*a0 to have a q of 5 per A)
q=4.5
 
print("for q=",q," per a0 (",q / a0," per A) The ratio of k=0, k=2 and k=4 transition strength is:",RjRdd(q))
 
-- define some transition operators
qtheta=0
qphi=0
T = TnIXS_dd(q,qtheta,qphi)
 
-- calculate the spectra
nIXSSpectra = CreateSpectra(H, T, psi)