#!/usr/bin env python3 from formulary import * from scipy.constants import Boltzmann as kB, hbar, electron_volt # DEVICES # metal sc: schottky barrier def schottky_barrier(): fig, axs = plt.subplots(1, 3, figsize=(width_formula, height_default*0.6)) WD = 5 q = 1 ND = 1 eps = 11 dx = WD/10 xs = np.linspace(-WD/5, 6/5*WD, 300) rho_S = q*ND Q = rho_S * WD rho_M = -Q/dx @np.vectorize def rho_approx(x): if x < -dx: return 0.0 if x < 0: return rho_M if x < WD: return rho_S return 0.0 rhos_approx = rho_approx(xs) @np.vectorize def E(x): if x < -dx: return 0.0 if x < 0: return -rho_M/eps * (-dx-x) if x < WD: return -rho_S/eps * (WD-x) return 0.0 Es = E(xs) @np.vectorize def phi(x): # if x < -dx: return 0.0 # if x < 0: return -rho_M/(2*eps) * (dx**2-(dx-x)**2) if x < 0: return 0.0 if x < WD: return rho_S/(2*eps) * (WD**2-(WD-x)**2) return q*ND/(2*eps) * WD**2 phis = phi(xs) for ax in axs: ax.set_xlabel("$x$") ax.set_xticks([0,WD]) ax.set_xticklabels(["0", r"$W_\text{D}$"]) ax.set_yticks([0]) ax.set_yticklabels(["0"]) axs[0].plot(xs, rhos_approx) axs[0].set_ylabel(r"$\rho(x)$") axs[1].plot(xs, Es) axs[1].set_ylabel(r"$\mathcal{E}(x)$") axs[2].plot(xs, phis) axs[2].set_ylabel(r"$\phi(x)$") return fig # Charge carrier density def fn_i(T, NV, NC, Egap): return np.sqrt(NV*NC) * np.exp(-Egap*electron_volt/(2*kB*T)) def fn_freeze(T, NV, NC, Egap): return np.sqrt(NV*NC) * np.exp(-0.07*Egap*electron_volt/(2*kB*T)) def charge_carrier_density(): fig, ax = plt.subplots(1, 1, figsize=size_formula_normal_default) # Ts = np.power(10, np.linspace(-4, 0, 100)) N = 500 ts = np.linspace(0.01, 20, N) Ts = 1000/ts # Ts = np.linspace(2000, 50, N) # ts = 1000/Ts # ts = 1/Ts NV = 1e23 NC = 1e21 Egap = 2.4 n_is = fn_i(Ts, NV, NC, Egap) n_sat = np.empty(Ts.shape) n_sat.fill(1e15) idx1 = np.argmin(np.abs(n_is-n_sat)) print(idx1, N) # TODO make quadratic simple n_total = blend_arrays(n_is, n_sat, idx1, idx1+10, "linear") n_freeze = fn_freeze(Ts, NV, NC, Egap) idx2 = np.argmin(np.abs(n_freeze-n_sat)) print(idx1, idx2, N) n_total = blend_arrays(n_total, n_freeze, idx2-10, idx2+10, "quadratic_simple") # ax.plot(ts, n_is, label="Intrinsic") # ax.plot(ts, n_sat, label="Saturation") # ax.plot(ts, n_freeze, label="Freeze-out") # ax.plot(ts, n_total, label="Total", linestyle="dashed", color="black") n_total2 = n_is.copy() n_total2[idx1:idx2] = n_sat[idx1:idx2] n_total2[idx2:] = n_freeze[idx2:] ax.plot(ts, n_is, label="Intrinsic", linestyle="dashed") ax.plot(ts, n_total2, label="Total", color="black") ax.set_yscale("log") ax.set_ylim(1e13, 1e17) idx = int(N*0.9) ax.text(ts[idx], n_freeze[idx], "Freeze-out", ha="right", va="top") idx = int(N*0.5) ax.text(ts[idx], n_sat[idx], "Saturation", ha="center", va="bottom") idx = np.argmin(np.abs(n_is-3e16)) ax.text(ts[idx+10], n_is[idx], "Intrinsic", ha="left", va="bottom") # ax.set_xlim(ts[0], ts[idx1+N//6]) # ax.legend() ax.set_xlabel(r"$1000/T\,[\si{\per\kelvin}]$") ax.set_ylabel(r"$n\,[\si{\per\cm^3}]$") return fig def test(): fig, ax = plt.subplots() N = 100 left = np.empty(N) left.fill(4.0) left = np.linspace(5.0, 0.0, N) right = np.empty(N) right.fill(2.5) right = np.linspace(3.0, 2.0, N) ax.plot(left, label="l",linestyle="dashed") ax.plot(right, label="r",linestyle="dashed") total_lin = blend_arrays(left, right, 40, 60, "linear") ax.plot(total_lin, label="lin") # total_tanh = blend_arrays(left, right, 40, 60, "tanh") # ax.plot(total_tanh) total_q = blend_arrays(left, right, 40, 60, "quadratic_simple") ax.plot(total_q, label="q") ax.legend() return fig # OPTICS @np.vectorize def fn(omega: float, omega_p:float, eps_infty:float, tau:float): eps_real = eps_infty * (1-(omega_p**2 * tau**2)/(1+omega**2 * tau**2)) eps_imag = eps_infty * omega_p**2 * tau/(omega*(1+omega**2 * tau**2)) eps_complex = (eps_real + 1j*eps_imag) n_complex: complex = np.sqrt(eps_complex) # n_real = n_complex.real # n_imag = n_complex.imag # return n_real, n_imag return n_complex def free_electrons_absorption(): # values taken from adv. sc. physics ex 10/1 fig, axs = plt.subplots(2, 2, figsize=size_formula_fill_default, sharex=True) omega_p = 30e12 # 30 THz eps_infty = 11.7 tau = 1e-6 omegas = omega_p * np.logspace(-1, 3, 1000, base=10, dtype=float) # omegas = omega_p * np.linspace(1e-1, 1e3, 1000) # print(omegas) n_complex = fn(omegas, omega_p, eps_infty, tau) n_real = n_complex.copy().real n_imag = n_complex.copy().imag R = np.abs((n_complex-1)/(n_complex+1)) alpha = 2*omegas*n_imag/3e8 for ax in axs[1,:]: ax.set_xlabel(r"$\omega/\omega_\text{p}$") ax.set_xscale("log") ax.set_xticks([1e-1,1e0,1e1,1e2,1e3]) ax.set_xticklabels([0.1, 1, 10, 100, 1000]) omegas_rel = omegas/omega_p axs[0,0].plot(omegas_rel, n_real) axs[0,0].set_yticks([0,1,2,np.sqrt(eps_infty)]) axs[0,0].set_yticklabels([0,1,2,r"$\epsilon_\infty$"]) axs[0,0].set_ylim(-0.2,0.2+np.sqrt(eps_infty)) axs[0,0].set_ylabel(r"$n^\prime(\omega)$") axs[1,0].plot(omegas_rel, n_imag) axs[1,0].set_ylabel(r"$n^{\prime\prime}(\omega)$") axs[1,0].set_yscale("log") axs[0,1].plot(omegas_rel, R) axs[0,1].set_ylabel(r"$R(\omega)$") axs[1,1].plot(omegas_rel, alpha) axs[1,1].set_ylabel(r"$\alpha(\omega)$") axs[1,1].set_yscale("log") return fig if __name__ == '__main__': # export(test(), "cm_sc_charge_carrier_density") export(schottky_barrier(), "cm_sc_devices_metal-n-sc_schottky") # export(charge_carrier_density(), "cm_sc_charge_carrier_density") # export(free_electrons_absorption(), "cm_sc_optics_free_electrons")