Time-domain Coupled Full Maxwell- and Drift-Diffusion-Solver for Simulating Scanning Microwave Microscopy of Semiconductors

This paper presents a transient electromagnetic Maxwell solver (MS), a transient semiconductor Poisson-Drift-Diffusion (PDD) solver, and their numerical coupling. The proposed numerical solution schemes are based on Finite Element Method (FEM). Due to the solvers complexity it is important to carefully examine the obtained initial results. Therefore, a simple 1D pnjunction diode illuminated by an external electromagnetic plane wave is considered, as the stationary solutions of this structure such as the depletion width, built-in voltage, and carrier concentration distribution can be analytically obtained. The presented initially obtained transient results converge well to the analytic stationary solutions. The electromagnetic waves reflected from a diode with zero bias and 0.8 V bias structure reveal a small signal difference within a wide frequency range, which is an encouraging initial step towards more realistic simulations of scanning microwave microscopy structures and arrangements. The extension of the presented field formulations and numerical methods to 2D and 3D problems is straightforward.


INTRODUCTION
Scanning Microwave Microscopy (SMM) is a powerful tool to predict material properties and quality without physically deforming the sample. In addition, SMM enables characterization of emerging materials and semiconductor devices based on their transient response. There have been many successful sample analyses based on SMM [1,2] mainly for bulk materials and dielectrics. Electromagnetic models that could simulate such setups with high accuracy are also available [3].
Predicting semiconductor properties and analyzing their behavior under microwave illumination requires a multi-physics simulation tool capable of simulating both Maxwell and semiconductor equations. In this work, a vector element based transient electromagnetic Maxwell solver nodal element based transient Poisson-Drift-Diffusion semiconductor solver, and their numerical coupling are presented. Both solvers are implemented by using Finite Element Method (FEM) [4,5]. A numerical approach to deal with one unknown vector function of electric field (E) and three scalar unknown functions [electric scalar potential (ϕ), hole-(p) and electron concentration (n)] within the coupling scheme is presented in Section 2. Section 3 presents the results obtained from defined exemplary cases in 1D with a plane-wave electromagnetic illumination. Section 4, describes further steps planned to be completed by the time of the PIERS Conference and concludes the paper.

PHYSICAL MODEL AND NUMERICAL METHOD
Simulation of electromagnetic behavior of a semiconductor device requires the computation of the electric field and the charge carrier distribution within the device. The set of equations modeling both the physics can be put together as follows [5,7,8]: where ϕ is the electric potential, E is the divergence free component of the electric field, E t is the total electric field, n(p) is the electron (hole) concentration, R is the net recombination rate, J is the current density, N + d and N − a are the ionized charge densities due to doping, µ p,n is the mobility, D p,n is the diffusion coefficient, ε and µ are the electric permittivity and magnetic permeability of semiconductor material respectively and q is elementary charge.
Equation (1) is the wave equation of the electric field, (2) is the wave equation of the electric potential, (4) and (5) are continuity equations for carriers. These are supplemented by (6) and (7), which represent the overall charge transport taking both drift and diffusion of carriers into account. For the recombination rate Read-Hall-Shockley recombination [7] is assumed.
This set of equations can be divided into two sets, where the semiconductor physics and the general electromagnetic propagation are defined and solved individually; and through current and total electric field, the solutions of each solvers could be coupled in an iterative manner.

Transient Poisson-Drift-Diffusion (PDD) Semiconductor Physics Solver
In order to lay out the structure of the transient semiconductor solver a pnjunction diode is defined as shown in Fig. 1 with abrupt doping profiles, and the boundary value problem (BVP) is constructed for 1D case for simplicity. The stationary solutions of this problem such as the depletion width, built-in voltage and carrier concentration distributions are known and can also be analytically obtained [7].
Equations (8)-(11) model the driftdiffusion based physics within the diode, and all the unknowns are scaled by using the normalization factors given in Table 1 in order to avoid scaling issues that might arise in numerical analysis. The boundary value problem is completed using the boundary conditions for the unknowns on the metal contacts (at x = 0 and x = 2L): A variant of the Gummel's method [6] has been used to discretize the coupled nonlinear equation system (8) iteratively until a convergence is reached. The weak forms of the system can be obtained by using the method of weighted residuals [4,5]: The semiconductor domain is discretized by using scalar linear 1D elements and the Galerkin choice of the weighting functions is used [4]. The matrices describing the system obtained from FEM discretization are where N i is the 1D linear shape function originating at the node i of the 1-D mesh.
The above discretization has yielded the following system of equations for semiconductor domain: The time discretization is based on the well-known backward difference approach (BDF) [4], which is applied on the Equations (26) and (27) as follows: Equations (25), (28), and (29) form a coupled nonlinear equation system. In every time step, starting from the initial values of the charge distribution, the linear Poisson Equation (25) is solved to obtain the potential distribution. The carrier conservation equations are then solved using an iterative under-relaxation procedure: where α is the relaxation factor, until a convergence for the time step is reached and all three unknowns of the system, φ, p, and n are obtained.

Transient Maxwell Solver
Time domain wave equation with respect to the electric field (1) is also discretized using FEM with the absorbing boundary conditions (ABC) given by Equation (32) at each end of the domain. The excitation is provided by E 0 on the right hand side of (32). In 1D arrangement the y-polarization of the electric field (E =â y E y ) is chosen.
Similar to Section 2.1, the weak form of the wave equation including the boundary conditions for the system defined in Fig. 2 becomes: (33) The computational domain is discretized using FEM with linear 1D elements and the Galerkin choice of the weighting functions is employed again (w i = N i ).

2L
ϕ=0V Si, p-type , Si, n-type , The matrix contributions of each element are In addition, for time discretization, the BDF approach is applied, which has yielded the following final equation in matrix form: (37) Current densities obtained from (6) and (7) then can be inserted into Equation (1) by also using (3) to ensure the coupling between the physics. Then solving both transient models in an iterative manner will yield a converged state for all unknowns.

RESULTS
As a first validation step, a transient behavior of a 1D pnjunction is simulated using the in-house developed PDD-solver presented in Section 2.1. Fig. 3 summarizes the obtained results. As visible in Fig. 3 the transient potential building up and depletion region formation at the junction can be seen when the p-and n-doped side are brought into contact at t = 0 without any external bias or field applied. With a bias of 0.8 V which is greater than the turn-on voltage of the diode, the evolution of built-in potential and depletion region can be seen in Fig. 4.
The final solutions presented in Fig. 3 and 4 are identical to the stationary analytical solutions available in the literature [7,8]. This confirms the accuracy of the developed transient PDD-solver.
Having confirmed the simulations of simple pnjunction the coupled case with Maxwell simulator in Fig. 2 is attempted. A Gaussian burst signal with central frequency of 50 GHz is applied in time domain from the electromagnetic input port at x = −L and the reflected electric field from the diode is collected back. By taking Fourier transform of the time domain signals, it is possible to obtain scattering parameters from the semiconductor sample in the frequency domain. When a bias is applied on the diode, the effect of the change in the charge distribution reveals itself as a shift in the reflected power as in Fig. 5.
To ensure a high accuracy level of the results in frequency domain after performing the FFTanalysis of the time domain results, much longer time than the time presented in Fig. 5(a) is simulated. The time scale in Fig. 5 is contracted in order to present the Gaussian burst with visible details.
The applied external E-field is polarized in the y-direction due to nature of 1D simulation (in order to obtain a wave propagating through the diode in the x-direction), and the current flow caused by carrier movement on the semiconductor domain takes place in x-direction. This then limits the effect of the semiconductor on reflected E-field, and only observable change in the reflected power in 1D simulation occurs when a different bias on the diode is applied resulting a different carrier distribution.  In order to be able to observe all the change caused by local induced current and charge distribution, a higher dimension analysis and simulations would be required. In that case, the effect of frequency could also be seen as the charges could only follow the external field up to some certain frequency.

CONCLUSION
In this work, a multiphysics formulation for coupling Maxwell equation with drift-diffusion based semiconductor physics is presented. Strongly coupled nonlinear transient Poisson-Drift-Diffusion equations to model semiconductor physics have been numerically decoupled and iteratively solved in time domain. Based on formulation presented in Section 2, a 1D solver is implemented using the FEM discretization and BDF time discretization. A 1D pnjunction is tested using the Poisson-Drift-Diffusion solver with different forward bias and transient formation of depletion region is observed. The obtained results are identical to the stationary analytical solutions available in the literature, which confirms the accuracy of the transient PDD-solver.
The effect of external E-field on the semiconductor in 1D is also demonstrated, which enables SMM simulations of semiconductor materials. The conventional parameters obtained by SMM, the frequency response and scattering parameters, are obtained by Fourier transform of time domain solutions during post-processing. The 2D implementation of the coupled formulation will also show the full effect of semiconductor on the external field as in 2D it will be possible to accommodate external field, in-plane with induced current in semiconductor. A transient 2D and 3-D Poisson-Drift-Diffusion solver coupled with the corresponding Maxwell solver is presently under development and the obtained 2-D results will be presented at the PIERS conference.