def _bohm_speed(Te_eV, m_ion):
"""Bohm velocity sqrt(e*Te / m_ion)."""
return np.sqrt(max(Te_eV, 0.01) * _E_CHARGE / m_ion)
def _k_rate(Te, A_coeff, E_thresh):
"""Electron-impact rate coefficient k = A * exp(-E_thresh / Te) [m^3/s]."""
return A_coeff * np.exp(-E_thresh / max(Te, 0.01))
def plasma_etch_0d_rhs(t, y):
n = np.maximum(y[:14], _FLOOR)
Te = max(y[14], 0.01)
# y[15] = wall_ions (cumulative), y[16] = theta_F, y[17] = theta_Si
theta_F = np.clip(y[16], 0.0, 1.0)
theta_Si = np.clip(y[17], 0.0, 1.0)
n_e = n[9]
# --- electron-impact reactions ---
# e + CF4 -> CF3 + F + e (dissociation)
k_diss_cf4 = _k_rate(Te, 1e-14, 12.5)
R_diss = k_diss_cf4 * n_e * n[0]
# e + CF3 -> CF2 + F + e
k_diss_cf3 = _k_rate(Te, 8e-15, 10.0)
R_diss2 = k_diss_cf3 * n_e * n[1]
# e + CF2 -> CF + F + e
k_diss_cf2 = _k_rate(Te, 5e-15, 11.0)
R_diss3 = k_diss_cf2 * n_e * n[2]
# e + CF -> C + F + e
k_diss_cf = _k_rate(Te, 3e-15, 13.0)
R_diss4 = k_diss_cf * n_e * n[3]
# e + Ar -> Ar* + e (excitation)
k_exc = _k_rate(Te, 5e-16, 11.55)
R_exc = k_exc * n_e * n[7]
# e + CF4 -> CF4+ + 2e (ionization)
k_ion_cf4 = _k_rate(Te, 2e-15, 16.0)
R_ion_cf4 = k_ion_cf4 * n_e * n[0]
# e + CF3 -> CF3+ + 2e
k_ion_cf3 = _k_rate(Te, 1.5e-15, 12.0)
R_ion_cf3 = k_ion_cf3 * n_e * n[1]
# e + Ar -> Ar+ + 2e
k_ion_ar = _k_rate(Te, 1e-15, 15.76)
R_ion_ar = k_ion_ar * n_e * n[7]
# Ar* + CF4 -> CF3 + F + Ar (Penning)
k_penning = 5e-16 # m^3/s
R_penning = k_penning * n[8] * n[0]
# F + F + M -> F2 + M (recombination, termolecular approx as bimolecular)
k_f_recomb = 1e-33 * (n[0] + n[7]) # effective second-order via third body
R_f_recomb = k_f_recomb * n[4] * n[4]
# --- wall losses ---
v_th_F = np.sqrt(8.0 * _K_B * _T_GAS / (np.pi * 19.0 * 1.661e-27))
wall_loss_rate = 0.25 * v_th_F * _A / _V # 1/s for neutrals
# radical wall loss (sticking ~ 0.01 for CFx, 0.1 for F on SiO2)
gamma_cfx = 0.01
gamma_f = 0.1
# Ion wall loss via Bohm flux
u_b_cf4p = _bohm_speed(Te, _M_CF4P)
u_b_cf3p = _bohm_speed(Te, _M_CF3P)
u_b_arp = _bohm_speed(Te, _M_ARP)
h_L = 0.5 # sheath edge-to-bulk density ratio (simplified)
ion_wall_cf4p = h_L * n[10] * u_b_cf4p * _A / _V
ion_wall_cf3p = h_L * n[11] * u_b_cf3p * _A / _V
ion_wall_arp = h_L * n[12] * u_b_arp * _A / _V
total_ion_wall = ion_wall_cf4p + ion_wall_cf3p + ion_wall_arp
# Ar* wall quenching
v_th_Ar = np.sqrt(8.0 * _K_B * _T_GAS / (np.pi * 40.0 * 1.661e-27))
ar_star_wall = 0.25 * v_th_Ar * _A / _V * n[8]
# --- surface etch ---
# F atom flux to wafer
flux_F = 0.25 * v_th_F * n[4] # /m^2/s
etch_rate_coeff = 1e-2 # reaction probability
R_etch = etch_rate_coeff * flux_F * _A_WAFER / _V * theta_Si
# 4 F consumed per SiF4 produced
R_F_consumed = 4.0 * R_etch * _V / max(_A_WAFER, 1e-30)
# --- species balance (dn/dt) ---
dy = np.zeros(18)
# CF4
dy[0] = -R_diss - R_ion_cf4 - R_penning
# CF3
dy[1] = R_diss + R_penning - R_diss2 - R_ion_cf3 - gamma_cfx * wall_loss_rate * n[1]
# CF2
dy[2] = R_diss2 - R_diss3 - gamma_cfx * wall_loss_rate * n[2]
# CF
dy[3] = R_diss3 - R_diss4 - gamma_cfx * wall_loss_rate * n[3]
# F
dy[4] = R_diss + R_diss2 + R_diss3 + R_diss4 + R_penning - 2.0 * R_f_recomb \
- gamma_f * wall_loss_rate * n[4] - R_F_consumed
# F2
dy[5] = R_f_recomb
# C
dy[6] = R_diss4 - gamma_cfx * wall_loss_rate * n[6]
# Ar
dy[7] = -R_exc - R_ion_ar + ar_star_wall + R_penning
# Ar*
dy[8] = R_exc - ar_star_wall - R_penning
# electrons (quasi-neutrality source)
dy[9] = R_ion_cf4 + R_ion_cf3 + R_ion_ar - total_ion_wall * _V / max(n_e, _FLOOR)
# Technically dn_e/dt = sum(ionization) - sum(wall loss), same total rate
dy[9] = R_ion_cf4 + R_ion_cf3 + R_ion_ar - (ion_wall_cf4p + ion_wall_cf3p + ion_wall_arp)
# CF4+
dy[10] = R_ion_cf4 - ion_wall_cf4p
# CF3+
dy[11] = R_ion_cf3 - ion_wall_cf3p
# Ar+
dy[12] = R_ion_ar - ion_wall_arp
# SiF4 (etch product)
dy[13] = R_etch
# --- electron energy balance ---
# P_abs / (n_e * V) heats electrons; collisions cool them
E_loss_per_ion = 50.0 # eV lost per ionization event (including excitation)
E_loss_coll = (12.5 * R_diss + 10.0 * R_diss2 + 11.0 * R_diss3
+ 13.0 * R_diss4 + 11.55 * R_exc
+ E_loss_per_ion * (R_ion_cf4 + R_ion_cf3 + R_ion_ar))
# dTe/dt ~ (2/3) * [P_abs/(n_e*V*e) - Te * (ionization_freq) - E_loss / n_e]
P_per_electron = _P_ABS / (max(n_e, _FLOOR) * _V * _E_CHARGE)
cool_rate = E_loss_coll / max(n_e, _FLOOR)
dy[14] = (2.0 / 3.0) * (P_per_electron - cool_rate)
# --- cumulative ion flux (diagnostic) ---
dy[15] = (ion_wall_cf4p + ion_wall_cf3p + ion_wall_arp) * _V / _A
# --- surface coverage ---
# theta_F increases with F flux, decreases with ion-assisted desorption
F_ads_rate = gamma_f * flux_F * (1.0 - theta_F) / 1e19 # normalized
ion_flux_wall = h_L * (n[10] * u_b_cf4p + n[11] * u_b_cf3p + n[12] * u_b_arp)
ion_desorb = 1e-3 * ion_flux_wall * theta_F / 1e19
dy[16] = F_ads_rate - ion_desorb
# theta_Si: exposed Si decreases as fluorine covers it
dy[17] = -F_ads_rate + ion_desorb
return dy