Learning Compact Physics-Aware Delayed Photocurrent Models Using Dynamic Mode Decomposition
LLEARNING COMPACT PHYSICS-AWARE DELAYEDPHOTOCURRENT MODELS USING DYNAMIC MODEDECOMPOSITION ∗ JOSHUA HANSON † , PAVEL BOCHEV ‡ , AND
BILIANA PASKALEVA § Abstract.
Radiation-induced photocurrent in semiconductor devices can be simulated usingcomplex physics-based models, which are accurate, but computationally expensive. This presents achallenge for implementing device characteristics in high-level circuit simulations where it is com-putationally infeasible to evaluate detailed models for multiple individual circuit elements. In thiswork we demonstrate a procedure for learning compact delayed photocurrent models that are effi-cient enough to implement in large-scale circuit simulations, but remain faithful to the underlyingphysics. Our approach utilizes Dynamic Mode Decomposition (DMD), a system identification tech-nique for learning reduced order discrete-time dynamical systems from time series data based onsingular value decomposition. To obtain physics-aware device models, we simulate the excess carrierdensity induced by radiation pulses by solving numerically the Ambipolar Diffusion Equation, thenuse the simulated internal state as training data for the DMD algorithm. Our results show thatthe significantly reduced order delayed photocurrent models obtained via this method accuratelyapproximate the dynamics of the internal excess carrier density – which can be used to calculatethe induced current at the device boundaries – while remaining compact enough to incorporate intolarger circuit simulations.
Key words.
Dynamic Mode Decomposition, Ambipolar Diffusion Equation, delayed photocur-rent, machine learning, data-driven compact models.
1. Introduction.
Ionizing radiation can affect operation of electronics in mul-tiple application contexts including space, terrestrial and manmade nuclear environ-ments. Specifically, fluence of ionizing radiation generates electron-hole pairs withinsemiconductor devices, resulting in excess currents flowing through the devices. Suchexcess current is not present during device operation in a normal environment and itsimpact on electronic systems can be catastrophic, ranging from instantaneous inter-ruptions in service, loss of stored memory, and even burnout of the entire system. Forexample spacecraft depend on electronic components that must perform reliably overmissions measured in years and decades and space radiation is a primary source ofdegradation, reliability issues, and potentially failure for these electronic components.Physics-based modeling and simulation of radiation effects on electronic systems invarious radiation environments can facilitate understanding of the mechanisms gov-erning the radiation response of the electronic materials, parts, and systems, and canbe used to devise ways to mitigate radiation effects, and create new materials anddevices that are resilient to radiation exposure.Thus, computational analysis of radiation effects on electronic systems has utilityranging from guiding the initial designs of systems, setting up the design of exper-iments, and final qualification. At a device level the excess carrier behavior can be ∗ Sandia National Laboratories is a multimission laboratory managed and operated by NationalTechnology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of HoneywellInternational, Inc., for the U.S. Department of Energys National Nuclear Security Administrationunder contract DE-NA-0003525. † University of Illinois at Urbana-Champaign, [email protected] ‡ Sandia National Laboratories, [email protected] § Sandia National Laboratories, [email protected] a r X i v : . [ phy s i c s . c o m p - ph ] A ug Compact Photocurrent Models Using DMD accurately modeled by the Drift-Diffusion equations (DDE)[15] given by ∇ · ( (cid:15) E ) = q ( p − n + C ) (1.1) ∂n∂t = ∇ · ( nµ n E + D n ∇ n ) − R + g (1.2) ∂p∂t = ∇ · ( pµ p E + D p ∇ p ) − R + g (1.3)where n and p are the concentrations of the electrons and holes, respectively, µ n , µ p , D p , and D n are the carrier mobilities and diffusivities, C is the doping concentration, E is the electric field, q is the electron charge, R is the recombination rate, and g isthe generation rate.Numerical solution of (1.1)–(1.3) using, e.g., a multi-dimensional finite elementdiscretization of the device, forms the basis of the so-called Technology Computer-Aided Design (TCAD) simulators. However, (1.1)–(1.3) is a coupled system of non-linear Partial Differential Equations (PDEs) whose accurate numerical solution canbe very time consuming. As a result, although TCAD device simulators based onDDE could in principle be coupled to circuit simulators, their high computationalcost makes the analysis of all but very small circuits computationally intractable.At the other end of the spectrum are the so-called compact device models, whichare computationally efficient but rely on empirical approximations and/or simplifiedanalytic solutions to the semiconductor transport equations. Often such models mustbe recalibrated for different operating regimes and inclusion of new physics may forceredevelopment of the model from scratch.Data-driven techniques present an opportunity to automate the development ofcompact models by learning them directly from laboratory measurements and/orsuitable synthetic data. For such models to be practically useful though, they must beable to correctly predict the device response when the device is integrated into a circuitand exposed to a wide range of stimuli. The caveat is that in a laboratory settingsone can only apply a limited type of signals and directly measure the device response,leading to “sparse” training sets. In contrast, traditional Machine Learning (ML)applications, such as natural language processing and image classification [10, 9, 7]operate in “big data” environments. As a result, a compact device model learnedsolely from such “sparse” data may fail to generalize to all relevant analysis conditionsbecause it will learn salient patterns in the dataset rather than the causal physicsunderpinning the device operation.The latter, i.e., the fact that device behavior is governed by strict physics lawscan be used to counter the lack of big data by incorporating the physics knowledgeinto the model development. Recent work on scientific Machine Learning [13, 12, 2]suggests this strategy is effective and can lead to generalizable models.The main goal of this work is to establish the viability of physics-aware machinelearning for the development of compact data-driven photocurrent models that areefficient enough to allow large circuit simulations, while remaining faithful to thebasic physics principles embodied by models such as (1.1)–(1.3). To that end weadopt a setting that has been used extensively in the past five decades at SandiaNational Laboratories for the development of compact analytic photocurrent models.Although this setting uses a simplified version of the DDE equations, it provides astepping stone towards future developments of data-driven models based on the fullycoupled system (1.1)–(1.3). In so doing we are able to leverage a wealth of experiencesand physics knowledge accumulated through the use of the existing compact models, . Hanson, P. Bochev, and B. Paskaleva N D P ℓ n W ℓ p Light or radiation
Fig. 2.1 . Traditional compact photocurrent model development splits each
P N -junction within adevice into a depletion region D having width W , and quasi-neutral P and N regions with lengths (cid:96) n and (cid:96) p , respectively; see, e.g., [1].
2. Physics-Based Compact Analytic Photocurrent Models.
Althoughthe DDE system (1.1)-(1.3) accurately describes the behavior of the excess carriers,it does not lend itself to exact analytical solution. To enable analytic approximationof the governing equations, most photocurrent models in use today follow the samebasic approach as in the classic paper [17] and split each
P N -junction within a deviceinto a depletion region and quasi-neutral P - and N -regions; see Fig. 2.1. Carriersin the depletion region are quickly converted to photocurrent and yield the so-calledprompt photocurrent I prompt . Carriers in the P and N regions have a delayed re-sponse and produce the delayed photocurrents I p and I n , respectively. As a result,the total photocurrent is given by I total = I prompt + I p + I n . To calculate these currents one makes additional simplifying assumptions. The first isthe electrical neutrality approximation, which stipulates that the excess electron andhole densities are equal throughout the entire device. The second is the congruenceassumption, which states that the electron and hole fluxes into or out of any regionmust be equal; see, e.g., [8]. Under these assumptions the DDE model (1.1)-(1.3) inthe P and N regions can be replaced by the Ambipolar Diffusion Equation (ADE) ∂u∂t = D a ∇ u − µ a E · ∇ u − τ a u + g (2.1)where u , D a , µ a and τ a are the excess carrier density (electrons or holes), the ambipo-lar diffusion coefficient, the ambipolar mobility and the carrier lifetime, respectively.In general, these parameters may depend on the excess carrier density and, as a re-sult, (2.1) is still a nonlinear PDE. However, for moderate radiation dose rates thesecoefficients can be approximated by constant values and the ADE becomes a linearparabolic PDE. The final assumption is that the depletion region width W is not af-fected by the excess carriers and is a constant. Under these assumptions, the prompt Compact Photocurrent Models Using DMD photocurrent is modeled as I prompt = qgAW , where A is the effective area of the P N -junction, while I p and I n are modeled by theADE (2.1).At this point traditional compact photocurrent model development proceeds withderiving analytic approximations for the solutions of (2.1) in one dimension and usingthem to obtain expressions for the delayed photocurrents. Early work [17] consideredunbounded P and N regions and negligible electric fields. The resulting Wirth-Rogersmodel tends to overestimate the photocurrent as it neglects the effects of an ohmiccontact at a finite distance from the depletion region [1]. Subsequent work [4] relaxedthese conditions, assumed that E is constant, and used approximate Laplace trans-forms to obtain analytic expressions for I p and I n . However, approximation of theLaplace transform results in a model that yields unphysical current estimates when E exceeds roughly 10 V cm − ; see [18]. The latter work used Fourier analysis tech-niques to develop a more accurate photocurrent model that avoids these drawbacks.The model in [18] was further improved in [1] by using a transformation of (2.1) intoan inhomogeneous heat equation and solving the latter exactly by Fourier techniques.Another popular photocurrent model is the Fjeldly Model [5]. Unlike the Axness-Kerrmodel [1], which solves the time-dependent equation (2.1) exactly, the Fjeldly Modeluses a steady-state solution of ADE combined with an RC delay circuit to achievetime dependence.In this work we adopt the above setting and focus on the development of physics-aware compact data-driven models for the delayed photocurrents in the P - and N -regions that can be used as ”plug-and-play” substitutes for conventional models. Thecore idea is to replace the simplified analytic solutions of the ADE, comprising thebasis of most standard compact models, by a discrete-time dynamical system approxi-mating the flow map of the ADE, i.e. the “solution operator” that returns the internaldevice state for any given external input.To that end we apply a Dynamic Mode Decomposition (DMD) [14] approachto samples of the internal state (carrier density) of the device obtained by solving(2.1) numerically by a Finite Element Method (FEM). In this paper we assume thatall device parameters such as diffusion coefficients, carrier lifetimes and doping con-centrations are known. In this case ADE already contains all the necessary physicsinformation and development of physics-aware compact data driven models can bedone entirely from synthetic data .A more general setting occurs when one or more of the ADE parameters areeither unknown or have large uncertainties. In this case the development of thecompact model must also include a parameter identification step to refine the ADE,which requires laboratory data. Since the purpose of this paper is to demonstrate theviability of DMD as an effective tool for the generation of compact device models,detailed discussion of this more general setting is beyond the scope of this paper andwill be addressed in a forthcoming work.Regardless of the particular setting though the physics-awareness of our modelsstems from the fact that they represent approximations of the flow map engenderedby the physics model , and built from simulated internal states of this model. Thesestates contain physics information that cannot be obtained by laboratory instruments, In this setting the resulting compact model can be interpreted as a discrete reduced order modelof the flow map generated by the ADE.. Hanson, P. Bochev, and B. Paskaleva
3. Numerical solution of the ADE.
We consider the ADE on the space-timedomain Ω :=
X × T ⊂ R , where X := (0 , (cid:96) ) and T := (0 , t final ). Without loss ofgenerality, one may assume that X is the N -region of the device; see Fig. 2.1 andso we set (cid:96) = (cid:96) n . To obtain a well-posed problem, we augment equation (2.1) withhomogeneous initial and boundary conditions, i.e., u ( x,
0) = 0 ∀ x ∈ X and u (0 , t ) = u ( (cid:96), t ) = 0 ∀ t ∈ T . The homogenous initial condition corresponds to the fact that at t = 0 there are noexcess carriers present in the device. The boundary condition choice corresponds toassuming infinite carrier recombination velocities and ohmic contacts at x = 0 andthe boundary of the depletion region x = (cid:96) ; see [1]. More general boundary conditionsnot requiring these assumptions can also be considered but are not necessary for thepurpose of this work.We will simulate the internal state of the device by using the method of linesto solve the ADE numerically. For the spatial discretization we consider a stan-dard Galerkin finite element method and then solve the resulting system of OrdinaryDifferential Equations (ODEs) using an implicit numerical integration scheme. Forcompleteness we briefly review the discretization process below.Let X h denote a uniform partition of X into n + 1 elements κ i with vertices { x i } n +1 i =0 , i.e., κ i = [ x i , x i +1 ], i = 0 , . . ., n , x = 0 and x n +1 = (cid:96) . The mesh parameteris given by h = (cid:96)n +1 .As usual, L ( X ) denotes the space of all square integrable functions on X withnorm and inner product denoted by (cid:107) · (cid:107) and ( · , · ) , respectively, and H ( X ) is theSobolev space of order one whose elements are constrained to vanish at the boundarypoints. The weak variational form of the ADE is then given by seek u ∈ H ( X ) suchthat ( u t , v ) + Q ( u, v ) = ( g, v ) ∀ v ∈ H ( X ) . (3.1)The bilinear form Q ( · , · ) : H ( X ) × H ( X ) (cid:55)→ R is defined as Q ( u, v ) = D a ( u x , v x ) + µ a E ( u x , v ) + 1 τ a ( u, v ) . (3.2)To discretize (3.1) in space we consider a nodal (Lagrangian) conforming finite elementsubspace V h ⊂ H ( X ); see, e.g., [3]. Let { v i } i =1 ,n be the standard nodal basis havingthe property that v i ( x j ) = δ ij . We then seek an approximate solution of (2.1) as u h ( x, t ) = n (cid:88) i =1 u i ( t ) v i ( x ) , (3.3)where u ( t ) = ( u ( t ) , . . . , u n ( t )) ∈ R n is a vector of unknown solution coefficients.Inserting (3.3) into the weak form (3.1) and restricting the test space to V h thenyields the system of ODEs M ˙ u ( t ) + K u ( t ) = g (3.4) Utilizing a variable mesh spacing with increased node density near the boundary may providesome potential advantages such as more accurate gradient estimation at the edge points, however inthe context of this work we restrict our attention to uniform meshes.
Compact Photocurrent Models Using DMD where
M, K ∈ R n × n are the (consistent) finite element mass and stiffness matriceswith elements M ij = ( v i , v j ) and K ij = Q ( v i , v j ) , respectively, and g ( t ) ∈ R n is a discrete source term with g i ( t ) = ( g, v i ) .The ODE system (3.4) can be solved by any standard time-integration scheme.However, in general (3.4) is stiff and an implicit scheme is preferred. In this paper, weuse an implicit multi-step variable-order routine based on a backward differentitationformula (BDF) for approximating the state derivative; for more details see [16]. Thismethod is included by default in the scipy.integrate submodule within the SciPyv1.5.1 package for Python 3.
4. A Dynamic Mode Decomposition Compact Photocurrent Model.
Assume that an approximate solution of the ODE (3.4) is available at uniformlyspaced time steps t k := k ∆ t , k = 0 , . . . , m , for some sampling interval ∆ t > u and the input g , i.e., u k +1 = A u k + B g k (4.1)where u k := [ u ( k ∆ t ) · · · u n ( k ∆ t )] T (4.2) g k := g ( k ∆ t ) (4.3)and A, B ∈ R n × n are linear maps. Expressing the system in this form facilitates theuse of many familiar system identification techniques. In particular, Dynamic ModeDecomposition (DMD) is a data-driven method for learning the maps A and B from atime series of state and input measurements { u k , g k } mk =0 , with the goal of identifying asmall number of dominant dynamic modes. Below we summarize the DMD algorithmadapted for control inputs as described in [11]. By organizing the samples into thefollowing matrices X = | | u · · · u m − | | ; X (cid:48) = | | u · · · u m | | ; and G = | | g · · · g m − | | , (4.4)we can express the linear relationships within the data as X (cid:48) = AX + BG = (cid:2) A B (cid:3) (cid:20) XG (cid:21) =: (cid:2) A B (cid:3) S. (4.5)Therefore the maps A and B can be approximated by (cid:2) A B (cid:3) ≈ (cid:2) ¯ A ¯ B (cid:3) := X (cid:48) S † (4.6)where † indicates the Moore-Penrose pseudoinverse. An efficient and accurate algo-rithm for estimating the pseudoinverse of a rectangular matrix is realized via truncatedsingular value decomposition. The matrix of samples S can be factored as S = U Σ V T = (cid:2) ˜ U ˜ U trun (cid:3) (cid:20) ˜Σ 00 ˜Σ trun (cid:21) (cid:20) ˜ V T ˜ V T trun (cid:21) ≈ ˜ U ˜Σ ˜ V T (4.7) . Hanson, P. Bochev, and B. Paskaleva U ∈ R n × n , ˜ U ∈ R n × p , Σ ∈ R n × m , ˜Σ ∈ R p × p , V T ∈ R m × m , and ˜ V T ∈ R p × m .Here trun denotes the m − p truncated singular values, so that the p greatest singularvalues are kept. In practice, one sets an error tolerance and truncates all singularvalues below this threshold. Now the pseudoinverse of S is naturally approximatedby S † = ˜ V ˜Σ − ˜ U T (4.8)Substituting equation (4.8) into (4.6) gives A ≈ ¯ A ≈ X (cid:48) ˜ V ˜Σ − ˜ U T ∈ R n × n (4.9) B ≈ ¯ B ≈ X (cid:48) ˜ V ˜Σ − ˜ U T ∈ R n × n , (4.10)where ˜ U = (cid:2) ˜ U T ˜ U T (cid:3) T with ˜ U ∈ R n × p , ˜ U ∈ R n × p .We can achieve a more compact model by incorporating the projection ˜ u = P u of the state onto the canonical dynamic mode coordinates. In the same manner asequation (4.7), we factor the matrix X (cid:48) ≈ ˆ U ˆΣ ˆ V T , where the truncations are chosento preserve the r greatest singular values. Then the projection onto dynamic modecoordinates is given simply by P = ˆ U T ∈ R r × n , hence the desired reduced ordermodel obtained via DMD is ˜ u k +1 = ˜ A ˜ u k + ˜ B g k (4.11)˜ u = ˆ U T u (4.12) u k ≈ ˆ U ˜ u k (4.13)where ˜ A := ˆ U T ¯ A ˆ U = ˆ U T X (cid:48) ˜ V ˜Σ − ˜ U T ˆ U ∈ R r × r (4.14)˜ B := ˆ U T ¯ B = ˆ U T X (cid:48) ˜ V ˜Σ − ˜ U T ∈ R r × n . (4.15)The parameter p represents the number of dynamic modes to fit to the data,which controls the model precision. The parameter r represents the number of modesto project onto, that is, the order of the final reduced-order model, which controls themodel compactness. The case where r > p usually results in diminished performance; r = p retains exactly the same number of modes fit to the data in the compactifiedmodel; r < p results in a more compact model, but ignores the p − r least significantmodes fit to the data, which can result in slightly reduced accuracy. In this work wewill use r = p . Construction of the DMD model (4.11)requires training samples representing time series of state and input measurements { u k , g k } mk =0 . To generate such samples we recall the assumption in Section 2 that alldevice parameters are known. In particular, here we consider a generic P N -junctiondevice characterized by the parameters in Table 4.1. The values in this table areadapted from [1] to enable a direct comparison with a published compact photocurrentmodel.We then select a suitable set of generation density functions { g k train } Mk =1 and usethe computational scheme in Section 3 to solve (2.1) numerically with homogeneousinitial and Dirichlet boundary conditions. Selection of the inputs g k train depends on Compact Photocurrent Models Using DMD
Table 4.1
ADE parameters for a generic
P N -junction device
Parameter Value Units Description (cid:96) . × − cm N -region length D a . × cm s − diffusion coefficient µ a . × cm V − s − typical hole mobility for Si(cid:96) n . × − cm s − diffusion length: √ D a τ a τ a . × − s typical hole lifetime E −
20 or 0 V cm − electric field C × cm − doping concentration (cid:98) g . × cm − s − max generation density u ( x, t ) variable cm − excess carrier concentration (ADE solution) the type of the anticipated testing input(s) g test for the model and will be revisitedin Section 5.Obtaining a numerical solution for the ADE requires proper scaling and non-dimensionalization of the governing equations. For convenience we scale the compu-tational domain X so that the length of the N -region, i.e., (cid:96) , becomes a unit of lengthand t final becomes a unit of time. After rescaling the domain and the equations thecomputational domain X becomes the unit square, and the non-dimensional ADEcoefficients are given by D a = 0 . µ a = 0 .
063 ; L p = 0 . τ a = 3 .
92 ; and E = − .
84 or 0 (cid:98) g = 2 . . (4.16)We highlight that the numerical solution of the ADE constitutes the ”physics-based” element of our procedure. Specifically, the physics information is incorporatedby using trusted a priori dynamics models – which are calibrated or driven by exper-imental measurements – to generate training data from simulating the unobservableinternal state of a device using robust numerical techniques. Remark 1.
In this work we use the ADE as the physics basis for the data-drivenmodel because it has been used to develop almost all compact photocurrent models inuse today, i.e., it is an example of model that is trusted based on decades-long practicalexperiences. However, we emphasize that the DMD algorithm is also suitable for morecomplex physics-based models, such as the full drift-diffusion equations or detailedmolecular dynamics simulations, with the capability of producing dramatically reducedorder approximations that are feasible to implement in high-level circuit simulatorsbut remain faithful to the underlying physics.
5. Simulation Results.
Assuming that a DMD model (4.11) has been trainedaccording to the procedure in Section 4.1, we test its performance as follows. Let g ( x, t ) be a target generation density for which we seek the response of our device.We sample g ( x, t ) in space using the vertices defining the finite element mesh X h ,and in time using a desired time step ∆ t for a total of m time steps. This samplingproduces the inputs g k to the DMD model. We then set u = and use (4.11) torecover the internal state of the device, i.e., the excess carrier concentration, at themesh nodes { x i } ni =1 for every t k = k ∆ t :˜ u k +1 = ˜ A ˜ u k + ˜ B g k , k = 0 , . . . , m − . Each vector ˜ u k induces a C finite element function u DMD h ( x, t k ) = n (cid:88) i =1 ˜ u i ( t k ) v i ( x ) (5.1) . Hanson, P. Bochev, and B. Paskaleva u DMD h and taking into accountthe homogeneous Dirichlet boundary condition u DMD h ( x , t k ) = u DMD h ( x n +1 , t k ) = 0,we define the approximate boundary photocurrent at t = t k as J DMD h ( x, t k ) = D a ∂ x u DMD h ( x, t k ) = D a n (cid:88) i =1 ˜ u i ( t k ) ∂ x v i ( x ) , (5.2)for x = x and x = x n +1 . Note that owing to the local support of the basis functions v i ( x ), J DMD h ( x , t k ) and J DMD h ( x , t k ) only include contributions from the basis func-tions supported on elements κ and κ n , respectively. We then compare the predictedDMD flux with simulated experimental measurements of the current out of the deviceterminals. These measurements are obtained using the same numerical procedure asin Section 3, i.e., by computing a finite element solution u FEM h ( x, t k ) and setting J FEM h ( x , t k ) = D a ∂ x u FEM h ( x , t k ) (5.3) J FEM h ( x n +1 , t k ) = D a ∂ x u FEM h ( x n +1 , t k ) . (5.4) We first evaluate the ability of the DMDmodel to reproduce an artificially manufactured solution: u MNF ( x, t ) := te − t sin (cid:0) πx(cid:96) (cid:1) . (5.5)Observe that this solution satisfies homogenous initial and boundary conditions. Since(5.5) is expressed in closed-form, we can directly compute the boundary photocurrentof the manufactured solution: J MNF ( x, t ) = D a ∂ x u MNF h ( x, t ) = D a te − t (cid:0) π(cid:96) (cid:1) cos (cid:0) πx(cid:96) (cid:1) . (5.6)Thus, for this test we will compare J DMD with the known manufactured solution cur-rent J MNF , rather than with the simulated current J FEM in (5.3)–(5.4). Substitutionof u MNF ( x, t ) into the governing equation (2.1) yields a generation density g MNF ( x, t ) = ∂u MNF ∂t ( x, t ) − D a ∂ u MNF ∂x ( x, t ) + µ a E ( x ) ∂u MNF ∂x ( x, t ) + 1 τ a u MNF ( x, t )(5.7)= (1 − t ) e − t sin (cid:0) πx(cid:96) (cid:1) + D a te − t (cid:0) π(cid:96) (cid:1) sin (cid:0) πx(cid:96) (cid:1) (5.8)+ µ a E ( x ) te − t (cid:0) π(cid:96) (cid:1) cos (cid:0) πx(cid:96) (cid:1) + 1 τ a te − t sin (cid:0) πx(cid:96) (cid:1) (5.9)such that when the ADE (2.1) is driven by g MNF ( x, t ), its solution will exactly matchthe desired manufactured solution u MNF ( x, t ). The inputs g k to the DMD model areobtained by sampling g MNF ( x, t ) according to the method described earlier.Since the input g MNF ( x, t ) corresponding to the manufactured solution is spatiallyirregular, a spatially uniform training input will generally result in poor performance.To address this, we design a sequence of localized pulses which will excite differentregions of the device. This sequence will then be used as the training input for theDMD model. We choose to use a Gaussian profile that has been windowed by a cosine0 Compact Photocurrent Models Using DMD function as the spatial envelope for the input pulses, where the window function isapplied to restrict the support of the envelope to a compact interval. This profile isconsistent with experimentally viable radiation doses; other reasonable choices include
Lorenz or Voigt profiles, which reflect different radiation broadening mechanisms. Thewindowed Gaussian profile with center x i and support [ x i − w , x i + w ] is given by ρ i ( x ) := (cid:40) cos (cid:0) π x − x i w (cid:1) exp (cid:0) − (cid:0) x − x i w (cid:1) (cid:1) if x i − w ≤ x ≤ x i + w x < x i − w or x > x i + w (5.10)where x i = i (cid:96)N pulses − for i = 0 , . . . , N pulses −
1, and w = (cid:96)N pulses − . The profiles ρ i ( x )are illustrated in Figure 5.1. Combining the spatial envelopes defined above with asquare temporal envelope gives g i ( x, t ) = (cid:40)(cid:98) gρ i ( x ) if 0 ≤ t ≤ . µ s0 otherwise , (5.11)where the value of (cid:98) g is defined in (4.16). Now we can express the training input as g train = N pulses − (cid:88) i =0 g i ( x, t − t i ) (5.12)where t i = i (5 . µ s) for i = 0 , . . . , N pulses −
1. For the manufactured solution (5.5),we choose N pulses = 10.To generate the training samples, we solve (2.1) with the source (5.12) and theparameters (4.16) on a mesh X h comprising 512 uniform elements, and sample thesolution in time at ∆ t = 0 . µ s increments. Training input spatial envelopes .
000 0 .
005 0 .
010 0 .
015 0 .
020 0 .
025 0 . Position x [ µ m ] . . . . . . I npu t en v e l ope ρ i ( x ) [ ∼ ] Fig. 5.1 . Generation density function g train used to obtain the training set, and g test used to verifythe model performance. The square edges in the training input are included to excite a wider range ofdynamic modes, which is due to the high-bandwidth content in the sharp transitions. Using the the manufactured photocurrent (5.6) and carrier density (5.5) andapproximate photocurrent (5.2) and approximate carrier density (5.1) derived fromthe DMD model, we define the following error quantities: E J ( t ) := | J DMD h ( t ) − J MNF ( t ) | (5.13) E u ( t ) := | u DMD h ( t ) − u MNF ( t ) | (5.14) . Hanson, P. Bochev, and B. Paskaleva Manufactured solution flux(no electric field) . . . . N o r m a li z ed fl u x J ( t ) [ A / c m ] Left boundary − . − . − . . Right boundary p=13p=14p=15p=16p=17FEMMNF
Time t [ µ s ] − − − N o r m a li z ed fl u x e rr o r E J ( t ) [ A / c m ] Time t [ µ s ] − − − p=13p=14p=15p=16p=17FEM Fig. 5.2 . Normalized photocurrent due tothe manufactured input (5.7) with no electric fieldapplied. The top two plots show the manufac-tured solution (5.6) , FEM solution (5.3) - (5.4) , andDMD solution (5.2) , and the bottom two plotsshow the DMD and FEM error (5.13) . The lefttwo plots correspond to the flux out of the left sideof the N -region, and similarly on the right. Manufactured solution flux(with electric field) . . . . N o r m a li z ed fl u x J ( t ) [ A / c m ] Left boundary − . − . − . . Right boundary p=13p=14p=15p=16p=17FEMMNF
Time t [ µ s ] − − − N o r m a li z ed fl u x e rr o r E J ( t ) [ A / c m ] Time t [ µ s ] − − − p=13p=14p=15p=16p=17FEM Fig. 5.3 . Normalized photocurrent due tothe manufactured input (5.7) with an electric field (4.16) applied. The top two plots show the manu-factured solution (5.6) , FEM solution (5.3) - (5.4) ,and DMD solution (5.2) , and the bottom two plotsshow the DMD and FEM error (5.13) . The lefttwo plots correspond to the flux out of the left sideof the N -region, and similarly on the right. which will characterize the ability of the DMD model to reproduce manually selecteddynamic modes. For simplicity, in this work the electric field is not considered aninput to the system, so the DMD model must be trained separately for each uniqueelectric field strength.We report results for the DMD model using several different values for the param-eter p , which defines the number of the greatest singular values kept in the truncatedSVD decomposition (4.8). Figure 5.3 compares the boundary photocurrent from themanufactured solution and from the FEM and DMD models subject to the manufac-tured input (5.7). Figures 5.5 and 5.7 show the carrier density from the manufacturedsolution and from the FEM and DMD models, with respect to time t or position x ,respectively. Figures 5.2, 5.4, and 5.6 show the same, but in the absence of an electricfield. In this test we compare our compact DMD model withthe Axness-Kerr [1] compact analytic model. We consider the case of a lightly dopeddiode as described in [1, Section B, p.2650] for which the length of the N -region equals2 diffusion lengths (case ξ p = 2 in [1, Figure 3, p.2651].) The device is irradiated bya 1 . µ s step pulse g test ( x, t ) = (cid:40)(cid:98) g if 0 ≤ t ≤ . µ s0 otherwise (5.15)with the value of (cid:98) g rescaled as in (4.16). This example from [1] corresponds to theparameters in 4.16 with the exception of a few corrections to account for typographicalerrors in that paper.In contrast to the manufactured solution test where the desired input g MNF ( x, t )to the model is spatially irregular, now the target generation density, defined in (5.15),is spatially constant. As a result, the training input does not have to excite different2 Compact Photocurrent Models Using DMD
Manufactured solution density (noelectric field) − . . . C a rr i e r den s i t y u ( x j , t ) [ / c m ] x = 0 . ℓ x = ℓ x = 1 . ℓ p=13p=14p=15p=16p=17FEMMNF − − − − C a rr i e r den s i t y e rr o r E u ( x j , t ) [ / c m ] Time t [ µ s ] − − − − − − − − p=13p=14p=15p=16p=17FEM Fig. 5.4 . Simulated excess carrier densitydue to the manufactured input (5.7) with no elec-tric field applied. The top three plots show themanufactured solution (5.5) , FEM solution (3.3) ,and DMD solution (5.1) , and the bottom threeplots show the DMD and FEM error (5.14) . Manufactured solution density(with electric field) − . . . C a rr i e r den s i t y u ( x j , t ) [ / c m ] x = 0 . ℓ x = ℓ x = 1 . ℓ p=13p=14p=15p=16p=17FEMMNF − − − − C a rr i e r den s i t y e rr o r E u ( x j , t ) [ / c m ] Time t [ µ s ] − − − − − − − − p=13p=14p=15p=16p=17FEM Fig. 5.5 . Simulated excess carrier densitydue to the manufactured input (5.7) with an elec-tric field (4.16) applied. The top three plots showthe manufactured solution (5.6) , FEM solution (3.3) , and DMD solution (5.1) , and the bottomthree plots show the DMD and FEM error (5.14) . Manufactured solution densitysnapshot (no electric field) .
000 0 .
008 0 .
016 0 .
024 0 . − . . . C a rr i e r den s i t y u ( x , . µ s ) [ / c m ] p=13p=14p=15p=16p=17FEMMNF .
000 0 .
008 0 .
016 0 .
024 0 . Position x [ µ m ] − − − − C a rr i e r den s i t y e rr o r E u ( x , . µ s ) [ / c m ] p=13p=14p=15p=16p=17FEM Fig. 5.6 . Snapshot of the manufactured ex-cess carrier density at t = 5 . µ s due to the man-ufactured input (5.7) with no electric field ap-plied. The top plot shows the manufactured so-lution (5.6) , FEM solution (3.3) , and DMD so-lution (5.1) , and the bottom plot shows the DMDand FEM error (5.14) . Manufactured solution densitysnapshot (with electric field) .
000 0 .
008 0 .
016 0 .
024 0 . − . . . C a rr i e r den s i t y u ( x , . µ s ) [ / c m ] p=13p=14p=15p=16p=17FEMMNF .
000 0 .
008 0 .
016 0 .
024 0 . Position x [ µ m ] − − − − C a rr i e r den s i t y e rr o r E u ( x , . µ s ) [ / c m ] p=13p=14p=15p=16p=17FEM Fig. 5.7 . Snapshot of the manufactured ex-cess carrier density at t = 5 . µ s due to the man-ufactured input (5.7) with an electric field (4.16) applied. The top plot shows the manufactured so-lution (5.6) , FEM solution (3.3) , and DMD so-lution (5.1) , and the bottom plot shows the DMDand FEM error (5.14) . regions of the device and can be chosen to be spatially constant as well. Thus, forthis example we choose the training input g train ( x, t ) to be a constant in space and a . Hanson, P. Bochev, and B. Paskaleva Training input and test input functions . . . . . . Time t [ µ s ] . . . . . G ene r a t i onden s i t y g ( ℓ , t ) g train ( ℓ,t ) g test ( ℓ,t ) Fig. 5.8 . Generation density function g train used to obtain the training set, and g test used to verifythe model performance. The square edges in the training input are included to excite a wider range ofdynamic modes, which is due to the high-bandwidth content in the sharp transitions. discontinuous piecewise linear in time function such that g train ( x, t ) = t < µ s or 2 ≤ t < µ s or 4 t ≥ µ s (cid:98) g (4 t ) if 0 ≤ t < µ s (cid:98) g if 1 ≤ t < µ s or 5 ≤ t < µ s (cid:98) g (7 − t ) if 6 ≤ t < µ s , (5.16)See Figure 5.8 for an illustration of g train and g test .We then solve (2.1) with the source (5.16) and the parameters (4.16) on a mesh X h comprising 1024 uniform elements, and sample the solution in time at ∆ t = 0 . µ sincrements. The compact DMD photocurrent model for the device is now definedaccording to (4.11). Remark 2.
Since the input is applied uniformly across the entire device (andthus the entire state space), we could reduce the input dimension of the DMD model (4.11) from N to , however the higher dimensional input permits spatially irregu-lar excitations such as localized radiation pulses or non-transversal plane waves asrequired for the manufactured solution test in Section 5.1. To determine an appropriate dimension for the reduced-order DMD model, it isinformative to inspect the relative magnitudes of the singular values from the decom-positions of the sample matrices. For the training input (5.16), Figure 5.9 illustratesthe magnitude roll-off in the singular values for the state and state-input sample ma-trix decompositions. Observe that even for tight error thresholds, only a few modesare necessary to construct an accurate approximation of the state transition and inputmatrices.We wish to emphasize that including too many modes in the reduced-order modeloften leads to a realization which is unstable over a long time horizon. This phe-nomenon occurs due to unstable modes corresponding to small singular values. In-cluding these low-magnitude modes will typically lead to a better fit for the statedynamics resulting from the training input over the original training time horizon,but may lead to model divergence over a longer time horizon, signifying an exampleof overfitting .4 Compact Photocurrent Models Using DMD
Singular value decay
Index i − − − − − M agn i t ude σ i ˜ΣˆΣ Singular value decay
Fig. 5.9 . When listed from greatest to least, the singular values in the decomposition of the state-input sample matrix S ( ˜Σ ) and state sample matrix X (cid:48) ( ˆΣ ) demonstrate a nearly exponential decay inmagnitude. Based on the photocurrent (5.3)-(5.4) and carrier density (3.3) derived from thefinite-elements model and the approximate photocurrent (5.2) and approximate carrierdensity (5.1) derived from the DMD model, we redefine the following error quantities: E J ( t ) := | J DMD h ( t ) − J FEM h ( t ) | (5.17) E u ( t ) := | u DMD h ( t ) − u FEM h ( t ) | (5.18)which will facilitate the performance evaluation of the DMD model for typical, spa-tially uniform input functions. As before, the DMD model must be trained separatelyfor each electric field strength.Again, we show the results for the DMD model using multiple values for theparameter p . However, observe that for the case where the generation density isspatially uniform, fewer modes are necessary to achieve good performance than whenthe input magnitude varies along the length of the device, as in the manufacturedsolution test. Figures 5.10 and 5.11 compare the boundary photocurrent producedby the FEM and DMD models for the training input (5.16) and test input (5.15),respectively. Figures 5.12 and 5.13 show the same, but in the absence of an electricfield. Figures 5.14 and 5.15 illustrate the simulated excess carrier density (i.e., theinternal state) from the FEM and DMD models for both the training and test inputs.Several conclusions can be drawn from these results. First, the photocurrentplots at the left boundary of the N -region shown in Figures 5.11 and 5.13 are inan excellent agreement with the results reported in [1, Figure 3, p.2651] for ξ p = 2.Second, the error plots on the bottom rows of the figures quantify the differencesbetween the reference FEM solution and its DMD approximation as a function ofthe number p of selected dynamical modes. These results reveal that, as expected,the error decreases with increase of the number of dynamic modes; however, evenwith just 6 modes selected the DMD photocurrent model yields excellent accuracy.Overall, these results suggest that a data-driven approach is indeed a viable andeffective alternative to traditional analytic model development that can be used toquickly develop accurate and computationally efficient photocurrent models directlyfrom data.
6. Conclusions.
We have developed a compact data-driven delayed photocur-rent model given by a low-dimensional discrete-time dynamical system, which approxi- . Hanson, P. Bochev, and B. Paskaleva FEM vs. DMD flux(training input, with electric field) . . . . N o r m a li z ed fl u x J ( t ) [ A / c m ] Left boundary − . − . − . . Right boundary p=6p=7p=8p=9p=10FEMInput . . . . . Time t [ µ s ] . . . . N o r m a li z ed fl u x e rr o r E J ( t ) [ A / c m ] . . . . . Time t [ µ s ] . . . . p=6p=7p=8p=9p=10Input Fig. 5.10 . Normalized photocurrent due tothe training input (5.16) with an electric field (4.16) applied. The top two plots show the FEMsolution (5.3) - (5.4) and DMD solution (5.2) , andthe bottom two plots show the DMD training error (5.17) . The left two plots correspond to the fluxout of the left side of the N -region, and similarlyon the right. FEM vs. DMD flux(test input, with electric field) N o r m a li z ed fl u x J ( t ) [ A / c m ] Left boundary − − Right boundary p=6p=7p=8p=9p=10FEMInput . . . . . Time t [ µ s ] . . . N o r m a li z ed fl u x e rr o r E J ( t ) [ A / c m ] . . . . . Time t [ µ s ] . . . p=6p=7p=8p=9p=10Input Fig. 5.11 . Normalized photocurrent due tothe test input (5.15) with an electric field (4.16) applied. The top two plots show the FEM solu-tion (5.3) - (5.4) and DMD solution (5.2) , and thebottom two plots show the DMD training error (5.17) . The left two plots correspond to the fluxout of the left side of the N -region, and similarlyon the right. FEM vs. DMD flux(training input, no electric field) N o r m a li z ed fl u x J ( t ) [ A / c m ] Left boundary − − − Right boundary p=6p=7p=8p=9p=10FEMInput . . . . . Time t [ µ s ] . . N o r m a li z ed fl u x e rr o r E J ( t ) [ A / c m ] . . . . . Time t [ µ s ] . . p=6p=7p=8p=9p=10Input Fig. 5.12 . Normalized photocurrent due tothe training input (5.16) with no electric field ap-plied. The top two plots show the FEM solu-tion (5.3) - (5.4) and DMD solution (5.2) , and thebottom two plots show the DMD training error (5.17) . The left two plots correspond to the fluxout of the left side of the N -region, and similarlyon the right. FEM vs. DMD flux(test input, no electric field) N o r m a li z ed fl u x J ( t ) [ A / c m ] Left boundary − − Right boundary p=6p=7p=8p=9p=10FEMInput . . . . . Time t [ µ s ] . . . N o r m a li z ed fl u x e rr o r E J ( t ) [ A / c m ] . . . . . Time t [ µ s ] . . . p=6p=7p=8p=9p=10Input Fig. 5.13 . Normalized photocurrent due tothe test input (5.15) with no electric field applied.The top two plots show the FEM solution (5.3) - (5.4) and DMD solution (5.2) , and the bottom twoplots show the DMD training error (5.17) . Theleft two plots correspond to the photocurrent outof the left side of the N -region, and similarly onthe right. mates the flow map of the Ambipolar Diffusion Equation. To obtain the approximateflow map we use the Ambipolar Diffusion Equation to reconstruct numerically theinternal state of the device, which is not directly observable through a laboratorymeasurement, and then apply Dynamic Mode Decomposition to the simulated in-ternal state samples. In doing so physics knowledge is incorporated into the model6 Compact Photocurrent Models Using DMD
FEM vs. DMD density(training input, with electric field) . . . . C a rr i e r den s i t y u ( x j , t ) [ / c m ] x = 0 . ℓ x = ℓ x = 1 . ℓ p=6p=7p=8p=9p=10FEMInput . . . . . . . . C a rr i e r den s i t y e rr o r E u ( x j , t ) [ / c m ] . . . . . Time t [ µ s ] . . . . . p=6p=7p=8p=9p=10Input Fig. 5.14 . Simulated excess carrier densitydue to the training input (5.16) with an electricfield (4.16) applied. The top three plots show theFEM solution (3.3) and DMD solution (5.1) , andthe bottom three plots show the DMD training er-ror (5.18) . FEM vs. DMD density(test input, with electric field) . . . C a rr i e r den s i t y u ( x j , t ) [ / c m ] x = 0 . ℓ x = ℓ x = 1 . ℓ p=6p=7p=8p=9p=10FEMInput . . . . . . . . C a rr i e r den s i t y e rr o r E u ( x j , t ) [ / c m ] . . . . . Time t [ µ s ] . . . . . p=6p=7p=8p=9p=10Input Fig. 5.15 . Simulated excess carrier den-sity due to the test input (5.15) with an electricfield (4.16) applied. The top three plots show theFEM solution (3.3) and DMD solution (5.1) , andthe bottom three plots show the DMD test error (5.18) . FEM vs. DMD density(training input, no electric field) . . . . C a rr i e r den s i t y u ( x j , t ) [ / c m ] x = 0 . ℓ x = ℓ x = 1 . ℓ p=6p=7p=8p=9p=10FEMInput . . . . . . . . C a rr i e r den s i t y e rr o r E u ( x j , t ) [ / c m ] . . . . . Time t [ µ s ] . . . . . p=6p=7p=8p=9p=10Input Fig. 5.16 . Simulated excess carrier densitydue to the training input (5.16) with no electricfield applied. The top three plots show the FEMsolution (3.3) and DMD solution (5.1) , and thebottom three plots show the DMD training error (5.18) . FEM vs. DMD density(test input, no electric field) . . . C a rr i e r den s i t y u ( x j , t ) [ / c m ] x = 0 . ℓ x = ℓ x = 1 . ℓ p=6p=7p=8p=9p=10FEMInput . . . . . . . . . C a rr i e r den s i t y e rr o r E u ( x j , t ) [ / c m ] . . . . . Time t [ µ s ] . . . . . p=6p=7p=8p=9p=10Input Fig. 5.17 . Simulated excess carrier densitydue to the test input (5.15) with no electric fieldapplied. The top three plots show the FEM solu-tion (3.3) and DMD solution (5.1) , and the bottomthree plots show the DMD test error (5.18) . development, which allows us to obtain models from sparse data sets that accuratelyapproximate the dynamics of the excess carrier density. This in turn allows us to ac-curately estimate the induced current at the device boundaries, which is the quantityrequired for circuit simulations.Our results confirm that such physics-aware data-driven models are a viable alter-native to traditional analytic compact models that use simplified analytic solutions ofthe governing equations and often must undergo recalibration and/or redevelopmentto include new physics effects. . Hanson, P. Bochev, and B. Paskaleva Acknowledgments.
This work was supported by the Sandia National Laborato-ries (SNL) Laboratory-directed Research and Development (LDRD) program, and theU.S. Department of Energy, Office of Science, Office of Advanced Scientific Comput-ing Research under Award Number DE-SC-0000230927 and under the Collaboratoryon Mathematics and Physics-Informed Learning Machines for Multiscale and Multi-physics Problems (PhILMs) project. The work of J. Hanson was performed as partof Sandia Summer Student Program Internship at the Computer Science ResearchInstitute (CSRI).Sandia National Laboratories is a multimission laboratory managed and operatedby National Technology and Engineering Solutions of Sandia, LLC., a wholly ownedsubsidiary of Honeywell International, Inc., for the U.S. Department of Energys Na-tional Nuclear Security Administration contract number DE-NA0003525. This paperdescribes objective technical results and analysis. Any subjective views or opinionsthat might be expressed in the paper do not necessarily represent the views of theU.S. Department of Energy or the United States Government.This work benefited from numerous discussions and interactions with our col-leagues Eric Keiter, Suzey Gao and Larry Musson who shared their expertise withcompact model development and TCAD simulations. We are grateful for their helpand advice throughout the preparation of this paper.
REFERENCES[1]
C. L. Axness, B. Kerr, and T. F. Wunsch , Analytic light—or radiation—induced pn junctionphotocurrent solutions to the multidimensional ambipolar diffusion equation , Journal ofApplied Physics, 96 (2004), pp. 2646–2655.[2]
Y. Bar-Sinai, S. Hoyer, J. Hickey, and M. P. Brenner , Learning data-driven discretizationsfor partial differential equations , Proceedings of the National Academy of Sciences, 116(2019), pp. 15344–15349.[3]
P. Ciarlet , The Finite Element Method for Elliptic Problems , SIAM Classics in AppliedMathematics, SIAM, Philadelphia, 2002.[4]
E. W. Enlow and D. R. Alexander , Photocurrent modeling of modern microcircuit pn junc-tions , IEEE Transactions on Nuclear Science, 35 (1988), pp. 1467–1474.[5]
T. A. Fjeldly, Y. Deng, M. S. Shur, H. P. Hjalmarson, A. Muyshondt, and T. Ytterdal , Modeling of high-dose-rate transient ionizing radiation effects in bipolar devices , IEEETransactions on Nuclear Science, 48 (2001), pp. 1721–1730.[6]
J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton, and J. Nathan Kutz , On dynamic mode decomposition: Theory and applications , Journal of ComputationalDynamics, 1 (2014), pp. 391–421.[7]
K. He, X. Zhang, S. Ren, and J. Sun , Deep residual learning for image recognition , CoRR,abs/1512.03385 (2015).[8]
B. Kerr, C. L. Axness, J. C. Verley, C. E. Hembree, and E. R. Keiter , A new time-dependent analytic model for radiation-induced photocurrent in finite 1d epitaxial diodes ,Sandia Report SAND2012-2161, Sandia National Laboratories, 2012. Compact Photocurrent Models Using DMD[9]
A. Krizhevsky, I. Sutskever, and G. E. Hinton , Imagenet classification with deep convolu-tional neural networks , Commun. ACM, 60 (2017), pp. 84–90.[10]
Y. LeCun, Y. Bengio, and G. Hinton , Deep learning , Nature, 521 (2015), pp. 436 EP –.[11]
J. L. Proctor, S. L. Brunton, and J. N. Kutz , Dynamic mode decomposition with control ,SIAM Journal on Applied Dynamical Systems, 15 (2016), pp. 142–161.[12]
C. Rackauckas, Y. Ma, J. Martensen, C. Warner, K. Zubov, R. Supekar, D. Skinner,and A. Ramadhan , Universal differential equations for scientific machine learning , 2020.[13]
M. Raissi, P. Perdikaris, and G. Karniadakis , Physics-informed neural networks: A deeplearning framework for solving forward and inverse problems involving nonlinear partialdifferential equations , Journal of Computational Physics, 378 (2019), pp. 686 – 707.[14]
P. J. SCHMID , Dynamic mode decomposition of numerical and experimental data , Journal ofFluid Mechanics, 656 (2010), pp. 5–28.[15]
S. Selberherr , Analysis and simulation of semiconductor devices , Springer-Verlag, Berlin,1984.[16]
L. F. Shampine and M. W. Reichelt , The matlab ode suite , SIAM Journal on ScientificComputing, 18 (1997), pp. 1–22.[17]
J. L. Wirth and S. C. Rogers , The transient response of transistors and diodes to ionizingradiation , IEEE Transactions on Nuclear Science, 11 (1964), pp. 24–38.[18]