A Box Method Discretisation for the Drift-Diffusion Equations in a Magnetic Field

Introduction

The problem of discretising the semiconductor drift-diffusion 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 non-uniform 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₀ 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 Scharfetter-Gummel 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^-1is 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 Scharfetter-Gummel discretisation. This forms the basis of the VR-discretisation method of Patil [5], in a non-magnetic 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 non-magnetic 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 Scharfetter-Gummel discretisation by magnetic terms. We extended this to the Box Method and reduce the discretisation to the usual Scharfetter-Gummel discretisation, but with two modifications. The first is that there is a prefactor equal to (1 + a² + b² + c²)^-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 quasi-Fermi level for electron current and the hole quasi-Fermi 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 n-channel 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 least-squares current based method, but the Silvaco method avoids these.

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 split-drain device, with a uniform donor density of 10¹⁵cm^-3, was simulated using Magnetic2D with a source drain applied bias of 0.5 Volts. From the log-log 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 n-type doping under the split drains, p-type doping under the source and a uniform n-type background doping. The ATLAS syntax for this N-n-P 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 split-drains. 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 n-type and p-type 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.

The next structure considered is a 3-dimensional Hall cross structure. This consists of an n-type Greek cross shape embedded in a p-type substrate. Current is passed between the source and drain in the y-direction, a magnetic field of 1 Tesla is applied in the z-direction, and the Hall Voltage is measured between the Probe contacts in the x-direction. It is usual to add either an isolated conductive top plate to the n-type active region, or to have a shallow p+ layer diffused into the top of the n-type 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¹⁸ cm^-3 and is 0.1 microns deep.

Current is passed between the anode and cathode, and the open-circuited 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 current-voltage 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.

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 two-dimensional 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 quasi-Fermi level. The first subscript is the row number and the second subscript is the column number. We notate these values as Δ₁ and Δ₂ 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₀, 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 quasi-Fermi 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 non-uniform, the value of expression 7 will not usually be zero because it becomes a weighted sum of quasi-Fermi level differences. Only the sum of the differences in quasi-Fermi 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 non-magnetic 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

1. L. Andor, H. P. Baltes, A. Nathan and H.G. Schmidt-Weinmar, ‘Numerical Modelling of Magnetic-Field-Sensitive Semiconductor Devices’, IEEE Transactions on Electron Devices, Vol. ED-32, No. 7, pp. 1224-1230, (1985).
2. W. Allegretto, A. Nathan, H. Baltes, ‘Numerical Analysis of Magnetic-Field-Sensitive Bipolar Devices, IEEE Transactions on Computer Aided Design, Vol. 10, No 4, pp. 501-511, (1991).
3. C. Riccobene, K. Gartner, G. Wachutka, H. Baltes and W. Fichtner, ‘First Three-Dimensional Numerical Analysis of Magnetic Vector Probe’, Proceeding of IEDM94, pp. 727-730, (1994).
4. Rodrigo Rodriguez-Torres, Robert Klima and Siegfried Selberherr,’Three-Dimensional Analysis of a MAGFET at 300 K and 77 K’, ESSDERC 2002, pp. 151-154, (2002).
5. M. B. Patil, ‘New Discretization Scheme for 2-dimensional Semiconductor Device Simulation on Triangular Grid’, IEEE Transactions on CAD of Integrated Circuits and Systems, Vol. 17, No. 11, pp. 1160-1165, (1998).
6. Yie He, ‘Discretization method for Current Continuity Equation Containing Magnetic Field’,IEEE Transactions on Electron devices, Vol. 39, No. 4, pp. 820-824 (1994).
7. E. Jovanovic, T. Pesic and D. Pantic,’3D Simulation of Cross-Shaped Hall Sensor and its Equivalent Circuit Model’, Proceedings 24th International Conference on Microelectronics, Vol. 1, pp. 235-238, (2004).

Download PDF version of this article