mattisbardeen.py

"""Calculate the electrical properties of superconductors using Mattis-Bardeen
theory.

Warning: The integrals are pretty inefficient... be careful and check your
results!

References:

    [1] A. R. Kerr, “Surface Impedance of Superconductors and Normal
    Conductors in EM Simulators,” Green Bank, West Virginia, 1996.

    [2] V. Y. Belitsky and E. L. Kollberg, “Superconductor–insulator–
    superconductor tunnel strip line: Features and applications,” J.
    Appl. Phys., vol. 80, no. 8, pp. 4741–4748, Oct. 1996.

"""

from copy import deepcopy

import numpy as np
import scipy.constants as sc
# import numba


# Constants (all energies in eV)
kb = sc.k / sc.e
h = sc.h / sc.e
hbar = sc.h / (2 * sc.pi) / sc.e


# Default parameters (for niobium, from ref. [2])
PARAM = dict(Tc=9.0,              # critical temperature [K]
             Vgap0=2.862e-3,      # gap voltage at T=0K [V]
             sigma_n=1.565e7,     # normal-state conductivity [1/ohm*m]
             lambda0=85*sc.nano,  # London penetration depth at T=0K [m]
             )


# Conductance ----------------------------------------------------------------

def conductance(f, T, Vgap, **kwargs):
    """Calculate conductance from Mattis-Bardeen theory.

    Args:
        f (ndarray): frequency, in units [Hz]
        T (float): temperature, in units [K]

    Keyword arguments:
        Tc (float): critical temperature, in units [K]
        Vgap0 (float): gap voltage at T=0K, in units [V]
        sigma_n (float): normal-state conductance, in units [1/ohm*m]
        lambda0 (float): London penetration depth at T=0K, in units [m]

    Returns:
        conductance

    """

    # Import keyword arguments
    tmp = deepcopy(PARAM)
    tmp.update(kwargs)
    kw = tmp

    if isinstance(f, float) or isinstance(f, int):
        sigma = _conductance(f, T, Vgap, **kw)
        return sigma
    elif isinstance(f, np.ndarray):
        sigma = np.empty_like(f, dtype=complex)
        for i in range(len(f)):
            sigma[i] = _conductance(f[i], T, Vgap, **kw)
        return sigma
    else:
        raise ValueError


def _conductance(f, T, Vgap, **kwargs):
    """Calculate conductance from Mattis-Bardeen theory at a single frequency
    point.

    Args:
        f (float): frequency, in units [Hz]
        T (float): temperature, in units [K]

    Keyword arguments:
        Tc (float): critical temperature, in units [K]
        Vgap0 (float): gap voltage at T=0K, in units [V]
        sigma_n (float): normal-state conductance, in units [1/ohm*m]
        lambda0 (float): London penetration depth at T=0K, in units [m]

    Returns:
        conductance

    """

    # Unpack keyword arguments
    # tc = kwargs['Tc']
    sigma_n = kwargs['sigma_n']

    # Other parameters
    delta = Vgap / 2   # binding energy
    fgap = Vgap / h    # gap frequency
    w = 2 * sc.pi * f  # angular frequency of signal

    # It's a superconductor -> infinite conductivity at f=0Hz
    if f == 0:
        return sigma_n * (1e10 - 1j * 1e10)

    def _g(energy, wt, deltat):
        """The g(E) function in the Mattis Bardeen paper."""

        e0 = np.abs(energy ** 2 + deltat ** 2 + hbar * wt * energy)
        e1 = (energy ** 2 - deltat ** 2) ** 0.5
        e2 = ((energy + hbar * wt) ** 2 - deltat ** 2) ** 0.5

        return e0 / (e1 * e2)

    # Thermal contribution ---------------------------------------------------
    # Contribution from thermally excited quasi-particles

    # Play with these values to optimize the integration...
    int1_de = 1e-8
    int1_npts = 10000
    int1_max = 10.

    start = np.log10((delta+int1_de)/delta)
    stop = np.log10(int1_max)

    # Calculate intergral
    e = np.logspace(start, stop, int1_npts) * delta
    # e = np.linspace(int1_de+delta, 10*delta, int1_npts)
    func = (fermi(e, T) - fermi(e + hbar * w, T)) * _g(e, w, delta)
    thermal = 2 / hbar / w * np.trapz(func, e)

    # plt.figure()
    # plt.plot(e / delta, func)
    # plt.axvline(1)
    # plt.show()

    # Radiation contribution -------------------------------------------------
    # Contribution from quasi-particles generated by radiation f>fgap

    # TODO: switch to log spacing

    # Play with these values to optimize the integration...
    int2_de = 1e-8

    if f > fgap:
        e = np.arange(delta-hbar*w+int2_de, -delta, int2_de)
        func = (1 - 2 * fermi(e + hbar * w, T)) * _g(e, w, delta)
        radiation = 1 / hbar / w * np.trapz(func, e)
    else:
        radiation = 0

    # plt.figure()
    # plt.plot(e / delta, func)
    # plt.show()

    # Cooper pair contribution -----------------------------------------------

    # TODO: switch to log spacing

    # Play with these values to optimize the integration...
    int3_de = 1e-8

    if f > fgap:
        e = np.arange(-delta+int3_de, delta, int3_de)
    else:
        e = np.arange(delta-hbar*w+int3_de, delta, int3_de)

    line1 = 1 - 2*fermi(e + hbar * w, T)
    line2 = e**2 + delta**2 + hbar*w*e
    line3 = (delta**2 - e**2)**0.5 * ((e + hbar*w)**2 - delta**2)**0.5
    func = line1 * line2 / line3
    cooper = 1 / hbar / w * np.trapz(func, np.real(e))

    # plt.figure()
    # plt.plot(e/delta, func)
    # plt.show()

    # Conductance ------------------------------------------------------------

    sigma = sigma_n * (thermal + radiation - 1j*cooper)

    return sigma


# Surface Impedance ----------------------------------------------------------

def surface_impedance(f, d, T, Vgap, method='MB', **kw):
    """Calculate the surface impedance of a superconductor using Mattis-
    Bardeen theory.

    Args:
        f (ndarray): frequency, in units [Hz]
        d (float): superconductor thickness, in units [m]
        T (float): temperature, in units [K]
        Vgap (float): gap voltage, in units [V]
        method (str): 'MB' for Mattis-Bardeen, or 'simple' for the simplified
            surface impedance from Kerr (1996).

    Keyword arguments:
        Tc (float): critical temperature, in units [K]
        Vgap0 (float): gap voltage at T=0K, in units [V]
        sigma_n (float): normal-state conductance, in units [1/ohm*m]
        lambda0 (float): London penetration depth at T=0K, in units [m]

    Returns:
        surface impedance

    """

    # Import keyword arguments (use defaults if not specified)
    tmp = deepcopy(PARAM)
    tmp.update(kw)
    kw = tmp

    # Using Mattis-Bardeen theory...
    if method.lower() == 'mb':
        if isinstance(f, float) or isinstance(f, int):
            return _surface_impedance(f, d, T, Vgap, **kw)
        elif isinstance(f, np.ndarray):
            zs = np.empty_like(f, dtype=complex)
            for i in range(len(f)):
                zs[i] = _surface_impedance(f[i], d, T, Vgap, **kw)
            return zs

    # Using "simple" method from [1]...
    elif method == 'simple':
        dl = _lambda(T, **kw)
        return _surface_impedance_simple(f, d, dl)

    # Method not recognized...
    else:
        raise ValueError("Method not recognized")


def _surface_impedance(f, d, T, Vgap, **kw):
    """Calculate the surface impedance of a superconductor at a single
    frequency point.

    Args:
        f (float): frequency, in units [Hz]
        d (float): superconductor thickness, in units [m]
        T (float): temperature, in units [K]
        Vgap (float): gap voltage, in units [V]
        method (str): 'MB' for Mattis-Bardeen, or 'simple' for the simplified
            surface impedance from Kerr (1996).

    Keyword arguments:
        Tc (float): critical temperature, in units [K]
        Vgap0 (float): gap voltage at T=0K, in units [V]
        sigma_n (float): normal-state conductance, in units [1/ohm*m]
        lambda0 (float): London penetration depth at T=0K, in units [m]

    Returns:
        surface impedance

    """

    # Angular frequency
    w = 2 * sc.pi * f

    # Conductance (using Mattis-Bardeen theory)
    sigma = _conductance(f, T, Vgap, **kw)

    # Surface impedance (Eqn. 5 in [2])
    zs = (np.sqrt(1j * w * sc.mu_0 / sigma) / np.tanh(np.sqrt(1j * w * sc.mu_0 * sigma) * d))

    return zs


def _surface_impedance_simple(f, d, dl):
    """Calculate the surface impedance of a superconductor using the simple
    model presented in [1]. This techniques assumes that the frequency much
    lower than the gap frequency.

    Args:
        f (ndarray): frequency, in units [Hz]
        d (float): superconductor thickness, in units [m]
        dl (float): London penetration depth

    Returns:
        surfance impedance

    """

    # Angular frequency
    w = 2 * sc.pi * f

    # Surface impedance (see Appendix A5 in [1])
    zs = 1j * w * sc.mu_0 * dl / np.tanh(d / dl)

    return zs


# Penetration depths ---------------------------------------------------------

def efield_penetration(zs, f):
    """Calculate the electric field penetration depth.

    Args:
        zs (ndarray): surface impedance
        f (ndarray): frequency, in units [Hz]

    Returns:
        E-field penetration depth

    """

    # Angular frequency
    w = 2 * sc.pi * f

    return zs.imag / (w * sc.mu_0)


def hfield_penetration(zs, f):
    """Calculate the magnetic field penetration depth.

    Args:
        zs (ndarray): surface impedance
        f (ndarray): frequency, in units [Hz]

    Returns:
        H-field penetration depth

    """

    # Angular frequency
    w = 2 * sc.pi * f

    return zs.real / (w * sc.mu_0)


# Lambda ---------------------------------------------------------------------

def _lambda0(Vgap0=None, sigma_n=None, **kw):

    delta = Vgap0 / 2

    return (hbar / sc.pi / sc.mu_0 / sigma_n / delta) ** 0.5


def _lambda(T, **kw):

    t = T / kw['Tc']

    return _lambda0(**kw) / np.sqrt(1 - t ** 4)


# Fermi function -------------------------------------------------------------

def fermi(energy, T):
    """Fermi function.

    Args:
        energy (ndarray): energy in units [eV]
        T (float): temperature in units [K]

    Returns:
        fermi function

    """

    if T != 0.:
        beta = 1 / (kb * T)
        ff = 0.5 * (1 - np.tanh(0.5 * beta * energy))
        return ff
    else:
        ff = np.zeros_like(energy, dtype=float)
        ff[energy == 0] = 0.5
        ff[energy > 0] = 1.
        return ff