I'm currently trying to plot a graph wich describes a photoionization cross section as a function of incident photon energy for optical transition in a semiconductor for different values of the $\gamma$ factor, wich inlvolves a double integral. $$ \sigma= \left[ \left(\frac{\xi_{eff}}{\xi_{0}}\right)\frac{n_{r}}{\varepsilon}\right]\alpha_{fs}h\nu\sum_{f}\left|\left<\psi_{i}\right|\overrightarrow{r}\left|\psi_{f}\right>\right|^2 \delta\left(E_{f} - E_{i}-h\nu\right)$$

where $$\delta\left(E_{f} - E_{i}-h\nu\right)=\frac{1}{\pi}\frac{\hbar\Gamma_{f}}{\left(E_{f} - E_{i}-h\nu\right)²+\left(\hbar\Gamma_{f}\right)²}$$

and the wave function is given by $$\psi_{nm}(r,\theta)=\frac{1}{\sqrt{2\pi}}\sqrt{\frac{\Gamma(n+1)}{2^{\beta}\Gamma(\beta+n+1)}}\left(\gamma\beta\right)^{\beta}e^{-\frac{1}{4}\gamma²\rho²}L_{n}^{\beta}\left(\frac{1}{2}\gamma²\rho²\right)$$

where $$\beta=\sqrt{m-\Phi+\gamma/4}$$

To calculete the matrix element, I chose the states $\psi_{00}$ and $\psi_{01}$. Here's my attempt:

from scipy.integrate import nquad
import numpy as np
from scipy.special import genlaguerre, gamma
from scipy.constants import alpha
import matplotlib.pyplot as plt
import cmath
epsilon = 13.1 #dielectric constant of the material
gamma_C = 0.5 # donor impurity linewidth 
nr = 3.2 #refractive index of semiconductor
flux = 0  # Phi in eqn 8 magnetic flux
R = 5.0  #radius of the quantum ring in nm
r = np.linspace(0, 6 * R)
rho = r / R
m_0 = 0.0067*0.511 # electron effective mass
h = 4.13e-15  # Planck constant in eV
hbar =  6.58e-16  # reduced Planck constant in eV
#Photon energy
hnu = np.linspace(0, 100) #in eV

#Function that calculates the integrand
def func(rho, theta):
    betai = gama**2/2
    betaf = np.sqrt(1+gama**4/2)
    return R/np.pi*((gama * rho)**(betai + betaf) *
            np.exp(-1/2*(gama * rho)**2) *
          (gama*rho)**2/2  *np.cos(theta) * cmath.exp(-1j*theta))

def cross_section(hnu, gama):
    #function that calculates the photoionisation cross section
    betai = np.sqrt( gama**4/4)
    betaf = np.sqrt(1+gama**4/2)
    Ei = gama**2*(1+betai)-gama**4/2
    Ef = gama**2*(3+betaf)-gama**4/2
    delta = hbar * gamma_C/(Ef - Ei - hnu)**2 + ( hbar * gamma_C)**2    
    return (nr/epsilon * 4*np.pi/3 * alpha * hnu *
            (abs( np.sqrt(1/2**betai*gamma(betai + 1))*
            np.sqrt(gamma(2)/2**betaf*gamma(betaf + 2)) *
            nquad(func, [[0, np.infty], [0, np.pi]])[0])**2) 


for gama in [1.0, 1.5, 2.0]:
    plt.plot(hnu, np.real(cross_section(hnu, gama)))
    plt.plot(hnu, np.imag(cross_section(hnu, gama)))
plt.legend(['$\gamma = 1.0$', '$\gamma = 1.5$', '$\gamma = 2.0$'] )
plt.ylabel('Photoionization cross\n section $\sigma (10^{-14}cm^{2}$)')
plt.xlabel('Photon energy $h\\nu (meV)$ ')    

The plot lines did not come out as expected. I think that my problem is twofold: first, $\beta$, wich is a parameter of the wave, function vary for to each value $m$ and I dont know how to implement this correctly; secondly, there is a complex angular phase and I believe dblquad can solve only real integrals.

Any help?

  • $\begingroup$ I suppose plt.legend(['$\gamma = 1.0$ (real)', '$\gamma = 1.0$ (imag)', ETC...] ) would be more correct. But check out the label parameter of plot(...). $\endgroup$ – user66081 May 20 at 20:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.