"""Methods for the calculation of optical properties.
Provided functions:
* getPermittivity_1D(2D, 3D): for the calculation of the permittivity from the complex \
conductivity for 1-, 2- and 3D systems.
* getReflectivity: for the calculation of the reflectivity from the permittivity via the \
permittivity, calculated with the first convenience function.
* getAbsorptionCoefficient: for the calculation of the absorption coefficient from the \
permittivity, via the extinction coefficient, i.e. imaginary part of the refractive index.
Copyright © 2019-2023 HQS Quantum Simulations GmbH. All Rights Reserved.
"""
import numpy as np
from .qolossal_tools import (
_EPSILON_0,
_HBAR_in_EV_S,
_C_in_M_over_S,
dir2idx,
)
[docs]
def getPermittivity_3D(
conductivity: np.ndarray,
w: np.ndarray,
tensor_element: str,
) -> np.ndarray:
r"""Compute the permittivity for 3D systems.
The permittivity is calculated directly from the complex optical conductivity. By manipulating
the Maxwell's equations, trying to express the quantities as a function of external fields
only, one can find
.. math::
\varepsilon_{\alpha \beta}(\omega) = \delta_{\alpha \beta}
+ \frac{i \sigma(\omega)}{\varepsilon_0 \omega}
Reference: https://en.wikipedia.org/wiki/Ewald%E2%80%93Oseen_extinction_theorem
Args:
conductivity (np.ndarray): Complex conductivity in units of e^2/(h Å).
w (np.ndarray): Frequencies at which the conductivity has been evaluated.
tensor_element (str): Tensor element of the provided conductivity.
Raises:
ValueError: if the passed frequencies contain values too close to zero.
Returns:
np.ndarray: specified tensor element of the permittivity at the provided frequencies.
"""
tolerance = 1e-12
if np.min(np.abs(w)) < tolerance:
message = "Error: the passed frequencies contain elements too close to zero. "
message += f"Tolerance {tolerance: .2e}."
raise ValueError(message)
# use the helper tool since it checks that the tensor_element is in the correct format
directionL, directionR = dir2idx(tensor_element)
# epsilon_ab = delta_ab + i sigma / (eps_0 w)
# the factor of 2 pi dividing the conductivity is because we chose hbar=1 while the units of
# conductivity contain a 1/h. This factor converts that h into an hbar.
return (directionL == directionR) + 1.0j * conductivity / (2 * np.pi) / (_EPSILON_0 * w)
[docs]
def getPermittivity_2D(
conductivity: np.ndarray, w: np.ndarray, tensor_element: str, thickness: float
) -> np.ndarray:
r"""Compute the permittivity for 2D systems.
The permittivity is calculated directly from the complex optical conductivity. By manipulating
the Maxwell's equations, trying to express the quantities as a function of external fields
only, one can find
.. math::
\varepsilon_{\alpha \beta}(\omega) = \delta_{\alpha \beta}
+ \frac{i \sigma(\omega)}{\varepsilon_0 \omega}
Reference: https://en.wikipedia.org/wiki/Ewald%E2%80%93Oseen_extinction_theorem
Args:
conductivity (np.ndarray): 2D complex conductivity in units of e^2/h.
w (np.ndarray): Frequencies at which the conductivity has been evaluated.
tensor_element (str): Tensor element of the provided conductivity.
thickness (float): Thickness of the 2D sheet.
Returns:
np.ndarray: specified tensor element of the permittivity at the provided frequencies.
"""
# simply normalize the conductivity by the sheet thickness and then use the 3D permittivity
# This normalization effectively makes the conductivity a "3D" quantity.
conductivity_norm = conductivity / thickness
return getPermittivity_3D(conductivity_norm, w, tensor_element)
[docs]
def getPermittivity_1D(
conductivity: np.ndarray, w: np.ndarray, tensor_element: str, cross_section: float
) -> np.ndarray:
r"""Compute the permittivity for 1D systems.
The permittivity is calculated directly from the complex optical conductivity. By manipulating
the Maxwell's equations, trying to express the quantities as a function of external fields
only, one can find
.. math::
\varepsilon_{\alpha \beta}(\omega) = \delta_{\alpha \beta}
+ \frac{i \sigma(\omega)}{\varepsilon_0 \omega}
Reference: https://en.wikipedia.org/wiki/Ewald%E2%80%93Oseen_extinction_theorem
Args:
conductivity (np.ndarray): 1D complex conductivity in units of (e^2/h) Å.
w (np.ndarray): Frequencies at which the conductivity has been evaluated.
tensor_element (str): Tensor element of the provided conductivity.
cross_section (float): Cross section of the 1D wire.
Returns:
np.ndarray: specified tensor element of the permittivity at the provided frequencies.
"""
# simply normalize the conductivity by the cross section and then use the 3D permittivity
# This normalization effectively makes the conductivity a "3D" quantity.
conductivity_norm = conductivity / cross_section
return getPermittivity_3D(conductivity_norm, w, tensor_element)
[docs]
def getReflectivity(
permittivity: np.ndarray,
) -> np.ndarray:
r"""Compute the normal incidence reflectivity.
The reflectivity is computed as
.. math::
\frac{(1 - n)^2 + k^2}{(1 + n)^2 + k^2}
with :math:`N(\omega) = n(\omega) + i k(\omega) = \sqrt{\varepsilon(\omega)}` the complex
refractive index.
Reference: https://link.springer.com/article/10.1007/s11664-998-0094-3
Args:
permittivity (np.ndarray): permittivity of the system.
Returns:
np.ndarray: normal incidence reflectivity.
"""
# calculate complex index of refraction
N = np.sqrt(permittivity)
# simply according to literature.
return ((1 - N.real) ** 2 + N.imag**2) / ((1 + N.real) ** 2 + N.imag**2)
[docs]
def getAbsorptionCoefficient(
permittivity: np.ndarray,
w: np.ndarray,
) -> np.ndarray:
r"""Compute the absorption coefficient.
Defined as
.. math::
\alpha = \frac{2 \omega k(\omega)}{c}
with :math:`k(\omega) = {\rm Im}[\sqrt{\epsilon(\omega)}]` the extinction coefficient.
Reference: https://en.wikipedia.org/wiki/Ewald%E2%80%93Oseen_extinction_theorem
NOTE: the output is in units of m^-1.
Args:
permittivity (np.ndarray): permittivity of the system.
w (np.ndarray): Frequencies at which the permittivity has been evaluated.
Returns:
np.ndarray: absorption coefficient.
"""
# calculate complex index of refraction
N = np.sqrt(permittivity)
# simply according to literature.
# w / _HBAR is the frequency in s^-1. Dividing by _C in m/s we obtain the result in m^-1
return 2 * (w / _HBAR_in_EV_S) * N.imag / _C_in_M_over_S