计算科学 - 自洽 Poisson-Schrodinger 求解器的问题 - 吾爱随笔录

我目前正在为器件异质结构（高迁移率电子晶体管）编写自洽薛定谔-泊松求解器。该算法基于此期刊⁽¹⁾。然而，我遇到了一些困难。似乎我可以在基板（即结构的右端）有 Dirichlet BC 或 Neumann BC。下面分别是 Dirichlet 和 Neumann BC 的快照，其中 $x_n=0$ 首先，和 $\partial x/\partial n = 0$ 对于第二个情节。请注意，有问题的边界 $\partial \Omega$ 位于情节的右端。

我真的对有边界很感兴趣 $\partial\Omega$ 就像下面的情节。

_{(1) 使用非均匀网格的薛定谔-泊松方程的自洽解。谭，I-H. and Snider, GL and Chang, LD and Hu, EL, Journal of Applied Physics, 68, 4071-4076 (1990)}

下面概述了我正在使用的离散化方案。

我的计划是模拟一个 AlGaAs/GaAs HEMT 堆栈，它看起来与我上面引用的期刊中的堆栈几乎相同（堆栈在图 3 中）。在我的某些层中只有几埃的差异。我从一个简单的试验潜力开始 $V_{t r i a l} (x) = - q ϕ (x) + Δ E_{c} (x)$ $V_{trial}(x) = -q\phi(x) + \Delta E_c(x)$ 最初，这是一条简单的线，从 0.7 eV 开始，阶梯函数为 0.2 eV，跨越 AlGaAs 区域。
然后我使用有限差分方案求解薛定谔方程。

函数薛定谔( $V_{trial}$ ) 首先被调用。

function [prob_Density, eig_Vals] = schrodinger(phi)
p = parameters;
n = p.grid_points;
t0 = p.hbar^2 / (2*p.m_AlGaN*p.dz^2);

A = zeros(n,n);

% Define first and last rows
A(1,1) = 2*t0; A(1,2) = -2*t0; A(n-1,n) = -2*t0; A(n,n) = 2*t0;      

for i=2:(n-1)  % Create Hamiltonian
    A(i,i) = 2*t0 + phi(i);
    A(i+1,i) = -t0; A(i,i+1) = -t0; 
    A(i-1,i) = -t0; A(i,i-1) = -t0;
end

[eig_states, eig_Vals] = eig(A);  % Solve for eigensolutions

prob_Density = eig_states.*conj(eig_states);  % Calculate probability function

end  % End of function

使用步骤 2 中的结果“prob_Density”，然后我使用以下方法求解当前分布 $n (x) = \sum_{k = 1}^{m} ψ_{k}^{*} (x) ψ_{k} (x) n_{k}$ $n(x) = \sum_{k=1}^{m} \psi_k^* (x) \psi_k (x) n_k$ 在哪里 $n_{k} = \frac{m^{*}}{π ℏ^{2}} \int_{E_{k}}^{\infty} \frac{1}{1 + e^{(E - E_{F}) / k T}} d E$ $n_k = \frac{m^*}{\pi \hbar^2} \int_{E_k}^{\infty} \frac{1}{1+e^{(E-E_F)/kT}} dE$

接下来，我使用 $n(x)$ 并求解泊松方程。我再次将有限差分方案与泊松方程一起使用。我的 x=0 的狄利克雷边界条件是在肖特基势垒高度定义的。我的第二个边界是 Neumann 类型，其中导数为零。

“电荷密度”的代码（即 $n(x)$ ) 和泊松求解器分别写在下面。

%#### FIND CHARGE DISTRIBUTION ########
function charge_density = n_x(prob_Density, eig_Vals, Ef, final_Eigstate)

p = parameters;

n = p.grid_points;
n_array = zeros(n,1);
const = (p.m_AlGaAs ./ (pi*p.hbar^2));

fermi_dirac = @(E) 1./(1+exp((E-Ef)./p.Vt)); % Unit eV

for eig_state = 1:final_Eigstate

    Ek = eig_Vals(eig_state,eig_state); % First,second,third .. eigenenergy
    nk = const .* integral(fermi_dirac, Ek, inf)*p.q;   %  [m]^-2
    n_x_arr = prob_Density(:,eig_state).*nk;    % [m]^-1 * [m]^-2 -> [m]^-3
    n_array = n_array + n_x_arr;
end

charge_density = n_array; % Electron volumetric density [m]^-3


end % End program


%######### START POISSON SOLVER ###########
function phi = poisson(n_x)

first_boundary = 0.7; % eV
second_boundary = 0;    

p=parameters;
n = p.grid_points;

% Create Poisson matrix elements 
A = zeros(n,n);

for i = 2:n-1
   A(i,i) = -2; 
   A(i+1,i) = 1; 
   A(i-1,i) = 1; 
   A(i,i+1) = 1; 
   A(i,i-1) = 1;
end

A(1,1) = -2; 
A(1,2) = 1;
A(n,n) = -1; % -1 a result of Neumann BC
A(n,n-1) = 1;

F = zeros(n,1);

%% RHS of matrix equation (the 'B' n x 1 term) AU=B, U=A\B
h = p.dz;

for i = 2:n-1 
   f = p.q*(p.ND(i) - n_x(i))./p.eps_hemt(i);
   F(i) = f*h^2;
end

f1 = p.q*(p.ND(1) - n_x(1))./p.eps_hemt(1);
fn = p.q*(p.ND(n) - n_x(n))./p.eps_hemt(n);

F(1) = f1*h^2 - first_boundary;

F(n) = fn*h^2 - second_boundary;


phi = (A\F);                      %  Voltage profile

end
%##### END POISSON SOLVER

最后，我使用泊松求解器的结果并返回步骤 2。我再次运行函数 Schrodinger( $V_{new}$ ) 并遍历这个过程，直到找到一个自洽的解决方案。这种方法将我带到您在上面看到的第一个图，我的目标是模拟第三个图。