A Box Method Discretisation for the DriftDiffusion Equations in a Magnetic Field
Introduction
The problem of discretising the semiconductor driftdiffusion equations in the presence of a uniform magnetic field using the Box Method has previously been addressed in the literature. For the case of a rectangular nonuniform grid in 2D see[1], and for a general triangular mesh see[2]. These methods have also been extended to 3D [3],[4]. We have implemented and analysed these methods and found them all to be flawed in some way. In this paper we describe the established discretisations and elaborate on their flaws. We then briefly describe the new discretisation developed by Silvaco, and give some examples of its use.
Critique
The relationship between current in the presence of a magnetic field J_{B} and the current with no magnetic field J_{0} is given by
where µ* = rµ, with r a dimensionless scattering factor and µ the carrier mobility. If we let a = µ*B_{x},b = µ* B_{y}, c = µ*B_{z}, then we can write the inverse of M as
The domain of validity of approximation of equation (1) is that
The method of solution of Allegretto [2], is to calculate the current fluxes along the sides of the mesh elements (triangles in 2D, tetrahedra or prisms in 3D), using the usual ScharfetterGummel discretisation. For each apex of the mesh element, one constructs a matrix U, whose rows are the normalised edge vectors of the sides connected to the apex. This allows the array of edge current densities to be transformed to an orthonormal basis, the matrixM^{1}is applied to them and the resulting array is projected back to the sides of the mesh element. The contribution to the flux integral over the Voronoi cell from those edges is obtained by using the usual coupling areas, so the flux contribution is
where is the set of edge discretisations (not generally orthogonal), and is the set of coupling areas. The strength of this discretisation is that it has excellent convergence properties. It fatal weakness is that it fails to conserve current at the level of the Voronoi cells, leading to a lack of terminal current continuity. In the Appendix this lack of current continuity is proven for the 2d case.
An alternative way of doing the current discretisation is to obtain a current vector for a mesh element which gives the best fit in a least squares sense to edge current fluxes obtained using the usual ScharfetterGummel discretisation. This forms the basis of the VRdiscretisation method of Patil [5], in a nonmagnetic context. Having obtained the least squares currents in an orthonormal basis, one can then apply the matrix M^{1} and project the resulting current onto the mesh element edges in order to assemble the current flux using the Box method. Alternatively, one can decompose the matrix M^{1} into I + P, where I is the identity matrix and
is the matrix containing all the magnetic terms. This form transforms the discretisation into the usual nonmagnetic part that does not require the least squares current, plus a magnetic perturbation term which does. These least squares current methods have reasonable convergence properties, but can exhibit unphysical spatial oscillations in the converged values of current density. In extreme cases this results in negative carrier densities, which are clearly unphysical.
Although previous versions of ATLAS contain the discretisations described above, a new better discretisation has been developed by Silvaco engineers. The stability of discretisations of the current continuity equations in the presence of a Magnetic field has been discussed by He for Finite Element methods [6]. He introduces the idea of modifying the integrating factor used to obtain the ScharfetterGummel discretisation by magnetic terms. We extended this to the Box Method and reduce the discretisation to the usual ScharfetterGummel discretisation, but with two modifications. The first is that there is a prefactor equal to (1 + a^{2} + b^{2} + c^{2})^{1}, and the second is that the argument to the Bernoulli function has the following term added to it
where (Δx, Δy, Δz) is the spatial vector describing the edge, and K_{B} is the Boltzmann constant. The quantity φ is the electron quasiFermi level for electron current and the hole quasiFermi level for hole current. The Magnetic2D and Magnetic3D modules in the ATLAS device simulator use this discretisation. Magnetic3D also has an alternative discretisation, but we restict ourselves in this paper to the one described.
To illustrate the problems with the discretisations employing least squares estimates of current we created a simple 2D nchannel MOSFET structure biased into the saturation region with the channel in strong inversion. A magnetic field of intensity 2 Tesla was applied perpendicular to the plane of the device and the carrier concentrations obtained. The device, with the electron concentration obtained with the least squares method is shown in Figure 1a. A vertical cutline was taken through the source side depletion region and the carrier densities are shown in Figures 1b and 1c. There are some obvious unphysical oscillations with the leastsquares current based method, but the Silvaco method avoids these.


Figure 1: Comparison of carrier densities in a MOSFET with a 2T magnetic field. 
Simulation of Some Magnetic Devices
The first device simulated is a simple split drain device, also known as a MAGFET. It is a 100 x 125 micron rectangle of silicon with a contact extending along one short side, and two drains on the other short side with a separation of 10 microns. The input deck used to create the structure and define doping is in Figure 2a. If no magnetic field is applied, then clearly the drain currents are the same. If a magnetic field is applied perpendicular to the plane of the device, the current is deflected and a difference in the drain currents is obtained. The theory says that the Normalised Differential Current should be proportional to the applied magnetic field, B. The Normalised Differential Current is defined as
The splitdrain device, with a uniform donor density of 10^{15}cm^{3}, was simulated using Magnetic2D with a source drain applied bias of 0.5 Volts. From the loglog plot in Figure 2b, one can see that there is a linear relationship over the range of Magnetic field values considered. Another simulation was carried out on the same device structure, but with ntype doping under the split drains, ptype doping under the source and a uniform ntype background doping. The ATLAS syntax for this NnP diode structure is shown in Figure 2a. For this simulation the magnetic field was maintained at 1 Tesla, and a forward bias was applied between the source and the splitdrains. The normalised differential current as a function of applied bias is shown in Figure 2c, along with results obtained from the same device but with uniform ntype and ptype doping levels respectively. The normalised differential current decreases with applied bias, and at high bias it changes sign. The Allegretto discretisation also gives this effect, but is not reliable because the lack of current continuity is approximately 80% of the difference in currents between the two drains. The behaviour of the Normalised Differential Current arises from the bipolar nature of the carrier transport and the interaction between Hall effects and current deflections.


Figure 2: Simulations of a 2D MAGFET device. 
The next structure considered is a 3dimensional Hall cross structure. This consists of an ntype Greek cross shape embedded in a ptype substrate. Current is passed between the source and drain in the ydirection, a magnetic field of 1 Tesla is applied in the zdirection, and the Hall Voltage is measured between the Probe contacts in the xdirection. It is usual to add either an isolated conductive top plate to the ntype active region, or to have a shallow p+ layer diffused into the top of the ntype region [7]. We simulate the effect of having a shallow doped p+ layer on the top surface. The structure which we simulated is shown in Figure 3a, with the diffused p+ layer shown as the red square. The layer has an acceptor concentration of 10^{18} cm^{3} and is 0.1 microns deep.
Current is passed between the anode and cathode, and the opencircuited Hall sensors measure the potential difference arising from the Hall effect. In Figure 3b we see the Voltages at the probe electrodes as a function of anode bias, and in Figure 3c the Hall Voltage generated between the probe electrodes. The effect of the p+ diffused layer is to restrict the values of the probe potential. This is because much of applied bias (25 volts for this example) is dropped across the depletion region at the diffused layer edge at y=5 microns. This is shown in Figure 3d, where it can also be seen that the structure without the diffused layer has a higher electric field over most of the length of the active region, giving a higher current for the same applied bias. It also means that the currentvoltage relationship is more linear for the structure without the p+ diffused structure. This explains the differences seen in Figures 3b and 3c at higher values of anode bias.
A notable feature of the Figure 3c is that the device with the p+ diffused region has a larger Hall Voltage at lower values of applied bias. Additionally, the current is lower and so the device sensitivity in terms of Hall Voltage per Amp is greater as shown in Figure 3e. This effect depends on the details of the geometry and doping profiles, and MAGNETIC3D is an ideal package to try to optimise device sensitivity and linearity.
Figure 3: Simulations of a 3D Hall Cross Device. 
Conclusions
We have developed a new discretisation scheme for Magnetic problems. We have shown that it compares favourably with more established discretisations, and given some examples of its use. It is available for 2D and 3D simulations in ATLAS versions later than 3.16.3.R.
Appendix: Proof that Allegretto’s method fails to preserve current continuity
We consider a twodimensional rectangular piece of semiconductor with contacts at each end. A regular mesh for the discretisation of the current continuity and Poisson’s equations is created, having three horizontal grid lines perpendicular to the contacts and an arbitrary number greater than three parallel to the contacts. This is shown in the upper part of Figure 4. We consider the Voronoi cells surrounding any three mesh points in the same column, with label i, as shown in the lower part of Figure 4. From the Appendix of [2], the integral of current flux over a Voronoi cell can be partitioned into a normal component and a tangential component
The tangential component can be viewed as a line integral around the Voronoi cell. Allegretto endeavours to prove that this line integral is zero for any mesh points not lying on an exterior surface of the device. Calculations with specific meshes seem to support this. For mesh points on the surface, the value of the line integral can be calculated for the top surface to be
and for the bottom surface it is
where n is the carrier density at the surface node, µ is the mobility , and φ is the quasiFermi level. The first subscript is the row number and the second subscript is the column number. We notate these values as Δ_{1} and Δ_{2} respectively. For the three Voronoi cells shown in the lower part of Figure 4, we can explicitly impose the condition that the flux integral of J_{B} over each cell is zero. Thus we have for the cells from top to bottom
where a,b, and v are normal fluxes of J_{0}, that is from the first term on the right hand side of equation 4. Thus the difference between current flux entering and leaving the cells is
If we consider the entire device of N columns, the difference between anode and cathode currents will be
In the special case of constant carrier density, n, and mobility µ, then this simplifies to
and because the quasiFermi levels on the contacts are the value of the applied bias there, we have
Thus the flux difference vanishes and we have current continuity.
In the case where the carrier density is nonuniform, the value of expression 7 will not usually be zero because it becomes a weighted sum of quasiFermi level differences. Only the sum of the differences in quasiFermi levels is a conserved quantity. So we lose overall current continuity. This is exactly the behaviour we observed with the Allegretto discretisation.
It appears that the current continuity problem can be overcome by specifying an additional normal current flux boundary condition. But in order to be guaranteed to work, it would need to be equal and opposite to the Δ terms. In this case the discretisation would become a nonmagnetic one, modified by the term , which merely introduces quadratic magnetoresistance terms and does not give a Hall field.
Figure 4: Schematic 2D device for continuity analysis. 
References
 L. Andor, H. P. Baltes, A. Nathan and H.G. SchmidtWeinmar, ‘Numerical Modelling of MagneticFieldSensitive Semiconductor Devices’, IEEE Transactions on Electron Devices, Vol. ED32, No. 7, pp. 12241230, (1985).
 W. Allegretto, A. Nathan, H. Baltes, ‘Numerical Analysis of MagneticFieldSensitive Bipolar Devices, IEEE Transactions on Computer Aided Design, Vol. 10, No 4, pp. 501511, (1991).
 C. Riccobene, K. Gartner, G. Wachutka, H. Baltes and W. Fichtner, ‘First ThreeDimensional Numerical Analysis of Magnetic Vector Probe’, Proceeding of IEDM94, pp. 727730, (1994).
 Rodrigo RodriguezTorres, Robert Klima and Siegfried Selberherr,’ThreeDimensional Analysis of a MAGFET at 300 K and 77 K’, ESSDERC 2002, pp. 151154, (2002).
 M. B. Patil, ‘New Discretization Scheme for 2dimensional Semiconductor Device Simulation on Triangular Grid’, IEEE Transactions on CAD of Integrated Circuits and Systems, Vol. 17, No. 11, pp. 11601165, (1998).
 Yie He, ‘Discretization method for Current Continuity Equation Containing Magnetic Field’,IEEE Transactions on Electron devices, Vol. 39, No. 4, pp. 820824 (1994).
 E. Jovanovic, T. Pesic and D. Pantic,’3D Simulation of CrossShaped Hall Sensor and its Equivalent Circuit Model’, Proceedings 24th International Conference on Microelectronics, Vol. 1, pp. 235238, (2004).