SEMITIP, VERSION 1, DOCUMENTATION

A program for computing the electric potential around a probe tip in proximity to a semiconductor, with circular symmetry. Prolate spheroidal coordinates are used in the vacuum, and a carefully chosen updating scheme is used to ensure stability of the iterative solution.
Version 1.0 - written by R. M. Feenstra, Carnegie Mellon University, Dec 2002

All routines are written in standard FORTRAN.

A complete description of the background theory of this program is contained in Ref. 1. As described there, the routine uses an iterative technique to solve Poisson's equation on grids in the vacuum and the semiconductor, using the boundary condition of continuity of the electric displacement field across the semiconductor surface. The values of the potential in the program are contained in arrays:
SEM(2,NR,NS) - potential in semiconductor
VAC(2,NR,NV) - potential in vacuum
VSINT(2,NR) - potential on semiconductor surface
where NR is the number of radial points, NS is the number of z-points in the semiconductor, NV is the number of z-points in the vacuum, and the dimension 2 is for two arrays, one with current value of the potential and the other with updated value.

The array dimensions are variables passed to SEMITIP2: NRDIM, NSDIM, and NVDIM are the array sizes in their declaration, on input NR, NS, and NV are the starting size of the arrays, and on output they are the final size of the arrays.

The grid sizes are automatically doubled after a solution to the problem at a given size is found to the desired accuracy. For example starting with grids of 32x32 in both the vacuum and the semiconductor, subsequent grids would be 64x64, 128x128, etc. For the Lth scaling step in the grid size the iterative solution to Poisson's equation is continued until all changes in the potential values are less than a user specified value of EP(L) or until the number of iterations in that scaling step exceeds ITMAX(L). The number of steps performed in the scaling is IPMAX.

CAUTION: The grid sizes are user specified parameters. An INCORRECT final result may be obtained is these values are too small. The calling program, SEMITIP_V1, uses as an estimate for the size of the grids in the lateral (DELR*NR) and vertical (DELS*NS) directions the value W0, where W0 is an estimate of the 1D depletion width for the problem (see Ref. 1). This value is appropriate for many depletion problems so long as the tip radius or the tip-sample separation are not too large. If those values are large, or for problems of accumulation, it may be necessary to use larger values for the simulation size. An appropriate value can always be determined by observing the change in the final potential values as the grid sizes (ISIZE parameter in SEMITIP_V1) is increased.

In the semiconductor, cartesian coordinates (r,z) are used, with:
r_i = (i - 0.5)*delr and z_j = -j*delz (z<0 in the semiconductor).
In the vacuum, prolate spheroidal coordinates are used consisting of families of confocal hyperbolas and ellipses. The hyperbolas are labelled by eta (eta varies between 0 and 1), and the ellipses by xsi (values of xsi are greater than 1). In the simplest case the probe tip corresponds to a particular hyperbola, having eta = eta_T. The (xsi,eta) coordinates are related to (r,z) coordinates by:
r = a*sqrt((xsi**2-1)*(1-eta**2)) and z = a*xsi*eta.
The constant 'a' is the distance from the surface to the focal point of the hyperbolas, as described in Ref. 1. The xsi grid is chosen to equal the r grid at the semiconductor surface
xsi_i =sqrt(1+(r_i/a)**2)
and the eta spacing is uniform, eta_j=j*deleta.

A sample calling program, SEMITIP_V1, is provided, which takes all its input from the file FORT.9. The first action within the calling program is to invoke EFFIND to locate the Fermi-level position in the semiconductor. Then the one-dimensional depletion width W0 and band bending PHI00 are computed. Then, a table of charge density values are complied, contained in the array 'rhoc' with user specified size 'ne'. The potential distribution is then computed by a call to SEMITIP2. Examples of outputting the potential values is included in SEMITIP_V1, as follows (output depends on the value of the input parameter IWRIT as specified in the input file FORT.9):

FORT.10 - gives the numerical results for the following quantities:
- tip radius of curvature
- tip-sample separation
- tip potential
- opening angle of tip shank
- PHI0 - thesurface potential at a point directly opposite the tip apex
- ratio of PHI0 to PHI00
FORT.11 - provides the potential along the central axis, as a function of z-distance (output for IWRIT>0)
FORT.12 - provides the potential along the surface, as a function of the radial distance from the central axis (output for IWRIT>0)
FORT.13 - gives the entire array of potential values (output for IWRIT>2)
FORT.16 - gives an exact copy of the output to the console
FORT.20, FORT.21, ... - contour lines of the potential (output for IWRIT>1)

Auxiliary routines for computing the charge density in the semiconductor are contained in SEMIRHO.

In SEMITIP2 the iterative solution to Poisson's equation takes place. Intermediate output of PHI0, the potential at a point on the surface directly below the tip apex, is provided if a value of IWRIT=1 is specified. Once the desired accuracy of EP(1) is reached for the solution, or the iterations exceed ITMAX(1), the arrays are all doubled and the iterative solution continues. This process continues until the number of doubling exceeds IPMAX. At that point the SEMITIP2 routine is exited, with the output potentials contained in SEM, VAC, and VSINT. Additional output values are:

PHI0 = the surface potential directly below the tip. Quadratic interpolation is used to form this value; it equals (9*PHI1-PHI2)/8 where PHI1 and PHI2 are the potentials at the first and second r-grid points on the surface
ETAT = eta-value for the actual tip
AT = a-value for the actual tip
ETA1 = eta-value for the tip used to construct the coordinates
A1 = a-value for the tip used to construct the coordinates
DELV = array of z-spacing values in the vacuum, corresponding to the first eta-value (i.e. distance from semiconductor surface to first grid line in the vacuum).

As discussed in Ref. 1, when reconstructing potential profiles in the vacuum after execution of SEMITIP2, the values of ETA1 and A1 should be used in the transformation equations to transform from (xsi,eta) to (r,z). On the other hand, to define the precise boundary of the probe tip itself, the values of ETAT and AT (together with SEP) should be used.

A user specified protrusion extending out from the hyperbolic probe tip is possible, and is defined by a P(R) function which must be declared in the calling program. In SEMITIP_V1 this function is declared as:

      REAL FUNCTION P(R)
      COMMON/PROTRU/RAD2
      P=0.
      IF (R.LT.RAD2) P=SQRT(RAD2**2-R**2)
      RETURN
      END

where the value of RAD2 is passed from the main program by the named COMMON block PROTRU. To use some shape for a protrusion this declaration of P(R) would have to be changed in the source code. It should be noted that the protrusion defined in P(R) is NOT included in the SEP variable which defines the tip-sample separation. Thus for the above example of the protrusion with radius RAD2, the actual tip-sample separation at the apex would be SEP-RAD2. Also, in order to obtain meaningful results from the program, the final value of DELR must be significantly less than RAD2 (unless RAD2=0, in which case no protrusion exists), and actually even the initial value of DELR should be less than RAD2 to ensure good convergence of the result. of DELR

Parameters involving the semiconductor doping are contained in a labelled COMMON block:

      COMMON/SEMI/EGAP,ED,EA,ACB,AVB,CD,CA,EPSIL,TK,IDEG

where the all the following parameters must be specified in the calling program:

      EGAP=BAND GAP (EV)
      ED=DONOR BINDING ENERGY (EV)
      EA=ACCEPTOR BINDING ENERGY (EV)
      ACB=CONDUCTION BAND EFFECTIVE MASS
      AVB=VALENCE BAND EFFECTIVE MASS
      CD=DONOR DENSITY (CM^-3)
      CA=ACCEPTOR DENSITY (CM^-3)
      EPSIL=DIELECTRIC CONSTANT OF SEMICONDUCTOR
      TK=TEMPERATURE (EV)
      IDEG=1 IF DEGENERATE DOPING, 0 OTHERWISE

SEMITIP2 is invoked by a call of the form:

      CALL SEMITIP2(SEP,RAD,SLOPE,ETAT,AT,ETA1,A1,VAC,SEM,VSINT,
     &DELV,DELR,DELS,NRDIM,NVDIM,NSDIM,NR,NV,NS,BIAS,IWRIT,
     &NE,ESTART,DELE,RHOC,EF,ITMAX,EP,IPMAX,PHI0,IERR)

where the INPUT parameters are:

      SEP=SEPARATION BETWEEN SEMICONDUCTOR AND END OF TIP (NM)
      RAD=RADIUS OF HYPERBOLOID WHICH FORMS TIP (NM)
      SLOPE=SLOPE OF TIP SHANK
      NRDIM=X-DIMENSION FOR VAC, SEM, AND VSINT ARRAYS
      NVDIM=Z-DIMENSION FOR VAC ARRAY
      NSDIM=Z-DIMENSION FOR SEM ARRAY
      NR=INITIAL VALUE OF X-DIMENSION FOR VAC, SEM, AND VSINT ARRAYS
      NV=INITIAL VALUE OF Z-DIMENSION FOR VAC ARRAY
      NS=INITIAL VALUE OF Z-DIMENSION FOR SEM ARRAY
      DELR=INITIAL R-SPACING IN VACUUM AND SEMICONDUCTOR
      DELS=INITIAL Z-SPACING IN SEMICONDUCTOR
      BIAS=BIAS OF TIP RELATIVE TO A POINT FAR INSIDE THE SEMICONDUCTOR
      IWRIT=1 IF WANT WRITTEN OUTPUT, 0 OTHERWISE
      NE=NUMBER OF POINTS IN CHARGE DENSITY TABLE
      ESTART=STARTING ENERGY IN TABLE
      DELE=ENERGY STEP SIZE IN TABLE
      RHOC=CHARGE DENSITY TABLE
      EF=FERMI-LEVEL POSITION
      ITMAX=ARRAY OF MAXIMUM ITERATION SPECIFICATIONS
      EP=ARRAY OF STOPPING CRITERIONS
      IPMAX=MAXIMUM NUMBER OF SCALING LOOPS FOR GRID DOUBLING

and the OUTPUT parameters are:

      ETAT=ETA PARAMETER FOR THE HYPERBOLIC TIP
      AT=LOCATION OF FOCUS FOR THE HYPERBOLIC TIP
      ETA1=ETA PARAMETER FOR THE HYPERBOLIC COORDINATE SYSTEM
      A1=A2=LOCATION OF FOCUS FOR THE HYPERBOLIC COORDINATE SYSTEM
      VAC=ARRAY OF POTENTIAL VALUES IN VACUUM
      SEM=ARRAY OF POTENTIAL VALUES IN SEMICONDUCTOR
      VSINT=ARRAY OF POTENTIAL VALUES AT SEMICONDUCTOR SURFACE
      DELV=ARRAY OF FINAL Z-SPACING VALUES IN VACUUM
      DELR=FINAL R-SPACING IN VACUUM AND SEMICONDUCTOR
      DELS=FINAL Z-SPACING IN SEMICONDUCTOR
      NR=FINAL VALUE OF X-DIMENSION FOR VAC, SEM, AND VSINT ARRAYS
      NV=FINAL VALUE OF Z-DIMENSION FOR VAC ARRAY
      NS=FINAL VALUE OF Z-DIMENSION FOR SEM ARRAY
      PHI0=SURFACE BAND BENDING BELOW TIP APEX
      IERR=ERROR INDICATION (=1 IF GRID SIZE TOO SMALL, =2 IF UNSTABILE SOLN)

References:
1. R. M. Feenstra, Electrostatic Potential for a Hyperbolic Probe Tip near a Semiconductor, published in J. Vac. Sci. Technol. B 21, 2080 (2003). For preprint, see http://www.cmu.edu/physics/stm/publ/52/.