Simulateur numérique MXCAD1D

MXCAD1D est un code de simulation numérique développé par F. Martinez résolvant simultanément en une dimension les équations de Poisson et de continuité, en prenant en compte un taux de génération optique calculé via un solveur optique basé sur une technique de transfert matriciel. Le transport est décrit par un modèle drift diffusion. Ce code est implanté dans Matlab, offrant ainsi une grande souplesse en termes de modification des modèles physiques implantés, de l’interfaçage avec d’autres outils de simulation, d’interprétation de données expérimentales, etc...

Il n'existe pas pour l'instant d'interface "user friendly" et l'utilisation de mxcad1d nécessite une maitrise basique de Matlab. Si vous souhaitez utiliser ou adapter mxcad/atol à d'autres projets de simulation, n'hésitez pas à prendre contact : Cette adresse e-mail est protégée contre les robots spammeurs. Vous devez activer le JavaScript pour la visualiser. ou Cette adresse e-mail est protégée contre les robots spammeurs. Vous devez activer le JavaScript pour la visualiser.

Un exemple de simulation d'une cellule GaSb (Dark, 1 sun spectre AM1.5 et EQE/IQE) est disponible au téléchargement. Il suffit de décompresser l'ensemble des fichiers dans un dossier puis d'exécuter le script demo_gasb.m

 

télécharger la démo mxcad1d

 

La Figure 1 présente les variables utilisées dans la formulation des équations implantées dans MXCAD1D régissant le fonctionnement du dispositif :

Figure 1 : Définition des variables utilisées dans la formulation des équations.

Le potentiel électrostatique Ψ est défini par rapport l’énergie du vide Evac. c est l’affinité électronique du matériau. Le potentiel intrinsèque Ψint est relatif au niveau de Fermi intrinsèque Ei. Les potentiel des quasi-niveaux de Fermi des électrons (Φn) et des trous(Φp) sont relatifs aux quasi-niveaux de Fermi Efn et Efp.

En utilisant ces notations, le système d’équations différentielles à résoudre est donné par :

 

Jn et Jp sont les densités de courant des électrons et des trous, q la charge élémentaire, ND et NA les concentrations des impuretés ionisées, R le taux net de recombinaison, µn et µp la mobilité des trous et des électrons. En considérant des semi-conducteurs non dégénérés, la densité d’électrons n et de trous p et s’exprime par :

 

La densité intrinsèque de porteurs ni est défini par rapport aux densités d’état dans la bande conduction NC et de valence NV et du gap Eg du matériau :

 

Les conditions aux limites du système d’équation sont définies par les contacts aux extrémités du composant(x=0 et x=L). Pour un contact ohmique avec des vitesses de recombinaison infinies, sur lequel on applique une tension Vapp, les quasi-niveaux de Fermi sont égaux (le contact est à l’équilibre) :

 

Le contact étant à l’équilibre, les densités d’électrons n et de trous p vérifient :

 

En utilisant les relations , et , on obtient l’expression de Ψint au contact en fonction de Vapp :

 

Compte tenu de la relation entre Ψ et Ψint, Ψ au contact s’exprime par :

 

Discrétisation des équations Poisson/Drift Diffusion

Les équations sont discrétisés suivant la méthode des volumes finis. Les interfaces ont été modélisées en utilisant la méthode du double point [1], permettant ainsi la simulation d’hétérojonctions avec transport de porteurs par effet tunnel ou thermoïonique aux interfaces.

Le maillage de l’espace est présenté sur le Figure 2. Le dispositif à simuler est composé de domaines où le matériau est homogène.

 

Figure 2 : Maillage d’une hétérojonction pour la discrétisation des équations de Poisson/Drift Diffusion.

Pour un point i hors interface, dans un matériau homogène, l’équation de Poisson se discrétise suivant l’expression :

 
                                   

L’équation de continuité des électrons se discrétise sous la forme :


L’évaluation des courants aux demis intervalles est obtenue suivant l’approximation de Scharfetter-Gummel [2]. En utilisant les expressions utilisant le potentiel intrinsèque, on obtient [3] :


                                                

Une forme plus stable numériquement, notamment pour les faibles valeurs de yi+1-yi, est obtenue en réarrangeant l’expression  :


avec zcsch(z)=z∙csch(z).

Une expression similaire est utilisée pour le courant de trou :

 
         

Aux interfaces, les points k et k+1 ont la même abscisse xk=xk+1.Le potentiel Ψ(x) est continu aux interfaces d’où Ψkk+1. En l’absence de courant tunnel ou thermoïoniques aux interfaces, dans l’hypothèse du modèle de diffusion pour les hétérojonctions[5], les quasi niveaux de Fermi sont égaux.

Solveur optique

Le composant 1D défini pour la simulation électrique est constitué d’un empilement de couches dites « électriques », c’est-à-dire semi-conductrices ou isolantes. On peut ajouter des couches dites « optiques », notamment pour la simulation de couches anti-reflet. On obtient alors la structure sur laquelle on utilise un solveur optique basé sur une approche de transfert matriciel afin de calculer la génération optique (terme G dans l’équation ).

Le code Atol, développé par Frédéric Martinez et Stéphanie Parola, correspond à une implémentation dans Matlab d’un solveur basé sur une approche de transfert matriciel similaire à celui disponible dans le simulateur commercial sdevice [6]. Les paramètres matériaux sont l’indice de réfraction et le coefficient d’extinction en fonction de la longueur d’onde.

Exemple de simulation MXCAD1D/Atol

La Figure 3 présente la simulation d’une cellule GaSb constituée d’une jonction abrupte PN, avec un émetteur de 200 nm dopé NA = 1019 cm-3 et d’un absorbeur de 10 mm dopé ND = 1017 cm-3. Les paramètres et le script correspondant à cette simulation sont donnés en annexe A.

 

Figure 3 : Simulation MXCAD1D/Atol d’une cellule GaSb (a) sans couche anti-reflet, (b) avec couche anti-reflet SiNx de 115 nm d’épaisseur pour des durées de vie tn = tp = 10 ns et 100 ns.

On notera l’amélioration du photocourant de la cellule avec l’utilisation d’une couche anti-reflet SiNx de 115 nm d’épaisseur, ainsi que l’influence de la durée de vie des porteurs qui améliore la tension en circuit ouvert.

 

 

Annexe A

Script de simulation d’une cellule GaSb avec MXCAD1D

 

 q=1.6021765e-19; h=6.626069e-34; c=299792458; kb=1.38065e-23;

 clear options;

 options.infolevel=1;

 options.solver_cust=1;

 options.T=273.15+25;

 options.MB=0;

 T=options.T;

 

%% Paramètres matériaux du GaSb

 gasb.X=4.06;

 gasb.Eg=0.813-3.78e-4*T^2/(T+94);

 gasb.epsilon_r=15.7;

 gasb.Nc=4.0e13*T^1.5*1e6;

 gasb.Nv=3.5e15*T^1.5*1e6;

 gasb.taun=10e-9;   

 gasb.taup=600e-9;

 gasb.un=2600*1e-4;

 gasb.up=800*1e-4;

 gasb.Cn=5e-30*1e-12;

 gasb.Cp=5e-30*1e-12;

 gasb.Brad=1/10*8.5e-11*1e-6;

 load nk_GaSb.mat

 gasb.nk=nk;

 

%% Definition de la structure

 device={};

 k=1;

 a=1;

% Emetteur (type P)

 dev=gasb;

 dev.bound=[0 200e-9);

 dev.N=200;

 dev.Nd=3e16*1e6;

 dev.Na=1e19*1e6;

 [dev.un dev.up]=mobility_gasb(dev.Na, dev.Nd);

 device{k}=dev;k=k+1;

% Absorbeur (type N)

 dev=gasb;

 dev.bound=[0 10e-6]+max(device{k-1}.bound);

 dev.N=500;

 dev.Na=1e16*1e6;

 dev.Nd=1e17*1e6;

 [dev.un dev.up]=mobility_gasb(dev.Na, dev.Nd);

 xprobe=mean(dev.bound);

 device{k}=dev;k=k+1;

% Substrat (type N)

 dev=gasb;

 dev.bound=[0 500e-6]+max(device{k-1}.bound);

 dev.N=2000;

 dev.Na=3e16*1e6;

 dev.Nd=9e17*1e6;

 [dev.un dev.up]=mobility_gasb(dev.Na, dev.Nd);

 device{k}=dev;k=k+1;

 clear dev

 

%% Prise en compte d’une couche anti-reflet

 load AM15.mat

 lamb=350:5:1690;

 air.nk=[lamb' ones(length(lamb),1) zeros(length(lamb),1)];

 bound.first=air;

 bound.last=gasb;

 

 load nk_SiN_200C.mat

 siN.nk=max(0,nk);

 k=1;

 clear device_SiN

 dev=siN;

 dev.bound=[-115e-9 0];

 dev.N=100;

 device_SiN{k}=dev;k=k+1;

 

%% Solveur optique

 out=atol([device_SiN device],bound,lamb);

 R=out.R;

 

 options.Gatol=out;

 options.pvmode=0;

 options.xprobe=xprobe;

 options.pvmode=0;

 options.infolevel=0;

 options.Spect=[wavelength Global_tilt*0];

 options.pool=1;

 

%% Simulation sous obscurité

 [res_d0 jj0 vv0 jjprobe0]=MXCAD1Dcorev3(device,0,0,options);

 options.in=res_d0;

 vv=[-0.2:0.05:0.5];

 [res_d jj0 vv0 jjprobe0]=MXCAD1Dcorev3(device,vv,0,options);

 

%% Simulation sous 1 soleil

 vv=[-0.2:0.05:0.5];

 options.in=res_d0;

 options.Spect=[wavelength Global_tilt];

 [res jj vv1 jjprobe1]=MXCAD1Dcorev3(device,vv,0,options);

 

%% Simulation du rendement quantique

 clear JT

 options.pool=false;

 lrs=out.lamb(1:2:end);

 for i=1:length(lrs)

      options.lambda=lrs(i);

      options.in=res_d0;

     [res jj vv ]=MXCAD1Dcorev3(device,0,0,options);

     JT(i)=res.Jprobe;

 end

 R_lrs=spline(lamb,R,lrs);

 EQE=-JT*1240./lrs*100;

 IQE=EQE./(1-R_lrs);

 

Bibliographie

 

[1]     M. Grupen, K. Hess, G. Song, Simulation of transport over heterojunctions, Simul. Semicond. …. 4 (1991) 303–311.

[2]     D.L. Scharfetter, H.K. Gummel, Large-Signal Analysis of a Silicon Read Diode Oscillator, IEEE Trans. Electron Devices. 16 (1969) 64–77. doi:10.1109/T-ED.1969.16566.

[3]     W.L. Engl, H.K. Dirks, B. Meinerzhagen, Device modeling, Proc. IEEE. 71 (1983) 10–33. doi:10.1109/PROC.1983.12524.

[4]     K. Horio, H. Yanai, Numerical modeling of heterojunctions including the thermionic emission mechanism at the heterojunction interface, Electron Devices, IEEE Trans. 37 (1990).

[5]     Synopsys Inc., Senturus Device User Guide, in: 2013.