equilibrium_using_gummel.py

import numpy as np
import matplotlib.pyplot as plt

#defining derivative function using central difference, will need this to plot electric field
def gradient_over_linspace_central(arr, x):
    derivative_arr = np.zeros_like(x)
    for i in range(0,len(arr)):
        if(i==0):
            derivative_arr[i] = (arr[i+1] - arr[i])/(x[i+1]-x[i])
        elif(i == len(arr)-1):
            derivative_arr[i] = (arr[i-1] - arr[i])/(x[i-1]-x[i])
        else:
            derivative_arr[i] = (arr[i+1]-arr[i-1])/(x[i+1]-x[i-1])
    return derivative_arr

# Fundamental Constants
q = 1.602e-19        # C
kb = 1.38e-23        # J/K
T = 300              # K
ni = 1.5e10          # Intrinsic carrier concentration [1/cm^3]
Vt = kb * T / q      # Thermal voltage [V]

# Doping Values
Na = 1e16          # p-side doping [1/cm^3]
Nd = 1e17            # n-side doping [1/cm^3] keep in mind this behaves quite wierdly for higher doping values, but upon inspection its as if the functions are just sharper, which makes sense cuz if you have Na=1e16 and Nd=1e25, obviously the depletion width on the p side will be fatter than n, and by a lot, so yeah

# Built-in potential
Vbi = Vt * np.log(Na * Nd / (ni * ni))
print(f"Built-in potential: {Vbi:.3f} V")

# Simulation domain
x_min = -2e-4        # cm
x_max = 2e-4         # cm
n_points = 500
x = np.linspace(x_min, x_max, n_points) 
dx = x[1] - x[0]

# Set up spatially varying dielectric constant (epsilon)
# Example: Different epsilon for p-side and n-side with smooth transition couldnt really simulate this cuz the bandgap has some things to be taken care of
eps_Si = 1.05e-12    # F/cm for silicon
eps_p = eps_Si       # p-side dielectric constant
eps_n = eps_Si       # n-side dielectric constant

# Create spatially varying epsilon array
# Smooth transition using tanh function over ~0.1 micron
eps = np.where(x < 0, eps_p, eps_n)
# Optional: Add smooth transition region for simulation of heterojunctions
transition_width = 0.5e-5  # cm
eps = eps_p + (eps_n - eps_p) * 0.5 * (1 + np.tanh(x / transition_width))
#tanh transition for epsilon
#again, we havent done this simulation cuz it was too much working changing all constants like mu, tau bandgap, and what not
#also needed to account for the difference in bandgap at the heterojunction, so we just didnt do that part

# Calculate epsilon at cell interfaces (needed for finite differences)
eps_interface_right = np.zeros(n_points)
eps_interface_left = np.zeros(n_points)
for i in range(1, n_points-1):
    eps_interface_right[i] = 2 * eps[i] * eps[i+1] / (eps[i] + eps[i+1])  # Harmonic mean 
    eps_interface_left[i] = 2 * eps[i-1] * eps[i] / (eps[i-1] + eps[i])   # Harmonic mean
# Boundaries
eps_interface_right[0] = eps[0]
eps_interface_left[0] = eps[0]
eps_interface_right[-1] = eps[-1]
eps_interface_left[-1] = eps[-1]

# Set up doping profile
dop = np.where(x < 0, -Na, Nd)

# Initialize potential (simple initial guess)
fi = np.where(x < 0, -Vbi/2, Vbi/2)
#fi = (np.tanh(50*x)/2 + 0.5)*Vbi
#this tanh thing would have been a smoother guess but sir asked us to keep the initial guess as a inverted signum-like function 

# Normalize everything
dop_norm = dop / ni
fi_norm = fi / Vt

# Normalize epsilon and dx for each point
eps_norm = eps / (q * ni)
eps_right_norm = eps_interface_right / (q * ni)
eps_left_norm = eps_interface_left / (q * ni)
dx_sq = dx * dx

# Jacobi iteration solver
delta_acc = 1e-4
max_iter = 10000
omega = 0.15  # Relaxation parameter (SOR-like) for faster convergence

print("Starting Jacobi iteration...")

for iteration in range(max_iter):
    fi_old = fi_norm.copy()
    
    # Jacobi iteration: update all points simultaneously using old values
    for i in range(1, n_points-1):
        # Poisson equation with varying epsilon:
        # d/dx(eps * dfi/dx) = -q*ni*(exp(fi) - exp(-fi) - dop_norm)
        
        # Discretized form:
        # (eps_right * (fi[i+1] - fi[i]) - eps_left * (fi[i] - fi[i-1])) / dx^2 
        # = -q*ni*(exp(fi) - exp(-fi) - dop_norm)
        
        # Linearized for Newton-Raphson:
        # (eps_right + eps_left)/dx^2 * fi[i] = 
        #   (eps_right*fi[i+1] + eps_left*fi[i-1])/dx^2 
        #   - q*ni*(exp(fi) - exp(-fi) - dop_norm - fi*(exp(fi) + exp(-fi)))
        
        exp_fi = np.exp(fi_old[i])
        exp_minus_fi = np.exp(-fi_old[i])
        
        # Coefficient for fi[i]
        coeff_i = (eps_right_norm[i] + eps_left_norm[i]) / (Vt * dx_sq) + \
                  (exp_fi + exp_minus_fi)
        
        # RHS terms
        neighbor_term = (eps_right_norm[i] * fi_old[i+1] + \
                        eps_left_norm[i] * fi_old[i-1]) / (Vt * dx_sq)
        
        nonlinear_term = exp_fi - exp_minus_fi - dop_norm[i] - \
                        fi_old[i] * (exp_fi + exp_minus_fi)
        
        # Jacobi update
        fi_new = (neighbor_term - nonlinear_term) / coeff_i
        
        # Apply relaxation (optional, improves convergence)
        fi_norm[i] = omega * fi_new + (1 - omega) * fi_old[i]
    
    # Boundary conditions (Dirichlet - fixed potential at contacts)
    fi_norm[0] = fi_old[0]   # Left boundary
    fi_norm[-1] = fi_old[-1]  # Right boundary
    
    # Check convergence
    delta_max = np.max(np.abs(fi_norm - fi_old))
    
    if iteration % 1000 == 0:
        print(f"Iteration {iteration}: max change = {delta_max:.2e}")
    
    if delta_max < delta_acc:
        print(f"Converged after {iteration+1} iterations")
        break
else:
    print(f"Warning: Did not converge after {max_iter} iterations. Final change: {delta_max:.2e}")

# Calculate physical quantities
n_conc = ni * np.exp(fi_norm)  # 1/cm^3
p_conc = ni * np.exp(-fi_norm)  # 1/cm^3
#fi_norm += Vbi/(2*Vt)
potential = fi_norm * Vt  # V

E_field = -gradient_over_linspace_central(potential, x)  # V/cm

# Band diagram (simplified)
Eg = 1.12  # Silicon bandgap [eV]
Ec = -potential  # Conduction band
Ev = Ec - Eg  # Valence band
Ei = Ec - Eg/2  # Intrinsic Fermi level
Ef = 0  # Fermi level (reference)

# Convert x to microns for plotting
x_um = x * 1e4

# Create plots
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Plot 1: Potential
axes[0, 0].plot(x_um, potential + Vbi/2, 'b', linewidth=2)
axes[0, 0].set_xlabel('Position [μm]')
axes[0, 0].set_ylabel('Potential [V]')
axes[0, 0].set_title('Electrostatic Potential')
axes[0, 0].grid(True, alpha=0.3)
axes[0, 0].axvline(x=0, color='k', linestyle='--', alpha=0.3)

# Plot 2: Carrier Concentrations
axes[0, 1].semilogy(x_um, n_conc, 'r', linewidth=2, label='Electrons (n)')
axes[0, 1].semilogy(x_um, p_conc, 'b', linewidth=2, label='Holes (p)')
axes[0, 1].set_xlabel('Position [μm]')
axes[0, 1].set_ylabel('Carrier Concentration [1/cm³]')
axes[0, 1].set_title('Electron and Hole Concentrations')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
axes[0, 1].axvline(x=0, color='k', linestyle='--', alpha=0.3)

# Plot 3: Band Diagram
axes[1, 0].plot(x_um, Ec, 'r', linewidth=2, label='Ec (Conduction)')
axes[1, 0].plot(x_um, Ev, 'b', linewidth=2, label='Ev (Valence)')
axes[1, 0].plot(x_um, Ei, 'g--', linewidth=1.5, label='Ei (Intrinsic)')
axes[1, 0].axhline(y=Ef, color='k', linestyle='-', linewidth=2, label='Ef (Fermi)')
axes[1, 0].set_xlabel('Position [μm]')
axes[1, 0].set_ylabel('Energy [eV]')
axes[1, 0].set_title('Energy Band Diagram')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].axvline(x=0, color='k', linestyle='--', alpha=0.3)

# Plot 4: Electric Field and Epsilon
ax4 = axes[1, 1]
color1 = 'g'
ax4.plot(x_um, E_field, color=color1, linewidth=2)
ax4.set_xlabel('Position [μm]')
ax4.set_ylabel('Electric Field [V/cm]', color=color1)
ax4.set_title('Electric Field and Dielectric Constant')
ax4.tick_params(axis='y', labelcolor=color1)
ax4.grid(True, alpha=0.3)
ax4.axvline(x=0, color='k', linestyle='--', alpha=0.3)

# Add epsilon on secondary y-axis
ax4_twin = ax4.twinx()
color2 = 'm'
ax4_twin.plot(x_um, eps * 1e12, color=color2, linewidth=2, linestyle='--')
ax4_twin.set_ylabel('ε [pF/cm]', color=color2)
ax4_twin.tick_params(axis='y', labelcolor=color2)

plt.tight_layout()
plt.show()