Optical properties of metals

Optical properties of metals

In this notebook, let us evaluate the optical properties of metals by modelling the metal as consisting of sea of electrons that are not bound to the lattice.

Reference:

• MIT Open Course Ware 6.007

Since the electrons are decoupled from the lattice, there is no restoring force acting on the electron within the Lorentz oscillator model for the interaction of light with matter.

Hence,

\[\omega_0 = 0.\]

However, the possibility of collisions with adjacent electrons cannot be ruled out. This implies that there still acts as a damping force

Hence

\[\lambda \neq 0.\]

And due to the sea of electrons, there is a non zero plasma frequency associated with the number density of electrons.

\[\omega_p = \sqrt{ \frac{Ne^2}{m \epsilon_0} }.\]

We have shown earlier that the complex relative permitiviy is given by

\[\epsilon_r = \frac{\epsilon}{\epsilon_0} = \epsilon_{real} - \iota \epsilon_{imag},\]

where

\[\epsilon_{real} - 1 = \frac{\omega_p^2 \left( \omega_0^2 - \omega^2 \right)}{\left( \omega_0^2 - \omega^2 \right)^2 + \omega^2 \gamma^2 },\]

and

\[\epsilon_{imag} = \frac{\omega_p^2 \gamma \omega }{\left( \omega_0^2 - \omega^2 \right)^2 + \omega^2 \gamma^2 }.\]

In the case of metals,

\[\epsilon_{real} = 1 - \frac{\omega_p^2 }{ \omega^2 + \gamma^2 },\]

and

\[\epsilon_{imag} = \frac{\omega_p}{\omega} \frac{\omega_p \gamma }{ \omega^2 + \gamma^2 }.\]

Low frequency regime

For $ 0 \leq \omega \ll \gamma$,

\[\epsilon_{real} \sim - \frac{\omega_p^2 }{ \gamma^2 },\]

suggests large and negative values of real value of permitivity. and

\[\epsilon_{imag} \rightarrow \infty.\]

suggests large absorption of the radiation. The absorptions happends upto frequencies of the order of $\gamma$.

For $ \gamma \ll \omega \ll \omega_p$,

\[\epsilon_{real} \sim - \frac{\omega_p^2 }{ \omega^2 },\]

is large and negative and

\[\epsilon_{imag} \sim - \left( \frac{\omega_p }{ \omega } \right)^2 \frac{\gamma}{\omega}.\]

is still large suggesting absorption, however the absorption reducing with increasing frequency.

At $ \omega = \omega_p$,

\[\epsilon_{real} \sim 0, \quad \epsilon_{imag} \sim 0.\]

This suggests the refractive index of medium $n = \sqrt{\epsilon_r}$ tends to zero, so that velocity of light within the medium $v = c/n$ tends to infinity. Since the frequency of the incident radiation is maintained during the process, this suggests that the wavelength of the emitted radiation is $\lambda = 2\pi v/\omega$ tends to infinity! The physical interpreation of this behaviour is that the electrons across the whole medium oscillate in phase leading to a collective excitation called plasmon.

At $ \omega \gg \omega_p$,

\[\epsilon_{real} \sim 1 , \quad \epsilon_{imag} \sim 0.\]

This suggests the medium behaves as vacuum, nearly transmitting all the radiation.

Estimate of plasma frequency

Let us take the example of Zinc (Z = 30).

To estimate the plasma frequency, we need the number density Typically, the number density in metals can be esimated from the the average sphere of atom with radius $r_s$ given by

\[r_s = \left(\frac{3}{4\pi n}\right)^{1/3}\]

The typical size of the atom is of the order of Bohr radius, so it is convenient to express as

\[n = \frac{1}{\frac{4\pi a_0^3}{3}} \frac{1}{\left(\frac{r_s}{a_0}\right)^3},\] \[n = 1.612 \frac{1}{\left(\frac{r_s}{a_0}\right)^3} \mathrm{per A}^3\]

For Zn, (from Ashcroft and Mermin)

\[\frac{r_s}{a_0} = 2.30\] \[n = 13.2 \mathrm{per A}^3\]

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.style.use('../../../myMatplotlibStylesheet.mplstyle')
from scipy.constants import pi, epsilon_0, e, m_e, c, physical_constants

%matplotlib notebook

a_0 = physical_constants['Bohr radius'][0]

a_0

5.29177210903e-11

# Material properties of Zn

r_s = 2.30  # Radius of classical Sphere enclosing the atom. 

# Number density (in per m^3)

v_0 = 4.0 * pi * a_0**3 / 3
N = 1/ (v_0 * r_s**3 )

N

1.3241112753183164e+29

Hence, the plasma frequency is given by

# Estimate Plasma frequency (in rad s^-1)

w_p = np.sqrt(N * e**2 / (m_e * epsilon_0))

w_p

2.0528337200906964e+16

Therefore

$\omega_p = 2.05 \cdot 10^{16}\, \mathrm{rad s}^{-1} \quad = 3.26 \, \mathrm{PHz}$.

The visible spectrum corresponds to wavelength from [400, 700] nm. This corresponds to frequency in the range $[2.7, 4.7] \cdot 10^{15}$ rad per s.

lambdas = [400, 700] # wavelength range in nm
omegas = [c/el*2*pi*1e+9/1e+15 for el in lambdas] # omegas range in 10^15 rad s^-1

omegas

[4.709128918272133, 2.690930810441219]

Hence, in terms of the plasma frequency, the visible frequency range $[0.13\omega_p, 0.23\omega_p]$ is nearly an order of magnitude lower.

omegas_norm = [omega/w_p*1e+15 for omega in omegas]

omegas_norm

[0.22939651040338904, 0.13108372023050804]

The response in terms of these relative magnitudes can be interpreted as the visible spectrum getting mostly reflected due to Zinc, and since most of the metals have electron density of the same order, it is reasonable to infer that the plasma frequency for metals is also of the order of PHz.

# Variation of relative permitivity with normalized frequency (w/w_p)
ws = np.logspace(-3, +1, 100)
gamma = 0.1

def epsilon_real(w):
    return 1 - 1/(w**2 + gamma**2)

def epsilon_imag(w):
    return gamma/(w*(w**2 + gamma**2))

Plot the dispersion in both the high frequency and low frequency regime

Let us plot the dispersion in both the high frequency $\omega » \omega_p$ and low frequency regime $\omega » \omega_p$ simultaneously.

# Plot the dispersion in relative permitivity.

fig, ax = plt.subplots()

# kwargs common to the lines
#kwargs = dict(markevery=4)
ax.plot(ws, epsilon_real(ws), label=r'$\epsilon_\mathrm{real}(\omega)$', markevery=(0,4))
ax.plot(ws, epsilon_imag(ws), label=r'$\epsilon_\mathrm{imag}(\omega)$', markevery=(2,4))
#ax.axvline(1, color='C2')
ax.annotate(r'$\omega_p$', (1, 0), (0.5, 0.5),
            xycoords=ax.transData, textcoords=ax.transAxes,
            arrowprops=dict(arrowstyle='->'))
ax.legend()
ax.set_xscale('log')
# Set the ticks for the log scale.
major_ticks = mpl.ticker.LogLocator(base=10.0)
minor_ticks = mpl.ticker.LogLocator(base=10.0, subs=np.arange(1, 10), numticks=10)
ax.xaxis.set_major_locator(major_ticks)
ax.xaxis.set_minor_locator(minor_ticks)

ax.set_xlabel(r'$\omega/\omega_p$')
ax.set_ylabel(r'$\epsilon_r(\omega) = \epsilon_{\mathrm{real}}(\omega) + \iota \epsilon_{\mathrm{imag}}(\omega)$')

Text(0, 0.5, '$\\epsilon_r(\\omega) = \\epsilon_{\\mathrm{real}}(\\omega) + \\iota \\epsilon_{\\mathrm{imag}}(\\omega)$')

Plot the dispersion in the low and high frequency regimes as different plots.

Due to the differing scales of magnitudes of the relative permitivity for the real and imaginary, let us plot them as different plots with a common $\omega$ axis.

Let us limit the frequency range to values closer to the plasma frequency.

fig, axs = plt.subplots(2,1)
# Remove the clipping of the markers at the spines for proper rendering close the y=0.
kwargs = dict(clip_on=False)
# Also remove the facecolors of the suplot axes to avoid the markers getting masked by the face of the axes.
for ax in axs:
    ax.set_facecolor("None")

x_break = 3e-1 # break point for the right subplots.

axs[0].plot(ws[ws>=1], epsilon_real(ws[ws>=1]), 'C0o', label=r'$\epsilon_\mathrm{real}(\omega)$', markevery=(0,2), **kwargs)
axs[0].plot(ws[ws>=x_break], epsilon_imag(ws[ws>=x_break]), 'C1s', label=r'$\epsilon_\mathrm{imag}(\omega)$', markevery=(1,2), **kwargs)
axs[1].plot(ws[np.array([b1 and b2 for b1, b2 in zip(ws<1, ws>x_break)])], epsilon_real(ws[np.array([b1 and b2 for b1, b2 in zip(ws<1, ws>x_break)])]), 'C0o', label=r'$\epsilon_{real}$', **kwargs)

# Set the lims.
axs[0].set_xlim(x_break, 1e+1)
axs[1].set_xlim(x_break, 1e+1)
axs[0].set_ylim(bottom=0)
axs[1].set_ylim(top=0)

# Set the scales.
axs[0].set_xscale('log')
axs[1].set_xscale('log')

# Set the xticks  in the log scale.
major_ticks = mpl.ticker.LogLocator(base=10.0)
minor_ticks = mpl.ticker.LogLocator(base=10.0, subs=np.arange(1, 10), numticks=10)
ax.xaxis.set_major_locator(major_ticks)
ax.xaxis.set_minor_locator(minor_ticks)

# Set the labels and legend.
axs[1].set_xlabel(r'$\frac{\omega}{\omega_p}$')
kwargs_label = dict(#transform=fig.transFigure,
                     ha='center', va='center',
                     fontsize=12, fontweight=900,
                     rotation='vertical'
                    )
fig.text(0.1, 0.55, r'$\epsilon_r(\omega) = \epsilon_{\mathrm{real}}(\omega) + \iota \epsilon_{\mathrm{imag}}(\omega)$',
             **kwargs_label)
axs[0].legend()

# Highlight scale change
axs[1].annotate('scale change', (0.15, 1.0), xycoords=axs[1].transAxes, rotation='vertical', va='center')
axs[1].annotate("", (0.1, 1.25), (0.1, 0.75), xycoords=axs[1].transAxes, textcoords=axs[1].transAxes, arrowprops=dict(arrowstyle="fancy")) 


# Remove the padding between subplots
fig.subplots_adjust(hspace=0.0)
axs[0].spines["bottom"].set_visible(False)
axs[1].spines["top"].set_visible(False)



# Remove the tick labels and ticks using tick_params
axs[0].tick_params(which='both', bottom=False, labelbottom=False)
axs[1].tick_params(which='both', top=False)

# Set the scale change symbol, usually used for break in the axes range. However let us use cross symbol
# Reference: https://stackoverflow.com/questions/32185411/break-in-x-axis-of-matplotlib/32186074#32186074
d = .015 # how big to make the diagonal lines in axes coordinates
# arguments to pass plot, just so we don't keep repeating them
kwargs.update(transform=axs[0].transAxes)
axs[0].plot((0-d,0+d), (0-d,0+d), 'k-', **kwargs)
axs[0].plot((1-d,1+d),(0-d,0+d), 'k-', **kwargs)
kwargs.update(transform=axs[1].transAxes)
axs[1].plot((0-d,0+d), (1+d,1-d), 'k-', **kwargs)
axs[0].plot((1-d,1+d),(1+d,1-d), 'k-', **kwargs)

[<matplotlib.lines.Line2D at 0x1d356e92490>]

Low and high frequency regimes in different axes

fig, axs = plt.subplots(2,2)
# Remove the clipping of the markers at the edges for proper rendering close the y=0.
kwargs = dict(clip_on=False)

# Use twin axes.
tax0 = axs[0][1].twinx()
tax1 = axs[1][1].twinx()

# Also remove the facecolors of the suplot axes to avoid the markers getting masked by the face of the axes.
axs[0][1].set_facecolor("None")
axs[1][1].set_facecolor("None")
tax0.set_facecolor("None")
tax1.set_facecolor("None")

x_break = 3e-1 # break point for the right subplots.

# Plot the dispersion in high frequency regime.


tax0.plot(ws[ws>=1], epsilon_real(ws[ws>=1]), 'C0o', label=r'$\epsilon_{real}$', markevery=(0,2), **kwargs)
tax0.plot(ws[ws>=x_break], epsilon_imag(ws[ws>=x_break]), 'C1s', label=r'$\epsilon_{imag}$', markevery=(1,2), **kwargs)


tax1.plot(ws[np.array([b1 and b2 for b1, b2 in zip(ws<1, ws>x_break)])], epsilon_real(ws[np.array([b1 and b2 for b1, b2 in zip(ws<1, ws>x_break)])]), 'C0o', label=r'$\epsilon_{real}$', **kwargs)

# Set the xlims.
axs[0][1].set_xlim(x_break, 1e+1)
axs[1][1].set_xlim(x_break, 1e+1)

# Set the ylims.
tax0.set_ylim(bottom=0)
tax1.set_ylim(top=0)

# Set the scales.
tax0.set_xscale('log')
tax1.set_xscale('log')

# Set the labels.
axs[1][1].set_xlabel(r'$\frac{\omega}{\omega_p}$')

# Set the legend.
tax0.legend()

# Highlight scale change
axs[1][1].annotate('scale change', (0.15, 1.0), xycoords=axs[1][1].transAxes, rotation='vertical', va='center')
axs[1][1].annotate("", (0.1, 1.25), (0.1, 0.75), xycoords=axs[1][1].transAxes, textcoords=axs[1][1].transAxes, arrowprops=dict(arrowstyle="fancy")) 


# Remove the padding between subplots
fig.subplots_adjust(hspace=0.0)

# Remove the overlapping x-spines.
axs[0][1].spines["bottom"].set_visible(False)
axs[1][1].spines["top"].set_visible(False)

tax0.spines["bottom"].set_visible(False)
tax1.spines["top"].set_visible(False)


# Remove the tick labels using ticker

locator = mpl.ticker.NullLocator()
axs[0][1].xaxis.set_major_locator(locator)
axs[0][1].xaxis.set_minor_locator(locator)

#mpl.rcParams['xtick.top'] = False
#mpl.rcParams['xtick.minor.visible'] = False

# Set the scale change symbol, usually used for break in the axes range. However let us use cross symbol
# Reference:
# [Broken Axis tutorial](https://matplotlib.org/3.3.3/gallery/subplots_axes_and_figures/broken_axis.html#sphx-glr-gallery-subplots-axes-and-figures-broken-axis-py)
#https://stackoverflow.com/questions/32185411/break-in-x-axis-of-matplotlib/32186074#32186074
d = .03 # how big to make the diagonal lines in axes coordinates
# arguments to pass plot, just so we don't keep repeating them
kwargs.update(transform=axs[0][1].transAxes)
axs[0][1].plot((0-d,0+d), (0-d,0+d), 'k-', **kwargs)
axs[0][1].plot((1-d,1+d),(0-d,0+d), 'k-', **kwargs)
kwargs.update(transform=axs[1][1].transAxes)
axs[1][1].plot((0-d,0+d), (1+d,1-d), 'k-', **kwargs)
axs[0][1].plot((1-d,1+d),(1+d,1-d), 'k-', **kwargs)

# Plot the dispersion in the low frequency regime.
axs[0][0].set_facecolor("None")
axs[1][0].set_facecolor("None")

#axs[0][0].plot(ws[ws<=x_break], epsilon_real(ws[ws<=x_break]), 'C0o', label=r'$\epsilon_{real}$', markevery=(0,2), **kwargs)
#axs[0][0].plot(ws[ws<=x_break], epsilon_imag(ws[ws<=x_break]), 'C1s', label=r'$\epsilon_{imag}$', markevery=(1,2), **kwargs)
#axs[1][0].plot(ws[ws<=x_break], epsilon_real(ws[ws<=x_break]), 'C0o', label=r'$\epsilon_{real}$', markevery=(0,2), **kwargs)

axs[1][0].plot(ws[ws<=x_break], epsilon_real(ws[ws<=x_break]), markevery=(0,2))
axs[0][0].plot(ws[ws<=x_break], epsilon_imag(ws[ws<=x_break]), 'C1s',  markevery=(1,2))

# Set same xlims.
axs[0][0].set_xlim(left=1e-3, right=x_break)
axs[1][0].set_xlim(left=1e-3, right=x_break)

# Set ylims with matching zeros.
axs[0][0].set_ylim(bottom=0)
axs[1][0].set_ylim(top=0)

# Set log scales for x axis.
axs[0][0].set_xscale('log')
axs[1][0].set_xscale('log')


# Highlight plasma frequency
axs[0][1].axvline(1, color='C2', ls='--')
axs[1][1].axvline(1, color='C2', ls='--')
axs[1][1].annotate(r'$\omega_p$', (1, 0.0), xycoords=axs[1][1].transData, ha='center', color='C2')

# Highlight gamma
axs[0][0].axvline(gamma, color='C2', ls='--')
axs[1][0].axvline(gamma, color='C2', ls='--')
axs[1][0].annotate(r'$\gamma$', (0.1, 0.0), xycoords=axs[1][0].transData, ha='center', color='C2')


# Place break points.

kwargs.update(transform=axs[0][0].transAxes)
axs[0][0].plot((1-d,1+d), (1-d,1+d), 'k-', **kwargs)
kwargs.update(transform=axs[0][1].transAxes)
axs[0][1].plot((0-d,0+d),(1-d,1+d), 'k-', **kwargs)
kwargs.update(transform=axs[1][0].transAxes)
axs[1][0].plot((1-d,1+d), (0-d,0+d), 'k-', **kwargs)
kwargs.update(transform=axs[1][1].transAxes)
axs[1][1].plot((0-d,0+d),(0-d,0+d), 'k-', **kwargs)

# Reduce the padding between columnar subplots
fig.subplots_adjust(wspace=0.05)

# Remove the overlapping y-spines.
axs[0][1].spines["left"].set_visible(False)
axs[1][1].spines["left"].set_visible(False)
### Need to remove the spines of the twinaxes also.
tax0.spines["left"].set_visible(False)
tax1.spines["left"].set_visible(False)
axs[0][0].spines["right"].set_visible(False)
axs[1][0].spines["right"].set_visible(False)


# Remove the overlapping x-spines.
axs[0][0].spines["bottom"].set_visible(False)
axs[1][0].spines["top"].set_visible(False)

# Remove the ticklabels using ticker
locator = mpl.ticker.NullLocator()
axs[0][1].yaxis.set_major_locator(locator)
axs[0][1].yaxis.set_minor_locator(locator)
axs[1][1].yaxis.set_major_locator(locator)
axs[1][1].yaxis.set_minor_locator(locator)

axs[0][0].xaxis.set_major_locator(locator)
axs[0][0].xaxis.set_minor_locator(locator)

mpl.rcParams['ytick.right'] = False
mpl.rcParams['ytick.minor.visible'] = False
#axs[0][0].tick_params('y', width=0)
#axs[0][0].yaxis.set_major_locator(locator)
#axs[1][0].yaxis.set_minor_locator(locator)

# Revert back the rcParams
mpl.rcParams['ytick.right'] = False
mpl.rcParams['ytick.minor.visible'] = False

fig.savefig('relative_permitivity_metals.png')
fig.savefig('relative_permitivity_metals.pdf')

Using GridSpec

Let us plot the same dispersion in relative permitivity using gridspec

# Plot the dispersion using gridspec.


fig = plt.figure()

#Specify a grid and add axes to it.
gs = fig.add_gridspec(4, 2,
                      left=0.15, right=0.85,
                      width_ratios=(1, 1), height_ratios=(2, 1, 1, 2),
                      wspace=0.0, hspace=0.0
                     )

ax00 = fig.add_subplot(gs[0, 0])
ax11 = fig.add_subplot(gs[1, 1])
ax21 = fig.add_subplot(gs[2, 1])
ax30 = fig.add_subplot(gs[3, 0])

x_break = 3e-1 # break point for the right subplots.

# Plot the dispersion in high frequency regime.
kwargs = dict(clip_on=False)

ax11.plot(ws[ws>=1], epsilon_real(ws[ws>=1]), 'C0o', label=r'$\epsilon_{real}$', markevery=(0,4), **kwargs)
ax11.plot(ws[ws>=x_break], epsilon_imag(ws[ws>=x_break]), 'C1s', label=r'$\epsilon_{imag}$', markevery=(2,4), **kwargs)
ax21.plot(ws[np.array([b1 and b2 for b1, b2 in zip(ws<1, ws>x_break)])],
          epsilon_real(ws[np.array([b1 and b2 for b1, b2 in zip(ws<1, ws>x_break)])]),
          'C0o', markevery=(0,4), **kwargs)

ax30.plot(ws[ws<=x_break], epsilon_real(ws[ws<=x_break]), 'C0o', markevery=(0, 4), **kwargs)
ax00.plot(ws[ws<=x_break], epsilon_imag(ws[ws<=x_break]), 'C1s',  markevery=(2, 4), **kwargs)

# Set the legend.
fig.legend(loc='upper right')

# Remove the facecolors of the suplot axes to avoid the markers getting masked by the face of the axes.
for ax in [ax11, ax21, ax00, ax30]:
    ax.set_facecolor("None")

# Set the xlims.
ax11.set_xlim(x_break, 1e+1)
ax21.set_xlim(x_break, 1e+1)

ax00.set_xlim(left=1e-3, right=x_break)
ax30.set_xlim(left=1e-3, right=x_break)

# Set the ylims.
ax11.set_ylim(bottom=0, top=epsilon_imag(x_break))
ax21.set_ylim(top=0, bottom=epsilon_real(x_break))

ax00.set_ylim(bottom=epsilon_imag(x_break),top=np.max(epsilon_imag(ws)))
ax30.set_ylim(top=epsilon_real(x_break), bottom=np.min(epsilon_real(ws)))

# Set the scales.
for ax in [ax11, ax21, ax00, ax30]:
    ax.set_xscale('log')

# Set the xticks.
major_locator=mpl.ticker.LogLocator(base=10.0)

# NOTE** It is necessary to mention the number of ticks for the minorticks
#minor_locator=mpl.ticker.LogLocator(base=10.0)
minor_locator=mpl.ticker.LogLocator(base=10.0, subs=np.arange(1, 10), numticks=10)

for ax in [ax11, ax21, ax00, ax30]:
    ax.xaxis.set_major_locator(major_locator)
    ax.xaxis.set_minor_locator(minor_locator)

kwargs_tick_params = dict(
                            left=True,
                            bottom=True,
                            labelleft=False,
                            labelbottom=True,
                            top=False,
                            labeltop=False,
                            labelright=True
                            ) 
ax21.tick_params(which='both', **kwargs_tick_params)
kwargs_tick_params.update(bottom=False, labelbottom=False, top=True, labeltop=True)
ax11.tick_params(which='both', **kwargs_tick_params)

# Remove the spines.
ax11.spines['bottom'].set_visible(False)
ax21.spines['top'].set_visible(False)

# Set the axes labels.
kwargs_label = dict(transform=fig.transFigure,
                     ha='center', va='center',
                     fontsize=12, fontweight=900,
                    )
fig.text(0.5, 0.1, r'$\omega/\omega_p$', **kwargs_label)
fig.text(0.075, 0.55, r'$\epsilon_r, \epsilon_i$', rotation='vertical', **kwargs_label)

# Highlight important frequencies and scale changes.
ax21.annotate(r'$\omega_p$', (1, 1), (1, 0.3),
            xycoords=ax21.get_xaxis_transform(), textcoords=ax21.get_xaxis_transform(),
            ha = 'center', color='C2',
            arrowprops=dict(arrowstyle="fancy", color='C2')
           ) 
ax30.annotate(r'$\gamma$', (gamma, 0), (gamma, 0.4),
            xycoords=ax30.get_xaxis_transform(), textcoords=ax30.get_xaxis_transform(),
            ha = 'center', color='C2',
            arrowprops=dict(arrowstyle="fancy", color='C2')
           ) 
ax21.annotate('', (-0.1, 1.25), (-0.1, 0.75),
            xycoords=ax21.transAxes, textcoords=ax21.transAxes,
            ha = 'center', color='C2',
            arrowprops=dict(arrowstyle="fancy", color='blue')
           ) 
ax21.text(-0.2, 1.0, 'scale change',
            transform=ax21.transAxes,
            rotation='vertical', ha = 'center', va='center', color='blue',
           ) 

Text(-0.2, 1.0, 'scale change')

fig.savefig('dispersion_relative_permitivity_metals.png')
fig.savefig('dispersion_relative_permitivity_metals.pdf')

Modify ticklabels and spines using tick_params() and spines dictionary

The tick_params() is a convenient method of Axes to control the visibility of ticklabels. Similarly, spines, a dictionary, is an attribute of Axes to manage the spines.

def lorentzian(ws):
    del_w = 1 - ws
    return gamma/(2*np.pi)*1/(del_w**2 + gamma**2/4.0)

# Remove the bottom and left ticklabels and spines.


fig, ax = plt.subplots()
ax.plot(ws, lorentzian(ws), label='lorentzian')
ax.set_xscale('log')
# Set the logscale minor ticks.
major_ticks = mpl.ticker.LogLocator(base=10.0)
minor_ticks = mpl.ticker.LogLocator(base=10.0, subs=np.arange(1,10), numticks=10)
ax.xaxis.set_major_locator(major_ticks)
ax.xaxis.set_minor_locator(minor_ticks)

kwargs_tick_params = dict(
                            left=False,
                            bottom=False,
                            labelleft=False,
                            labelbottom=False,
                            labeltop=True,
                            labelright=True
                            ) 
ax.tick_params(which='both', **kwargs_tick_params)

# Remove the spines.

ax.spines['left'].set_visible(False)
ax.spines['bottom'].set_visible(False)

ax.legend()
x_label = ax.set_xlabel(r'$\omega/\omega_0$')
y_label = ax.set_ylabel('I')
y_label.set_position((1.0, 0.5))
#y_label.set_transform(ax.transAxes)

# Annotate the resonance peak. 
# Use blended transform.
# Ref : https://matplotlib.org/3.1.0/tutorials/advanced/transforms_tutorial.html

ax.annotate(r'$\omega_0$', (1, 0.1), (1, -0.1),
            xycoords=ax.get_xaxis_transform(), textcoords=ax.get_xaxis_transform(),
            ha = 'center',
            arrowprops=dict(arrowstyle="fancy")
           ) 
text_reso = ax.text(0.5, 0.5, 'resonance', ha='center')
text_reso.set_position((0.0, 0.0))
text_reso.set_transform(ax.transAxes)

fig.savefig('dispersion_relative_permitivity_metals.png')

Plot an x2-y2 plot using twinx() and twiny()

Earlier, we used the tick_params() to plot x2-y2. Now let us higher level methods, twinx() and twiny() of Axes to do the same. These methods create new instances of axes that are drawn on top of the original Axes instance.

fig, ax = plt.subplots()
tax = ax.twinx()
taxy = tax.twiny()

# Plot on x2 y2 axes.

taxy.plot(np.linspace(0, np.pi)/np.pi, np.sin(np.linspace(0, np.pi)), label='sin(x)')

# Set the legend and labels.
taxy.legend()
taxy.set_xlabel('x')
# NOTE: To set y_label, we need to set that of first twin.
tax.set_ylabel('y')

# Format the xticklabels in units of pi.
xtick_labels = mpl.ticker.FuncFormatter(lambda x, pos: "%.2f" % x + r"$\pi$" )
taxy.xaxis.set_major_formatter(xtick_labels)

# Remove the ticks of the original axes.

null_locator = mpl.ticker.NullLocator()
ax.yaxis.set_major_locator(null_locator)
tax.xaxis.set_major_locator(null_locator)

#ax.plot(np.linspace(0, 2*np.pi), 10*np.cos(np.linspace(0, 2*np.pi)), 'C1')

Dispersion in refractive index

The refractive index is directly connected to the relative permitivity. This connection is evident when we consider the propation of electromagnetic wave. The celebrated governing equation for the EM wave propagation is given by

\[\frac{\partial^2 }{\partial x^2} \left\{ \begin{array}{l} E \\ H \end{array} \right. = \frac{1}{c^2}\frac{\partial^2 }{\partial t^2} \left\{ \begin{array}{l} E \\ H \end{array} \right.,\]

so that for a harmonic wave

\[\left. \begin{array}{r} E \\ H \end{array} \right\} (\mathbf{r}, t) = \left. \begin{array}{l} E_0 \\ H_0 \end{array} \right\} e^{\iota \overline{ \mathbf{k}\cdot\mathbf{r} - \omega t} }\]

as a trial solution, we have a travelling wave with a velocity given by its frequency and wavevector.

\[c = \frac{\omega}{k}.\]

Further it can be shown that wavevector in a medium is modified by

\[k = k_0\sqrt{\epsilon_r}\]

that is the wavevector stretches'' in aheavier’’ medium so that the velocity in the medium is

\[c' = c / {\epsilon_r}.\]

Since, the refractive index is the ratio of velocity of medium to that in vacuum,

we have

\[n = \sqrt{\epsilon_r}\]

Since relative permitivity is a complex quantity, so is the refractive index.

# Define the complex refractive index from the relative permitivity.

def n_real(w):
    e_r = epsilon_real(w)
    e_i = epsilon_imag(w)
    two_n_square = e_r + np.sqrt(e_r**2 + e_i**2)
    return np.sqrt(two_n_square/2)

def n_imag(w):
    e_r = epsilon_real(w)
    e_i = epsilon_imag(w)
    two_n_square = -1*e_r + np.sqrt(e_r**2 + e_i**2)
    return np.sqrt(two_n_square/2)

# Variation of relative permitivity with normalized frequency (w/w_p)
ws = np.logspace(-4, +1, 100)
gamma = 0.1

fig1, ax1 = plt.subplots()

ax1.plot(ws, n_real(ws), 'C0o', label=r'$n_\mathrm{real}(\omega)$')
ax1.plot(ws, n_imag(ws), 'C1s', label=r'$n_\mathrm{imag}(\omega)$')
ax1.set_xlim(1e-4, 1e+1)
ax1.set_xscale('log')
ax1.set_yscale('log')
# Set the logscale minor ticks.
major_ticks = mpl.ticker.LogLocator(base=10.0)
minor_ticks = mpl.ticker.LogLocator(base=10.0, subs=np.arange(1,10), numticks=10)
ax1.xaxis.set_major_locator(major_ticks)
ax1.xaxis.set_minor_locator(minor_ticks)

ax1.set_xlabel(r'$\frac{\omega}{\omega_p}$')
ax1.set_ylabel(r'$n(\omega) = n_{\mathrm{real}}(\omega) + \iota n_{\mathrm{imag}}(\omega)$')

ax1.legend()

# Highlight the important frequencies.

ax1.axvline(1, color='C2')
ax1.axvline(gamma, color='C2')

kwargs_vspan = dict(color='k', ec=None)
ax1.axvspan(1e-4, gamma/2, alpha=0.7, **kwargs_vspan)
ax1.axvspan(gamma/2, 1, alpha=0.3, **kwargs_vspan)
ax1.axvspan(1, 1e+1, alpha=0.1, **kwargs_vspan)

# Annotate the regions of transmission, absorption and reflection.
# Use blended transform
ax1.text(1e-3, 0.3, 'A', transform=ax1.get_xaxis_transform())
ax1.text(2e-1, 0.3, 'R', transform=ax1.get_xaxis_transform())
ax1.text(2e-0, 0.3, 'T', transform=ax1.get_xaxis_transform())

# Annotate the visible spectrum.

ax1_cmap = ax1.inset_axes([0.13, 0.0, 0.1, 1.0], transform=ax1.get_xaxis_transform(), zorder=1)


#ax1_cmap.spines['left'].set_visible(False)
color_data = mpl.cm.ScalarMappable(cmap='rainbow_r')
cbar = fig1.colorbar(color_data, cax=ax1_cmap, orientation='horizontal', alpha=0.1)

#ax1_cmap.set_xscale('log')
ax1_cmap.tick_params(which='both', bottom=False, labelbottom=False)

for spine in ['top', 'bottom', 'left', 'right']:
    ax1_cmap.spines[spine].set_visible(False)


#ax1_rainbow = fig1.add_axes([0.01, 0.1, 0.04, 0.05], transform=ax1.transAxes)
                                #, transform=ax1.get_xaxis_transform())
#ax1_rainbow.tick_params(which='both', bottom=False, labelbottom=False)

fig1.savefig('dispersion_refractive_index_metals.png')
fig1.savefig('dispersion_refractive_index_metals.pdf')

Draw inset axes

Let us draw inset axes using the higher level method inset_axes() of Axes instance and add_axes() method of Figure instance.

fig, ax = plt.subplots()
ax1 = fig.add_axes([0.5, 0.5, 0.2, 0.2])
#ax1 = fig.add_axes([0.5, 0.5, 0.2, 0.2], transform=fig.transFigure)
ax1.spines['left'].set_visible(False)

ax2 = ax.inset_axes([0.5, 0.5, 0.2, 0.2])
ax2.spines['left'].set_visible(False)

data = np.random.randn(8, 8)
cax = ax.imshow(data, alpha=0.5)
cbar = fig.colorbar(cax, cax=ax2)

Estimate of damping factor

Let us estimate the order of magnitude of $\gamma$. For this we need to look at the sources for the reaction force.

• Natural source is the emission of radiation itself that generates a reaction onto the accelerating charge. The width in the frequency spectrum associated with the nearly monochromatic radiation is called natural line width
• The harmonic nature of the wave is deviated from the collisions among charges. This is called collisional broadening.

For a particle emitting radiation with frequency $\omega$, the self force $\mathbf{F}_s$ or the reaction force can be shown to be proportional to the velocity $\mathbf{v}$ Heitler1935, so that

\[\mathbf{F}_s = - \frac{2}{3} \frac{e^2 \omega^2}{c^3} \mathbf{v}\]

Since we have modelled the damping force as

\[\mathbf{F}_s = - m \gamma \mathbf{v},\]

we have

\[\gamma = \frac{2}{3} \frac{e^2 \omega^2}{m c^3}.\]

In terms of the classical electron radius

\[r_e = \frac{e^2}{m c^2} = 2.82 \, \mathrm{fm},\]

we have

\[\gamma = \frac{2}{3} \frac{e^2 \omega^2}{c} r_e = \frac{2}{3} \left( \frac{\omega r_e}{c} \right) \omega = \frac{2}{3} \omega \left( r_e k \right) .\]

For frequencies close to the visible spectrum, $\omega \sim 3 \cdot 10^{15} \,\mathrm{rads}^{-1}$ or $k \sim 10^{7}\,\mathrm{m}^{-1}$, so that

\[r_e k \ll 1,\]

so that

\[\gamma \ll \omega.\]

Collisional broadening

In the previous section, we have seen how the emission of radiation from a harmonic oscillator is no longer monochromatic but instead the decreasing amplitude of the field strength with time is reflected as a spread of intensity around the frequency $\omega_0$ of the oscillator.

Now let us consider the situation where not only the amplitude but also the phase of the emission changes due to changes in the phase of the oscillator. The change in phase can be brought about by collisions with the neighbouring oscillators.

Now the nature of collisions is stochastic. As a result, the effect of these collisions needs to be accounted statistically. Let categorize the events as per the time $\tau$ it takes for the next event. So the number of events with interval between $\tau$ and $\tau + d\,\tau$ be $w(\tau)d\, \tau$, where $w(\tau)$ is the number of collisions happening per unit time that correspond to the elapsed time $\tau$ between collisions. $w(\tau)$ is also called the collision frequency.

Lorentz assumed a particular form for the collision frequency as

\[w(\tau) = \frac{1}{2} \Gamma e^{-\Gamma \tau/2}\]

This form qualitatively captures the behaviour that events corresponding to larger time interval between collisions are exponentially rarer that that of smaller time intervals.

Due to such temporal distribution of the collisions, we can talk of the average time between collisions as

\[\overline{\tau} = \frac{\int_0^\infty t w(\tau) \, d \tau}{\int_0^\infty w(\tau) \, d \tau}\]

Using the Lorentzian distribution for the collision frequency, we have

\[\overline{\tau} = \frac{\int_0^\infty t e^{-\Gamma \tau/2} \, d \tau}{\int_0^\infty e^{-\Gamma \tau/2} \, d \tau}\] \[\therefore \overline{\tau} = \frac{2}{\Gamma}.\]

The relaxation time is twice the reciprocal time $1/\Gamma$.

Let us further assume that the effect of the collision is to shift the phase of the oscillator, and thereby the emitted wave by a random value.

So the temporal extent of the wave is not truncated by the dampind force due to radiation reaction, but instead the emitted wave has a finite temporal width.

So for a wave with frequency $\omega$ emanating for a time interval $\tau$ from an oscillator of natural frequency $\omega_0$ has a frequency distribution given by

\[E(t) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} E(\omega) e^{i\omega t} \, d \omega,\]

where

\begin{align} E(\omega) &= \frac{1}{2\pi} \int_{-\infty}^{+\infty} E(t) e^{-i\omega t} \, d t,
& = \frac{1}{2\pi} \int_{0}^{\tau} E_0 e^{+i\left(\omega_0 - \omega \right) t} \, d t, \end{align}

\[\therefore E(\omega) = \frac{E_0}{2\pi} \frac{e^{+i\left(\omega_0 - \omega \right) \tau} - 1}{i (\omega_0 - \omega)}\]

Now the contribution to the frequency $\omega_0$ from the wave of time period $\tau$ needs to be weighed by the collision frequency to get the average fourier amplitude of the mode $\omega$.

\[\overline{| E(\omega)|^2} = \int_0^\infty \omega(\tau) | E(\omega)|^2 \, d\tau\] \[\therefore \overline{| E(\omega)|^2} = \frac{E_0}{2\pi^2}\frac{1}{\left(\omega_0 - \omega \right)^2 + \Gamma^2/4}.\]

So the line shape is again a Lorentzian. However, the line breadth is given by $\Gamma$ where

\[\Gamma = \frac{2}{\overline{\tau}},\]

so the line breadth due to the coliisions is twice the reciprocal of average time between collisions.

Equivalently,

\[\frac{1}{\Gamma} = \frac{\overline{\tau}}{2}\]

is the lifetime due to collisions similar to $1/\gamma$ is the radiative lifetime.

Estimate of collisional broadening

One can estimate the mean time between collisions using Drude’s model of electrical transport in metals that describes the experimentally observed Ohm’s law.

Ohm’s law states that the current density is proportional to the electric field

\[\mathbf{j} = \sigma \mathbf{E}\]

If $\tau$ is the mean time between collisions, also called the relaxation time, the drift velocity $\mathbf{v}_d $ picked up by the electron in the presence of an electric field $\mathbf{E}$ is

\[\mathbf{v}_d = \frac{e\tau}{m}\mathbf{E}.\]

Let us relate the current density to the number density. Since current density is the flow of charge per unit area per unit time, in a time interval $\tau$, the number of electrons that cross a section of area A is

\[N = n v_d \tau A\]

where $n$ is the number density, so that

\[\mathbf{j} = \frac{Ne}{A \tau} = n e v_d = \frac{n e^2 \tau}{m}\mathbf{E}\]

Replacing this expression in the Ohm’s law, we have

\[\sigma = \frac{n e^2 \tau}{m}.\]

As a result, we can estimate $\tau$ from the experimentally observed $\sigma, n$.

From Ashcroft and Mermin,

\[n^\mathrm{Zn} = 13.2\cdot 10^{22} \, \mathrm{cc}^{-3},\] \[\rho^\mathrm{Zn, 273K} = 5.5 \, \mu\Omega\mathrm{cm}.\]

This gives relaxation time $\tau = 4.8 \cdot 10^{-15}\, \mathrm{s}$.

rho = 5.5e-8 ## resistivity in ohm m.
n_e = 1.32e29 ## number density in per m^3

def tau(rho, n_e):
    return m_e/n_e/e**2/rho

relaxation_time = tau(rho, n_e)

relaxation_time

4.888004391932886e-15

Using this relaxation time, the line breadth due to collisional broadening is $\Gamma = 2/\overline{\tau} = 4.1 \cdot 10^{14}\, \mathrm{rad s}^{-1} = 0.2 \omega_p$. We can see that for metals, line breadth is an order of magnitude lesser than the plasma frequency.

line_breadth = 2/relaxation_time

line_breadth/1e+14

4.091649351421983

line_breadth/w_p

0.019931713471859815

References

Heitler1935 – Section 4 of Chapter 1 of Heitler’s “Quantum theory of radiation” Oxford University Press 1935.