# How to plot a function for multiple values of a parameter in the same set of axes in Python?

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
#Constants
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)) *
*delta)

#Plot
plt.figure();plt.clf()

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?

• 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(...). – user66081 May 20 at 20:15