9.1 A Pair-Diffusion Model in One Dimension

A five species phosphorus diffusion model introduced by Richardson and Mulvaney [12] has been developed. The species behavior can be formulated by the following system of coupled reaction diffusion equations

$$\frac{\partial V}{\partial t}=D_{V}\: \nabla ^{2}V-k^{E}_{for}\: P\: V+k^{E}_{rev}\: E-k_{bi}\: \left( V\: I-V^{eq}\: I^{eq}\right)$$ (6)

$$\frac{\partial I}{\partial t}=D_{I}\: \nabla ^{2}I-k^{E}_{for}\: P\: I+k^{E}_{rev}\: F-k_{bi}\: \left( V\: I-V^{eq}\: I^{eq}\right)$$ (7)

$$\frac{\partial P}{\partial t}=-k^{E}_{for}\: P\: V+k^{E}_{rev}\: E-k^{F}_{for}\: P\: I+k^{F}_{rev}\: F$$ (8)

$$\frac{\partial E}{\partial t}=D_{E}\: \nabla ^{2}E+k^{E}_{for}\: P\: V-k^{E}_{rev}\: E$$ (9)

$$\frac{\partial F}{\partial t}=D_{F}\: \nabla ^{2}F+k^{F}_{for}\: P\: I-k^{E}_{rev}\: F$$ (10)

where the solution vector $u = \left[ V,I,P,E,F \right]^T$ represents, respectively, the concentration of vacancies, interstitials, substitutional phosphorus as well as phosphorus-vacancy pairs (E-centers) and phosphorus-interstitial pairs (F-centers).

Furthermore, zero flux boundary conditions for all species are enforced everywhere except at the exposed surfaces where $V$, $I$, $P$, $E$, and $F$ are specified with

$V = V^{eq}\:$ and $I = I^{eq}$
$P = C^{*}\:$ (ambient gas concentration)
$E = \frac{k^{E}_{for}}{k^{E}_{rev}}\: P\: V \:$
$F = \frac{k^{F}_{for}}{k^{F}_{rev}}\: P\: I$

Using a finite element discretization within the analytical model interface the complete formulation of the discretized differential equation system is written in the code segment as shown below. The predefined variables like X,Y,Z (coordinates of the element points) and t (time) are initialized during runtime and can be used like any other variable within the model definition language.

MODEL PairDiff = [V,I,E,F,P]; 
{ 

############## Discretization ################################ 

N(xsi) = [1-xsi,xsi];        # the element shape function defined as 
                             # a function of xsi 
dNdxsi = [ -1  , 1 ];        # derivative of the shape function 

dxdxsi = dNdxsi * X;         # X is a predefined vector [X1,X2,...,Xn] 
                             # getting real size and coordinates 
                             # during runtime 

detJ = |dxdxsi|; 
dxsidx = Inv(dxdxsi);        # calculates the inverse matrix 

L = dNdxsi*dxsidx;           # gradient operator 
K = L^T*L;                   # Laplace operator where (^T) is a 
                             # transpose operator 
C = N(1/2)*N(1/2)^T;         # time operator using lumping 
                             # N(1/2) -> calculates the value 
                             #           of the shape function 
                             #           for xsi=0.5 



############## The Parameters ################################# 

Param k_for_e = 1.0E-14;     # in cm^3/s 
Param k_for_f = 1.0E-14;     # initializing default values 
Param k_bi = 1.0E-10;        # for parameters 

Param k_rev_e = 10;          # in 1/s 
Param k_rev_f = 12; 

Param Dv  = 1.0E-10;         # in cm^2/s 
Param Di  = 1.0E-09; 
Param De  = 1.0E-13; 
Param Df  = 2.0E-13; 

Param Veq  = 1.0E14;         # in 1/cm^3 
Param Ieq = 1.0E14; 
  

############## The Differential Equations ##################### 

dt = t.t0-t.t1;              # access to several time steps 
dV = V.t0-V.t1;              # of unknown and predefined 
dI = I.t0-I.t1;              # variables 
dE = E.t0-E.t1; 
dF = F.t0-F.t1; 
dP = P.t0-P.t1; 

i = 1..2;                    # running variable from 1 to 2 

PV[i]  = P[i] * V[i]; 
PI[i]  = P[i] * I[i]; 
VI[i]  = I[i] * V[i]; 
VIeq[i] = Veq * Ieq; 

resV = detJ * (Dv*K*V + C * (dV/dt-k_for_e*PV+k_rev_e*E-k_bi*(VI-VIeq))); 
resI = detJ * (Di*K*I + C * (dI/dt-k_for_f*PI+k_rev_f*F-k_bi*(VI-VIeq))); 
resE = detJ * (De*K*E + C * (dE/dt+k_for_e*PV-k_rev_e*E)); 
resF = detJ * (Df*K*F + C * (dF/dt+k_for_f*PI-k_rev_f*F)); 
resP = detJ *           C * (dP/dt-k_for_e*PV+k_rev_e*E-k_for_f*PI+k_rev_f*F); 

############## The Residual and its Derivative ##################### 

residuum = [[resV][resI][resE][resF][resP]];     # the residual vector 
jacobian = D([[V][I][E][F][P]],residuum^T)^T;    # auto derivative 
                                                 # of residual 
}

To define the necessary simulation parameters and to map the developed analytical pair-diffusion model onto a definite simulation domain including its boundaries the specifications of the Input & Control Interface (Section 3) have to look like:

# Where to read the grid from
Source=Pif
{
    Physical = lin2.pbf    # the input filename
    Logical  = BareWafer   # the logical name of the mesh
}
# Where to write the grid and the result to
Output = Pif
{
    Physical = lin2res.pbf # the output filename 
    Logical  = BareWafer   # the logical name of the new generated 
                           # output mesh
}

Time = 60                  # the first 60 seconds 
{
   Step    = 1.0E-3        # initial timestep
   Epsilon = 1.0E-3        # maximum error
   RejectionRatio = 0.5    # reduce step-size in case of non
                           # converging or epsilon error
   StretchRatio   = 1.5    # increase step-size in case of success
   MaxNewton      = 15     # stop after 15 iterations and reduce step 
                           # size 
}



Grid = Si_Region  # the name of the grid in the default input file
{ 
    Model = PairDiff(V,I,E,F,P)   # use the model PairDiff
    {                             # and map the quantities
        Value(V) = 1.0E14         # initialize the quantity V
        Value(I) = 1.0E14         # initialize the quantity I
        Value(E) = 0.0E0          # initialize the quantity E
        Value(F) = 0.0E0          # initialize the quantity F
        Value(P) = 0.0E0          # initialize the quantity P
    }
}

Boundary = Backside         # any name because auto boundary detection 
                            # is used
{
   Interface(Si_Region,Top) # Choose the top boundary of the grid with
                            # name Si_Region

   # Use Dirichlet boundaries for the defined quantities

   Model = Dirichlet(V){Value(V) = 1.0E14}     
   Model = Dirichlet(I){Value(I) = 1.0E14}     
   Model = Dirichlet(P){Value(P) = 5.0E20}     
   Model = Dirichlet(E){Value(E) = 5.0E19}     
   Model = Dirichlet(F){Value(F) = 4.16666E19} 
}

Finally the simulation results that show the distribution of all five species after 60 seconds (Fig. 10) are written to the file defined in the output section. Fig. 11 shows the result for the same model calculating the distribution after 600 seconds.

Figure 10: Dopant distribution after 60 seconds

Figure 11: Dopant distribution after 600 seconds